**2. The definition and rigorous formulation of interfacial area concentration**

Interfacial area concentration is defined as interfacial area per unit volume of two-phase flow. Therefore, the term "interfacial area concentration" is usually used in the meaning of

Transport of Interfacial Area Concentration in Two-Phase Flow 89

Kataoka et al. (1986), Kataoka (1986) and Morel (2007) derived the local instant formulation

One considers the one dimensional case where only one plane interface exists at the position of x=x0, as shown in Fig.2. In the control volume which encloses the interface in the width of

1

When one takes the limit of x 0 , local interfacial area concentration, ai, is obtained. It takes the value of zero at 0 x x and infinity at x=x0 . This local interfacial area concentration

ai=(x- x0) (4)

This formulation can be easily extended to three dimensional case. As Shown in Fig.3, three

f(x,y,z,t)=0 (5)

f(x,y,z,t)>0 (gas phase), f(x,y,z,t)<0 (liquid phase) (6)

dimensional interface of gas and liquid is mathematically given by

<sup>x</sup> (3)

i

a

x, as shown in Fig.2, average interfacial area concentration is given by

of interfacial area concentration as follows.

is given in term of delta function by

Fig. 2. Plane interface at x=x0

volume averaged value and denoted by <sup>V</sup> <sup>i</sup> a . For example, one considers the interfacial area concentration in bubbly flow as shown in Fig.1**.** In this figure, Ai is instantaneous interfacial area included in volume, V. The volume averaged interfacial area concentration is given by

$$\overline{\mathbf{a}\_i}^{\mathbf{v}} = \frac{\mathbf{A}\_i}{\mathbf{V}} \tag{1}$$

For simplicity, bubbles are sphere of which diameter is db, interfacial area concentration is given by

$$\overline{\mathbf{a}\_i}^{\mathbf{V}} = \frac{\pi \mathbf{N} \mathbf{d}\_{\mathbf{b}}^{\mathbf{Z}}}{\mathbf{V}} = \frac{6a}{\mathbf{d}\_{\mathbf{b}}} \tag{2}$$

Here, N is number of bubbles in volume V, and is void fraction (volumetric fraction of bubbles in volume V).

Fig. 1. Interfacial area in bubbly flow

Similarly, time averaged interfacial area concentration, i a and statistical averaged interfacial area concentration, <sup>A</sup> <sup>i</sup> a can be defined. The transport equation of interfacial area concentration is usually given in time averaged form in terms of time averaged interfacial area concentration, i a . However, for the derivation of the transport equation, it is desirable to formulate interfacial area concentration and its transport equation in local instant form.

concentration in bubbly flow as shown in Fig.1**.** In this figure, Ai is instantaneous interfacial area included in volume, V. The volume averaged interfacial area concentration is given by

> <sup>V</sup> <sup>i</sup> i

For simplicity, bubbles are sphere of which diameter is db, interfacial area concentration is

Here, N is number of bubbles in volume V, and is void fraction (volumetric fraction of

Similarly, time averaged interfacial area concentration, i a and statistical averaged

concentration is usually given in time averaged form in terms of time averaged interfacial area concentration, i a . However, for the derivation of the transport equation, it is desirable to formulate interfacial area concentration and its transport equation in local instant form.

<sup>i</sup> a can be defined. The transport equation of interfacial area

Total Volume of

Two-Phase Flow, V

<sup>2</sup> <sup>V</sup> <sup>b</sup>

Total Surface Area

of Bubbles, Ai

i

a

a

A

Nd 6

V d

b

<sup>i</sup> a . For example, one considers the interfacial area

<sup>V</sup> (1)

(2)

volume averaged value and denoted by <sup>V</sup>

given by

bubbles in volume V).

Fig. 1. Interfacial area in bubbly flow

db

interfacial area concentration, <sup>A</sup>

Kataoka et al. (1986), Kataoka (1986) and Morel (2007) derived the local instant formulation of interfacial area concentration as follows.

One considers the one dimensional case where only one plane interface exists at the position of x=x0, as shown in Fig.2. In the control volume which encloses the interface in the width of x, as shown in Fig.2, average interfacial area concentration is given by

$$\overline{\mathbf{a}\_i} = \frac{1}{\Delta \mathbf{x}}\tag{3}$$

When one takes the limit of x 0 , local interfacial area concentration, ai, is obtained. It takes the value of zero at 0 x x and infinity at x=x0 . This local interfacial area concentration is given in term of delta function by

ai=(x- x0) (4)

$$\begin{array}{c|c} \begin{array}{c} \begin{array}{c} \Delta X\\ \stackrel{\Delta X}{\longrightarrow} \end{array} \end{array} & \begin{array}{c} \begin{array}{c} \Delta X\\ \stackrel{\Delta X}{\longrightarrow} \end{array} \\\\ \times \\\\ \times \\\\ \times \\\\ \hline \\\\ \text{Interface} \end{array} \end{array}$$

Fig. 2. Plane interface at x=x0

This formulation can be easily extended to three dimensional case. As Shown in Fig.3, three dimensional interface of gas and liquid is mathematically given by

$$\mathbf{f}(\mathbf{x}, \mathbf{y}, \mathbf{z}, \mathbf{t}) = \mathbf{0} \tag{5}$$

$$\mathbf{f}(\mathbf{x}, \mathbf{y}, \mathbf{z}, \mathbf{t}) \succeq \mathbf{0} \text{ (gas phase)}, \quad \mathbf{f}(\mathbf{x}, \mathbf{y}, \mathbf{z}, \mathbf{t}) \not\le \mathbf{0} \text{ (liquid phase)}\tag{6}$$

Transport of Interfacial Area Concentration in Two-Phase Flow 91

f/|grad f(x,y,z,,t)|

When one takes the limit of f 0 , local interfacial area concentration, ai, is obtained by

ai=|grad f(x,y,z,t)|(f(x,y,z,t)) (10)

<sup>g</sup> 0 0 (w) (w w )dw g(w )

In relation to local instant interfacial area concentration, characteristic function of each phase

where suffixes G and L denote gas and liquid phase respectively. k is the local instant void fraction of each phase and takes the value of unity when phase k exists and takes the value of zero when phase k doesn't exist. Here, h(w) is Heaviside function which is defined by

h(w) =1 (w>0)

dh(w) (w) dw

Using above equations, the derivatives of characteristic function are related to interfacial

<sup>k</sup> v n a (k G,L) i ki i <sup>t</sup>

Here, **n**ki is unit normal outward vector of phase k as shown in Fig.5 and **v**i is the velocity of

Using above-mentioned relations, it is shown that local instant interfacial area concentration is given in term of correlation function of characteristic function (Kataoka (2008)). As shown in

<sup>i</sup> a grad f(x,y,z,t) / f (9)

(11)

G =h(f(x,y,z,t)) (gas phase) (12a)

=0 (w<0) (13)

(14)

grad n a (k G,L) k ki i (15)

(16)

L =1-h(f(x,y,z,t)) (liquid phase) (12b)

By the differential geometry, the width of the control volume is given by

where (w) is the delta function which is defined by

where g(w) is an arbitrary continuous function..

Heaviside function and the delta function are related by

(denoted byk) is defined by

area concentration as follows.

interface.

Then, average interfacial area concentration in this control volume is given by

Fig. 3. Mathematical representation of three-dimensional interface

As shown in Fig.4, one considers the control volume which encloses the interface by following two surfaces.

$$\mathbf{f}(\mathbf{x}, \mathbf{y}, \mathbf{z}, \mathbf{t}) = \mathbf{\hat{n}} \mathbf{f} / \mathbf{2} \tag{7}$$

$$\mathbf{f}(\mathbf{x}, \mathbf{y}, \mathbf{z}, \mathbf{t}) = \mathbf{\therefore} \mathbf{\therefore} / \mathbf{\?} \tag{8}$$

Fig. 4. Control volume enclosing three-dimensional interface

Gas Phase

f(x,y,z,t)>0

As shown in Fig.4, one considers the control volume which encloses the interface by

Interface f(x,y,z,t)=0

f(x,y,z,t)=f/2 (7)

f(x,y,z,t)=-f/2 (8)

Fig. 3. Mathematical representation of three-dimensional interface

Liquid Phase

f(x,y,z,t)<0

Fig. 4. Control volume enclosing three-dimensional interface

following two surfaces.

By the differential geometry, the width of the control volume is given by

$$\left| \text{Af} \right\rangle \mid \text{grad } \mathbf{f}(\mathbf{x}\_{\prime} \mathbf{y}\_{\prime} \mathbf{z}\_{\prime\prime} \mathbf{t}) \mid \mathbf{z}\_{\prime}$$

Then, average interfacial area concentration in this control volume is given by

$$\overline{\mathbf{a}\_i} = \left| \text{grad } \mathbf{f}(\mathbf{x}, \mathbf{y}, \mathbf{z}, \mathbf{t}) \right| / \Delta \mathbf{f} \tag{9}$$

When one takes the limit of f 0 , local interfacial area concentration, ai, is obtained by

$$\mathbf{a} = \lfloor \text{grad } \mathbf{f}(\mathbf{x}, \mathbf{y}, \mathbf{z}, \mathbf{t}) \rfloor \, \delta(\mathbf{f}(\mathbf{x}, \mathbf{y}, \mathbf{z}, \mathbf{t})) \, \tag{10}$$

where (w) is the delta function which is defined by

$$\int\_{-\infty}^{\infty} \mathbf{g}(\mathbf{w}) \delta(\mathbf{w} - \mathbf{w}\_0) d\mathbf{w} = \mathbf{g}(\mathbf{w}\_0) \tag{11}$$

where g(w) is an arbitrary continuous function..

In relation to local instant interfacial area concentration, characteristic function of each phase (denoted byk) is defined by

$$\phi\_{\mathbb{C}} \triangleq \mathbf{h}(\mathbf{f}(\mathbf{x}, \mathbf{y}, \mathbf{z}, \mathbf{t})) \tag{12a} \\ \text{ (gas phase)} \\ \tag{12a}$$

$$\boldsymbol{\Phi}\_{\text{L}} = \mathbf{1} \text{-} \mathbf{h}(\mathbf{f}(\mathbf{x}, \mathbf{y}, \mathbf{z}, \mathbf{t})) \text{ (liquid phase)}\tag{12b}$$

where suffixes G and L denote gas and liquid phase respectively. k is the local instant void fraction of each phase and takes the value of unity when phase k exists and takes the value of zero when phase k doesn't exist. Here, h(w) is Heaviside function which is defined by

$$\begin{array}{c} \text{h(w)} = 1 \text{ (w > 0)}\\ = 0 \text{ (w < 0)} \end{array} \tag{13}$$

Heaviside function and the delta function are related by

$$\delta(\mathbf{w}) = \frac{\text{dh}(\mathbf{w})}{\text{dw}} \tag{14}$$

Using above equations, the derivatives of characteristic function are related to interfacial area concentration as follows.

$$\mathbf{k}\text{ grad }\phi\_{\mathbf{k}} = -\mathbf{n}\_{\text{ki}}\mathbf{a}\_{\text{i}} \qquad \text{(k} = \mathbf{G}\text{,L)}\tag{15}$$

$$\frac{\partial \phi\_{\mathbf{k}}}{\partial \mathbf{t}} = \mathbf{v}\_{\mathbf{i}} \bullet \mathbf{n}\_{\mathbf{k}i} \mathbf{a}\_{\mathbf{i}} \quad \text{(k} = \mathbf{G}\_{\prime}\text{L)}\tag{16}$$

Here, **n**ki is unit normal outward vector of phase k as shown in Fig.5 and **v**i is the velocity of interface.

Using above-mentioned relations, it is shown that local instant interfacial area concentration is given in term of correlation function of characteristic function (Kataoka (2008)). As shown in

Transport of Interfacial Area Concentration in Two-Phase Flow 93

(x r) (x) 0 (0 /2)

cos a ( / 2 )

k k i <sup>1</sup> (x r) (x) ( cos cos )a r 2

kk k i (x) 2 (x r) (x) cos a

kk k <sup>i</sup> <sup>i</sup> 0 0 0 0 { (x) 2 (x r) (x)}sin d d cos a sin d d 2 a

<sup>a</sup> { (x) 2 (x r) (x)}sin d d 2r r

(23)

k k

k k k kk k k

k k k k

(x r) (x) (x) r r

(x r) (x r) (x r) (x) (x r) (x r) (x) (x r) (x) rr r

(x){ (x r) 1} (x r){ (x) 1} (x) 2 (x r) (x)

( ){1 ( )} ( ){1 ( )} cos *kk k k*

 

(27)

**x xr xr x r**

(26)

i

(21)

(22)

(24)

is directional differentiation of characteristic function k(**x**) in **r**

(25)

(20)

In view of Fig.6, the product of k(**x+r**) and Eq(19) is given by

k k

Equation (19) is rewritten by

From Eqs.(19) and (21), one obtains

Rearranging Eq(23), one obtains

As stated above,

direction.

Integrating Eq.(22) for all **r** directions, one obtains

r r

r 

Eq.(24) in the interval of |**r**|, one obtains,

k k

2 2

2

1

Then, the integrated function in Eq.(24) can be given by

kk k

*i a*

r

r r

<sup>i</sup> kk k 0 0

When one approximates the directional differentiation of characteristic function in

k

rr r

Eq.(15), local instant interfacial concentration is related to the derivative of characteristic function of each phase. Here, directional differentiation of characteristic function is considered. Spatial coordinate (x,y,z) is denoted by vector **x** and displacement vector of **r** is defined.

Fig. 5. Unit normal outward vector of phase k

$$\mathbf{x} = (\mathbf{x}, \mathbf{y}, \mathbf{z})\tag{17}$$

$$\mathbf{r} = (\mathbf{r}\_{\times}, \mathbf{r}\_{\mathbf{y}}, \mathbf{r}\_{\mathbf{z}}) \tag{18}$$

At position **x**, directional differentiation of characteristic function k(**x**) in **r** direction (denoted by r ) is defined by

$$\begin{aligned} \frac{\partial}{\partial \mathbf{r}} \, \Phi\_{\mathbf{k}}(\mathbf{x}) &= \mathbf{n}\_{\mathbf{r}} \bullet \mathbf{grad} \Phi\_{\mathbf{k}}(\mathbf{x}) \\ &= -\mathbf{n}\_{\mathbf{r}} \bullet \mathbf{n}\_{\mathbf{k}i} \mathbf{a}\_{\mathbf{i}} \\ &= -\cos \theta \quad \mathbf{a}\_{\mathbf{i}} \end{aligned} \tag{19}$$

Here, **nr** is unit vector of **r** direction and is the angle between **nr** and **n**ki as shown in Fig. 6.

Fig. 6. Configuration of **nr** and **n**ki

In view of Fig.6, the product of k(**x+r**) and Eq(19) is given by

$$\begin{aligned} \left(\Phi\_{\mathbf{k}}(\mathbf{x} + \mathbf{r}) \frac{\partial}{\partial \mathbf{r}} \,\, \Phi\_{\mathbf{k}}(\mathbf{x}) = 0 & \quad \text{ (} 0 \le \theta \le \pi \text{ / 2)}\\ &= -\cos \theta \,\, \mathbf{a}\_{\mathbf{i}} \quad \text{ (}\pi / \, 2 \le \theta \le \pi\text{)}\end{aligned} \tag{20}$$

Equation (19) is rewritten by

92 Nuclear Reactors

Eq.(15), local instant interfacial concentration is related to the derivative of characteristic function of each phase. Here, directional differentiation of characteristic function is considered.

**n**ki

x=(x,y,z) (17)

Interface

 r=(rx,ry,rz) (18) At position **x**, directional differentiation of characteristic function k(**x**) in **r** direction

kr k

Here, **nr** is unit vector of **r** direction and is the angle between **nr** and **n**ki as shown in Fig. 6.

 n na cos a

 

r

**nr x**

Interface

Phase k

r ki i i

Phase k

**n**ki

k(**x**+**r**)=1

**x**

(19)

Interface

**nr**

(x) n grad (x)

Spatial coordinate (x,y,z) is denoted by vector **x** and displacement vector of **r** is defined.

Fig. 5. Unit normal outward vector of phase k

) is defined by

**n**ki

k(**x**+**r**)=0

(denoted by

r 

Fig. 6. Configuration of **nr** and **n**ki

Phase k

$$\left(\phi\_{\mathbf{k}}(\mathbf{x}+\mathbf{r})\frac{\partial}{\partial \mathbf{r}}\,\phi\_{\mathbf{k}}(\mathbf{x})\right) = \frac{1}{2}(-\cos\theta \,\, + \left|\cos\theta\right|) \mathbf{a}\_{\mathbf{i}} \tag{21}$$

From Eqs.(19) and (21), one obtains

$$\frac{\partial}{\partial \mathbf{r}} \left| \phi\_{\mathbf{k}}(\mathbf{x}) - 2\phi\_{\mathbf{k}}(\mathbf{x} + \mathbf{r}) \frac{\partial}{\partial \mathbf{r}} \left| \phi\_{\mathbf{k}}(\mathbf{x}) \right.\right. \\ = - \left| \cos \theta \right| \left. \mathbf{a}\_{\mathbf{i}} \right. \tag{22}$$

Integrating Eq.(22) for all **r** directions, one obtains

$$\int\_{0}^{2\pi} \int\_{0}^{\pi} \langle \frac{\partial}{\partial \mathbf{r}} \, \Phi\_{\mathbf{k}}(\mathbf{x}) - 2\phi\_{\mathbf{k}}(\mathbf{x} + \mathbf{r}) \frac{\partial}{\partial \mathbf{r}} \, \Phi\_{\mathbf{k}}(\mathbf{x}) \rangle \sin \theta \mathbf{d} \theta d\mathbf{r} \, \eta = \int\_{0}^{2\pi} \int\_{0}^{\pi} - \left| \cos \theta \right| \, \mathbf{a}\_{\mathbf{i}} \sin \theta d\theta d\eta = -2\pi \mathbf{a}\_{\mathbf{i}} \tag{23}$$

Rearranging Eq(23), one obtains

$$\mathbf{a}\_{\mathbf{i}} = -\frac{1}{2\pi} \int\_{0}^{2\pi} \int\_{0}^{\pi} \frac{\partial}{\partial \mathbf{r}} \, \Phi\_{\mathbf{k}}(\mathbf{x}) - 2\phi\_{\mathbf{k}}(\mathbf{x} + \mathbf{r}) \frac{\partial}{\partial \mathbf{r}} \, \Phi\_{\mathbf{k}}(\mathbf{x}) |\sin \theta d\theta d\eta \tag{24}$$

As stated above, r is directional differentiation of characteristic function k(**x**) in **r** direction.

When one approximates the directional differentiation of characteristic function in Eq.(24) in the interval of |**r**|, one obtains,

$$\frac{\partial}{\partial \mathbf{r}} \, \phi\_{\mathbf{k}}(\mathbf{x}) \approx \frac{\phi\_{\mathbf{k}}(\mathbf{x} + \mathbf{r}) - \phi\_{\mathbf{k}}(\mathbf{x})}{|\mathbf{r}|} \tag{25}$$

$$\phi\_{\mathbf{k}}(\mathbf{x}+\mathbf{r}) \, \frac{\partial}{\partial \mathbf{r}} \, \phi\_{\mathbf{k}}(\mathbf{x}) \approx \frac{\phi\_{\mathbf{k}}(\mathbf{x}+\mathbf{r})\phi\_{\mathbf{k}}(\mathbf{x}+\mathbf{r}) - \phi\_{\mathbf{k}}(\mathbf{x}+\mathbf{r})\phi\_{\mathbf{k}}(\mathbf{x})}{|\mathbf{r}|} = \frac{\phi\_{\mathbf{k}}(\mathbf{x}+\mathbf{r}) - \phi\_{\mathbf{k}}(\mathbf{x}+\mathbf{r})\phi\_{\mathbf{k}}(\mathbf{x})}{|\mathbf{r}|} \tag{26}$$

Then, the integrated function in Eq.(24) can be given by

$$
\frac{\partial}{\partial \mathbf{r}} \left. \phi\_{\mathbf{k}}(\mathbf{x}) - 2\phi\_{\mathbf{k}}(\mathbf{x} + \mathbf{r}) \right. \\
\left. \frac{\partial}{\partial \mathbf{r}} \left. \phi\_{\mathbf{k}}(\mathbf{x}) \approx \frac{\phi\_{\mathbf{k}}(\mathbf{x}) \{\phi\_{\mathbf{k}}(\mathbf{x} + \mathbf{r}) - 1\} + \phi\_{\mathbf{k}}(\mathbf{x} + \mathbf{r})\} \phi\_{\mathbf{k}}(\mathbf{x}) - 1\right] \tag{27}
$$

$$
\left( a\_{i} \left| \cos \theta \right| \approx \frac{\phi\_{\mathbf{k}}(\mathbf{x}) \{1 - \phi\_{\mathbf{k}}(\mathbf{x} + \mathbf{r})\} + \phi\_{\mathbf{k}}(\mathbf{x} + \mathbf{r}) \{1 - \phi\_{\mathbf{k}}(\mathbf{x})\}}{\left| \mathbf{r} \right|} \right)
$$

Transport of Interfacial Area Concentration in Two-Phase Flow 95

(denoted byk) given by Eqs.(11) or (12) is needed. Kataoka(1986) also derived local instant formulation of two-phase flow which gives the local instant conservation equations of mass momentum and energy in each phase. The conservation equation of characteristic function

k k k k k ki ki i ki i ( ) div( v ) (v V ) n a (k G,L)

Here, k, **v**k are density, velocity of each phase. Suffix ki is value of phase k at interface. Using Eqs.(24) and (30), Local instant conservation equation of interfacial area concentration

<sup>2</sup> k k

<sup>2</sup> k k

(32)

i i i i Gi Gi

V (a V )/a i ii i (33)

(34)

V a (V n )n / a i i i Gi Gi i (35)

(36)

Averaging Eq.(31), one obtains the conservation equation of averaged interfacial area

The right hand side of Eqs.(31) and(32) represent the source term of interfacial area concentration due to the deformation of interface. In the dispersed flows such as bubbly flow and droplet flow, this term correspond breakup or coalescence of bubbles and droplets. Morel (2007) derived the conservation equation of averaged interfacial area concentration

> <sup>a</sup> a V a (V n ) n <sup>t</sup>

Here, denotes time averaging and Vi is the time averaged velocity of interface which

The research group directed by Prof. Ishii in Purdue university derived the transport equation of interfacial area concentration of time averaged interfacial area concentration based on the transport equation of number density function of bubbles

4

i i j ph j 1

(Kocamustafaogullari and Ishii (1995), Hibiki and Ishii (2000a)). It is given by

a aV t

i

(31)

<sup>i</sup> i i k k <sup>k</sup> 0 0 <sup>1</sup> v v (a ) grad(a V ) { grad( ) 2 (x r) grad( )}sin d d

<sup>i</sup> i i k k <sup>k</sup> 0 0 <sup>1</sup> v v (a ) grad(a V ) { grad( ) 2 (x r) grad( )}sin d d

t 2r r 

t 2r r 

where averaged interfacial velocity Vi is defined by

based on the detailed geometrical consideration of interface.

i

(30)

of each phase (denoted byk) given by

t

is given by

concentration by

is given by

Using this relation , Eq.(24) can be rewritten by

$$\mathbf{a}\_{i} = \frac{1}{2\pi} \int\_{0}^{2\pi} \int\_{0}^{\pi} \mathbf{l} \lim\_{|\mathbf{r}| \to 0} \frac{\phi\_{\mathbf{k}}(\mathbf{x}) [1 - \phi\_{\mathbf{k}}(\mathbf{x} + \mathbf{r})] + \phi\_{\mathbf{k}}(\mathbf{x} + \mathbf{r}) [1 - \phi\_{\mathbf{k}}(\mathbf{x})]}{|\mathbf{r}|} \quad \text{(\sin\theta d\theta d\eta} \tag{28}$$

Averaging Eq.(28), one obtains,

$$\overline{\mathbf{a}\_{i}} = \frac{1}{2\pi} \int\_{0}^{2\pi} \int\_{0}^{\pi} \lim\_{|\mathbf{r}| \to 0} \overline{\phi\_{\mathbf{k}}(\mathbf{x}) \{1 - \phi\_{\mathbf{k}}(\mathbf{x} + \mathbf{r})\}} + \overline{\phi\_{\mathbf{k}}(\mathbf{x} + \mathbf{r}) \{1 - \phi\_{\mathbf{k}}(\mathbf{x})\}} \quad \text{(\sin\theta d\theta d\eta} \tag{29}$$

In the right hand side of Eq.(29), the term

$$\overline{\phi\_{\mathbf{k}}(\mathbf{x})\{1-\phi\_{\mathbf{k}}(\mathbf{x}+\mathbf{r})\}}+\overline{\phi\_{\mathbf{k}}(\mathbf{x}+\mathbf{r})\{1-\phi\_{\mathbf{k}}(\mathbf{x})\}}$$

represents the probability where gas-liquid interface exists between x and x+r as shown in Fig.7.

Fig. 7. The case where interface exists between **x** and **x**+**r**.

#### **3. Basic transport equations of interfacial area concentration**

Based on the rigorous formulation of interfacial area concentration, one can derive transport equation of interfacial area concentration. The transport equations of interfacial area concentration consist of two equations. One is the conservation equation of interfacial area concentration and the other is the conservation equation of interfacial velocity (velocity of interface), **V**i.

Kataoka (2008) derived the local instant conservation equation of interfacial area concentration based on the formulation given by Eq.(24). In order to obtain the local instant conservation equation of interfacial area concentration, characteristic function of each phase

<sup>1</sup> (x){1 (x r)} (x r){1 (x)} a l im { }sin d d 2 r

<sup>1</sup> (x){1 (x r)} (x r){1 (x)} a lim{ }sin d d 2 r

kk k k (x){1 (x r)} (x r){1 (x)}

represents the probability where gas-liquid interface exists between x and x+r as shown in

(28)

(29)

<sup>2</sup> kk k k

<sup>2</sup> kk k k

Using this relation , Eq.(24) can be rewritten by

Fig. 7. The case where interface exists between **x** and **x**+**r**.

k(**x+r**)=0

**3. Basic transport equations of interfacial area concentration** 

Interface

Based on the rigorous formulation of interfacial area concentration, one can derive transport equation of interfacial area concentration. The transport equations of interfacial area concentration consist of two equations. One is the conservation equation of interfacial area concentration and the other is the conservation equation of interfacial velocity (velocity of

Interface

Phase k **x**

k(**x+r**)=1

k(**x**)=0 **x+r**

Kataoka (2008) derived the local instant conservation equation of interfacial area concentration based on the formulation given by Eq.(24). In order to obtain the local instant conservation equation of interfacial area concentration, characteristic function of each phase

<sup>i</sup> 0 0 r 0

<sup>i</sup> 0 0 r 0

In the right hand side of Eq.(29), the term

Averaging Eq.(28), one obtains,

Phase k

**x** 

k(**x**)=1 **x+r**

Fig.7.

interface), **V**i.

(denoted byk) given by Eqs.(11) or (12) is needed. Kataoka(1986) also derived local instant formulation of two-phase flow which gives the local instant conservation equations of mass momentum and energy in each phase. The conservation equation of characteristic function of each phase (denoted byk) given by

$$\frac{\partial}{\partial t}(\phi\_{\mathbf{k}}\rho\_{\mathbf{k}}) + \text{div}(\phi\_{\mathbf{k}}\rho\_{\mathbf{k}}\mathbf{v}\_{\mathbf{k}}) = -\rho\_{\text{ki}}(\mathbf{v}\_{\mathbf{ki}} - \mathbf{V}\_{\mathbf{i}}) \bullet \mathbf{n}\_{\text{ki}}\mathbf{a}\_{\mathbf{i}} \quad \text{( $\mathbf{k} = \mathbf{G}$ , L)}\tag{30}$$

Here, k, **v**k are density, velocity of each phase. Suffix ki is value of phase k at interface. Using Eqs.(24) and (30), Local instant conservation equation of interfacial area concentration is given by

$$\frac{\partial}{\partial \mathbf{\hat{t}}} (\mathbf{a\_i}) + \text{grad} (\mathbf{a\_i} \mathbf{V\_i}) = \frac{1}{2\pi} \int\_0^{2\pi} \int\_0^{\pi} \langle \frac{\partial \mathbf{v\_k}}{\partial \mathbf{r}} \text{grad} (\mathbf{q\_k}) - 2\phi\_\mathbf{k} (\mathbf{x} + \mathbf{r}) \frac{\partial \mathbf{v\_k}}{\partial \mathbf{r}} \text{grad} (\mathbf{q\_k}) \rangle \sin \theta d\theta d\eta \tag{31}$$

Averaging Eq.(31), one obtains the conservation equation of averaged interfacial area concentration by

$$\frac{\partial}{\partial \mathbf{\hat{t}}} (\overline{\mathbf{a}\_{\mathbf{i}}}) + \operatorname{grad}(\overline{\mathbf{a}\_{\mathbf{i}} \mathbf{V\_{i}}}) = \frac{1}{2\pi} \int\_{0}^{2\pi} \int\_{0}^{\pi} \overline{\langle \frac{\partial \mathbf{v\_{k}}}{\partial \mathbf{r}} \operatorname{grad}(\phi\_{\mathbf{k}}) - 2\phi\_{\mathbf{k}}(\mathbf{x} + \mathbf{r}) \frac{\partial \mathbf{v\_{k}}}{\partial \mathbf{r}} \operatorname{grad}(\phi\_{\mathbf{k}}) \rangle \sin \theta d\theta d\eta} \tag{32}$$

where averaged interfacial velocity Vi is defined by

$$
\overline{\mathbf{V\_i}} = \left(\overline{\mathbf{a\_i V\_i}}\right) / \overline{\mathbf{a\_i}}\tag{33}
$$

The right hand side of Eqs.(31) and(32) represent the source term of interfacial area concentration due to the deformation of interface. In the dispersed flows such as bubbly flow and droplet flow, this term correspond breakup or coalescence of bubbles and droplets.

Morel (2007) derived the conservation equation of averaged interfacial area concentration based on the detailed geometrical consideration of interface.

$$\frac{\partial \overline{\mathbf{a}\_i}}{\partial \mathbf{t}} + \nabla \bullet \overline{\mathbf{a}\_i \overline{\mathbf{V}\_i}} = \overline{\mathbf{a}\_i (\mathbf{V}\_i \bullet \mathbf{n}\_{\rm Ci}) \nabla \bullet \mathbf{n}\_{\rm Ci}}\tag{34}$$

Here, denotes time averaging and Vi is the time averaged velocity of interface which is given by

$$\overline{\mathbf{V\_i}} = \overline{\mathbf{a\_i}(\mathbf{V\_i} \bullet \mathbf{n\_{Gi}}) \mathbf{n\_{Gi}}} / \overline{\mathbf{a\_i}} \tag{35}$$

The research group directed by Prof. Ishii in Purdue university derived the transport equation of interfacial area concentration of time averaged interfacial area concentration based on the transport equation of number density function of bubbles (Kocamustafaogullari and Ishii (1995), Hibiki and Ishii (2000a)). It is given by

$$\frac{\partial \mathbf{a}\_i}{\partial \mathbf{t}} + \nabla \bullet \overline{\mathbf{a}\_i \mathbf{V}\_i} = \sum\_{j=1}^4 \phi\_j + \phi\_{\text{ph}} \tag{36}$$

Transport of Interfacial Area Concentration in Two-Phase Flow 97

On the other hand, using Eq.(29), following relation can be obtained for averaged velocity of

<sup>2</sup> k k k kk k

2 r

<sup>1</sup> (x){1 (x r)}v (x) {1 (x)} (x r)v (x) v a lim sin d d

Using, Eqs.(42) and (43), the difference between time averaged interfacial velocity, i V and

<sup>1</sup> (x){1 (x r)}{v (x) v (x)} {1 (x)} (x r){v (x r) v (x)} lim sin d d

k k k k kk k k

are fluctuating terms of local instant volume fraction and velocity of

(47)

(48)

(x){1 (x r)}{v (x) v (x)} {1 (x)} (x r){v (x r) v (x)} lim r

v vv k kk

k kk

Equations (45) and (46) indicate that the difference between time averaged interfacial velocity, i V and time averaged velocity of phase k, vk is given in terms of correlations between fluctuating terms of local instant volume fraction and velocity of phase k which are

Then, it is important to derive the governing equation of the correlation term given by Eq.(46). In what follows, one derives the governing equation based on the local instant basic equations of mass conservation and momentum conservation of phase k which are given below (Kataoka (1986)). In these conservation equations, tensor representation is used. Einstein abbreviation rule is also applied. When the same suffix appear, summation for that

suffix is carried out except for the suffix k denoting gas and liquid phases.

<sup>2</sup> k k k k kk k k

Rearranging the term in integration in the right hand side of Eq.(45) one obtains

k k k kk k

r

(x) (x r)v (x r) (x) (x r)v (x) lim

(43)

v ( v )/ k kk k (44)

(45)

(46)

each phase, k v by

where k v is defined by

i ki

r 0

Here, k

and vk

phase k which are given by

(V v )a

0 0 r 0

r 0

related to turbulence terms of phase k.

k i 0 0 r 0

time average velocity of phase k, vk is given by

2 r

Here, the first term in right hand side of Eq.(36) represent the source and sink terms due to bubble coalescence and break up. Interfacial area decreases when bubbles coalescence and increases when bubbles break up. This term is quite important in interfacial area transport. Therefore, the constitutive equations of this term are given by Hibiki and Ishii (2000a, 2000b) and Ishii and Kim (2004) based on detailed mechanistic modeling. The second term in right hand side of Eq.(36) represent the source and sink terms due to phase change. Equation (36) is practical transport equation of interfacial area concentration.

As for the conservation equation of interfacial velocity, Kataoka et al. (2010,2011a) have derived rigorous formulation based on the local instant formulation of interfacial area concentration and interfacial velocity, which is shown below. Since interface has no mass, momentum equation of interface cannot be formulated. Therefore, the conservation equation of interfacial velocity (or governing equation of interfacial velocity) has to be derived in collaboration with the momentum equation of each phase. Since interfacial velocity is only defined at interface, local instant formulation of interfacial velocity must be expressed in the form of

#### i i a V

Using Eq.(22), interfacial velocity is expressed by

$$\left| \mathbf{a}\_{\mathrm{i}} \mathbf{V}\_{\mathrm{i}} \right| \cos \theta \left| \begin{array}{c} \cos \theta \\ \end{array} \right| = \left| \begin{array}{c} \text{V}\_{\mathrm{i}} \\ \text{\raisebox{0.0pt}{0}} \Phi\_{\mathrm{k}} (\mathbf{x}) + 2 \,\mathrm{V}\_{\mathrm{i}} \Phi\_{\mathrm{k}} (\mathbf{x} + \mathbf{r}) \frac{\partial}{\partial \mathbf{r}} \,\Phi\_{\mathrm{k}} (\mathbf{x}) \end{array} \right. \tag{37}$$

Considering Fig.7, interfacial velocity is approximated by velocity of each phase without phase change.

$$\mathbf{V}\_{\mathbf{i}} \approx \phi\_{\mathbf{k}}(\mathbf{x}) \{ 1 - \phi\_{\mathbf{k}}(\mathbf{x} + \mathbf{r}) \} \mathbf{v}\_{\mathbf{k}}(\mathbf{x}) + \phi\_{\mathbf{k}}(\mathbf{x} + \mathbf{r}) \{ 1 - \phi\_{\mathbf{k}}(\mathbf{x}) \} \mathbf{v}\_{\mathbf{k}}(\mathbf{x} + \mathbf{r}) \tag{38}$$

From Eqs.(27) and (38) with some rearrangements, one obtains following approximate expression.

$$\left| \mathbf{a}\_{i} \mathbf{V}\_{i} \middle| \cos \theta \right| \approx \frac{\phi\_{\mathbf{k}}(\mathbf{x}) \{1 - \phi\_{\mathbf{k}}(\mathbf{x} + \mathbf{r})\} \mathbf{v}\_{\mathbf{k}}(\mathbf{x}) + \phi\_{\mathbf{k}}(\mathbf{x} + \mathbf{r}) \{1 - \phi\_{\mathbf{k}}(\mathbf{x})\} \mathbf{v}\_{\mathbf{k}}(\mathbf{x} + \mathbf{r})}{|\mathbf{r}|} \tag{39}$$

When one takes the limit of r 0 , one obtains

$$\mathbf{a}\_{\mathbf{i}} \mathbf{V}\_{\mathbf{i}} \| \cos \theta \| = \lim\_{|\mathbf{r}| \to 0} \frac{\phi\_{\mathbf{k}}(\mathbf{x}) \{1 - \phi\_{\mathbf{k}}(\mathbf{x} + \mathbf{r})\} \mathbf{v}\_{\mathbf{k}}(\mathbf{x}) + \phi\_{\mathbf{k}}(\mathbf{x} + \mathbf{r}) \{1 - \phi\_{\mathbf{k}}(\mathbf{x})\} \mathbf{v}\_{\mathbf{k}}(\mathbf{x} + \mathbf{r})}{|\mathbf{r}|} \tag{40}$$

Integrating Eq.(40) in all direction, one finally obtains

$$\mathbf{a}\_{i}\mathbf{V}\_{i} = \frac{1}{2\pi} \int\_{0}^{2\pi} \int\_{0}^{\pi} \lim\_{\|\mathbf{r}\|\to 0} \frac{\phi\_{\mathbf{k}}(\mathbf{x}) \{1 - \phi\_{\mathbf{k}}(\mathbf{x} + \mathbf{r})\} \mathbf{v}\_{\mathbf{k}}(\mathbf{x}) + \{1 - \phi\_{\mathbf{k}}(\mathbf{x})\} \phi\_{\mathbf{k}}(\mathbf{x} + \mathbf{r}) \mathbf{v}\_{\mathbf{k}}(\mathbf{x} + \mathbf{r})}{|\mathbf{r}|} \sin\theta d\theta d\phi \quad (41)$$

Averaging Eq.(41), averaged interfacial velocity is given by

$$\overline{\mathbf{V}\_{i}\mathbf{a}\_{i}} = \frac{1}{2\pi} \int\_{0}^{2\pi} \int\_{0}^{\pi} \lim\_{\|\mathbf{r}\|\to 0} \frac{\overline{\phi\_{\mathbf{k}}(\mathbf{x}) \{1 - \phi\_{\mathbf{k}}(\mathbf{x} + \mathbf{r})\} \mathbf{v}\_{\mathbf{k}}(\mathbf{x})}{|\mathbf{r}|} + \frac{\overline{\{1 - \phi\_{\mathbf{k}}(\mathbf{x})\} \phi\_{\mathbf{k}}(\mathbf{x} + \mathbf{r}) \mathbf{v}\_{\mathbf{k}}(\mathbf{x} + \mathbf{r})}{\text{sinc}\,\theta d\theta d\phi} \sin\theta d\theta d\phi \tag{42}$$

On the other hand, using Eq.(29), following relation can be obtained for averaged velocity of each phase, k v by

$$\overline{\mathbf{v}\_{\mathbf{k}}}\overline{\mathbf{a}\_{\mathbf{i}}} = \frac{1}{2\pi} \int\_{0}^{2\pi} \int\_{0}^{\pi} \lim\_{\|\mathbf{r}\| \to 0} \overline{\frac{\phi\_{\mathbf{k}}(\mathbf{x}) \langle 1 - \phi\_{\mathbf{k}}(\mathbf{x} + \mathbf{r}) \rangle}{|\mathbf{r}|} \overline{\mathbf{v}\_{\mathbf{k}}(\mathbf{x})}} + \overline{\frac{\langle 1 - \phi\_{\mathbf{k}}(\mathbf{x}) \rangle \phi\_{\mathbf{k}}(\mathbf{x} + \mathbf{r})}{|\mathbf{r}|} \overline{\mathbf{v}\_{\mathbf{k}}(\mathbf{x})}} \sin \theta d\theta d\phi \tag{43}$$

where k v is defined by

96 Nuclear Reactors

Here, the first term in right hand side of Eq.(36) represent the source and sink terms due to bubble coalescence and break up. Interfacial area decreases when bubbles coalescence and increases when bubbles break up. This term is quite important in interfacial area transport. Therefore, the constitutive equations of this term are given by Hibiki and Ishii (2000a, 2000b) and Ishii and Kim (2004) based on detailed mechanistic modeling. The second term in right hand side of Eq.(36) represent the source and sink terms due to phase change. Equation (36)

As for the conservation equation of interfacial velocity, Kataoka et al. (2010,2011a) have derived rigorous formulation based on the local instant formulation of interfacial area concentration and interfacial velocity, which is shown below. Since interface has no mass, momentum equation of interface cannot be formulated. Therefore, the conservation equation of interfacial velocity (or governing equation of interfacial velocity) has to be derived in collaboration with the momentum equation of each phase. Since interfacial velocity is only defined at interface, local instant formulation of interfacial velocity must be

i i a V

i i i k ik <sup>k</sup> a V cos V (x) 2V (x r) (x)

Considering Fig.7, interfacial velocity is approximated by velocity of each phase without

From Eqs.(27) and (38) with some rearrangements, one obtains following approximate

(x){1 (x r)}v (x) (x r){1 (x)}v (x r) a V cos

k k k k kk i i r 0 (x){1 (x r)}v (x) (x r){1 (x)}v (x r) a V cos lim <sup>r</sup>

<sup>2</sup> k k k kk k

<sup>2</sup> k k k kk k

<sup>1</sup> (x){1 (x r)}v (x) {1 (x)} (x r)v (x r) V a lim sin d d

<sup>1</sup> (x){1 (x r)}v (x) {1 (x)} (x r)v (x r) a V lim sin d d

(41)

(42)

k k k k kk

r

r r

V (x){1 (x r)}v (x) (x r){1 (x)}v (x r) ik k k k k k (38)

(39)

(40)

(37)

is practical transport equation of interfacial area concentration.

Using Eq.(22), interfacial velocity is expressed by

expressed in the form of

phase change.

expression.

i i

i i 0 0 r 0

i i 0 0 r 0

When one takes the limit of r 0 , one obtains

Integrating Eq.(40) in all direction, one finally obtains

2 r

2 r

Averaging Eq.(41), averaged interfacial velocity is given by

$$
\overline{\mathbf{v\_k}} = \left(\overline{\phi\_\mathbf{k} \mathbf{v\_k}}\right) / \overline{\phi\_\mathbf{k}}\tag{44}
$$

Using, Eqs.(42) and (43), the difference between time averaged interfacial velocity, i V and time average velocity of phase k, vk is given by

$$\begin{split} \overline{\mathbf{u}\_{i}(\overline{\mathbf{v\_{i}}}-\overline{\mathbf{v\_{k}}})} & \overline{\mathbf{a}\_{i}} \\ = \frac{1}{2\pi} \int\_{0}^{2\pi} \int\_{0}^{\pi} \lim\_{|\mathbf{r}| \to 0} \frac{\phi\_{\mathbf{k}}(\mathbf{x})[1-\phi\_{\mathbf{k}}(\mathbf{x}+\mathbf{r})] \langle \mathbf{v\_{k}}(\mathbf{x}) - \overline{\mathbf{v\_{k}}(\mathbf{x})} \rangle}{|\mathbf{r}|} + \frac{\langle 1-\phi\_{\mathbf{k}}(\mathbf{x}) \rangle \phi\_{\mathbf{k}}(\mathbf{x}+\mathbf{r}) \langle \mathbf{v\_{k}}(\mathbf{x}+\mathbf{r}) - \overline{\mathbf{v\_{k}}(\mathbf{x})} \rangle}{|\mathbf{r}|} \sin\theta d\theta d\phi \end{split} (45)$$

Rearranging the term in integration in the right hand side of Eq.(45) one obtains

$$\begin{split} \min\_{\begin{subarray}{c} \|\mathbf{r}\|\rightarrow 0 \\ \|\mathbf{r}\|\rightarrow 0 \end{subarray}} & \frac{\overline{\phi\_{\mathbf{k}}(\mathbf{x}) [1-\phi\_{\mathbf{k}}(\mathbf{x}+\mathbf{r})] [\mathbf{v}\_{\mathbf{k}}(\mathbf{x})-\overline{\mathbf{v}\_{\mathbf{k}}(\mathbf{x})}]} + \overline{\{1-\phi\_{\mathbf{k}}(\mathbf{x})\}\phi\_{\mathbf{k}}(\mathbf{x}+\mathbf{r}) [\mathbf{v}\_{\mathbf{k}}(\mathbf{x}+\mathbf{r})-\overline{\mathbf{v}\_{\mathbf{k}}(\mathbf{x})}]}}{|\mathbf{r}|} \\ &= -\lim\_{|\mathbf{r}|\rightarrow 0} \frac{\overline{\phi\_{\mathbf{k}}^{'}(\mathbf{x})\phi\_{\mathbf{k}}(\mathbf{x}+\mathbf{r})\mathbf{v}\_{\mathbf{k}}^{'}(\mathbf{x}+\mathbf{r})} + \overline{\phi\_{\mathbf{k}}(\mathbf{x})\phi\_{\mathbf{k}}^{'}(\mathbf{x}+\mathbf{r})\mathbf{v}\_{\mathbf{k}}^{'}(\mathbf{x})}}{|\mathbf{r}|} \end{split} \tag{46}$$

Here, k and vk are fluctuating terms of local instant volume fraction and velocity of phase k which are given by

$$\mathbf{v\_k}' = \mathbf{v\_k} - \overline{\mathbf{v\_k}}\tag{47}$$

$$
\boldsymbol{\phi\_{k}}' = \boldsymbol{\phi\_{k}} - \overline{\boldsymbol{\phi\_{k}}} \tag{48}
$$

Equations (45) and (46) indicate that the difference between time averaged interfacial velocity, i V and time averaged velocity of phase k, vk is given in terms of correlations between fluctuating terms of local instant volume fraction and velocity of phase k which are related to turbulence terms of phase k.

Then, it is important to derive the governing equation of the correlation term given by Eq.(46). In what follows, one derives the governing equation based on the local instant basic equations of mass conservation and momentum conservation of phase k which are given below (Kataoka (1986)). In these conservation equations, tensor representation is used. Einstein abbreviation rule is also applied. When the same suffix appear, summation for that suffix is carried out except for the suffix k denoting gas and liquid phases.

Transport of Interfacial Area Concentration in Two-Phase Flow 99

( v v )}sin d d

<sup>1</sup> ( F F )sin d d

1 11 1 lim ( P n a P n a )sin d d 2 r

1 11 1 lim ( n a n a )sin d d 2 r

1 1 lim { ( v v ) ( v v )}sin d d 2 rx <sup>x</sup>

k r k r k i ir k k ki i <sup>0</sup> kr <sup>k</sup>

<sup>1</sup> ( v v n a v v n a )sin d d

(55)

 

i k i i k ik

{ V v ) a } ( V v ) a v )

k kr k kr 0 0 r 0 <sup>k</sup>

2

t x

2 0 0 r

2

2

 

<sup>1</sup> lim

2 0 0 r

<sup>1</sup> lim

2

2

2

<sup>2</sup> kr <sup>k</sup>

11 1 P P lim ( )sin d d 2 rx x

ki k i ir ki k i i 0 0 r 0 kr k k <sup>k</sup>

ki k i ir ki k i i 0 0 r 0 kr k k <sup>k</sup>

k kr k r k r k kr kr k 0 0 r 0

lim ( v v )sin d d 2rx x

2 k r k k kr k r k kr k 0 0 r 0

 

 

> {v (v v ) }sin d d x

v ( v v )}sin d d

As shown above, the formulation of governing equation of interfacial velocity is derived. Then, the most strict formulation of transport equations of interfacial area concentration is given by conservation equation of interfacial area concentration (Eq.(32) , Eq(34), or Eq.(36)) and conservation equation of interfacial velocity (Eq.(55)). As shown in Eq.(55), the conservation equation of interfacial velocity consists of various correlation terms of fluctuating terms of velocity and local instant volume fractions. These correlation terms represent the turbulent transport of interfacial area, which reflects the interactions between gas liquid interface and turbulence of gas and liquid phases. Equation (55) represents such turbulence transport terms of interfacial area concentration. Accurate predictions of interfacial area transport can be possible by solving the transport equations derived here. However, Eq.(55)

 

k kr k r k k r k r 0 0 r 0 <sup>k</sup>

<sup>2</sup> k kr k kr

<sup>2</sup> k kr k kr

k kr k kr

1 1 v v

k r k r k r kr k 0 0 r 0

1 1 lim {v (v v ) 2r x

k k k k kr

k r kr r k k k r 0 0 r 0

1 1 lim {v ( v v ) 2 rx

k k kr k r kr k

x

11 1 lim { ( v v ) 2 rx

k kr k k k k

k kr k r k kr k <sup>0</sup>

x

r

r

(Mass conservation)

$$
\boldsymbol{\Phi}\_{\mathbf{k}} \frac{\partial \mathbf{v}\_{\mathbf{k}\boldsymbol{\beta}}}{\partial \mathbf{x}\_{\boldsymbol{\beta}}} = \mathbf{0} \tag{49}
$$

(Momentum conservation)

$$
\boldsymbol{\Phi}\_{\mathbf{k}} \frac{\partial \mathbf{v}\_{\mathbf{k}\alpha}}{\partial \mathbf{t}} + \boldsymbol{\Phi}\_{\mathbf{k}} \frac{\partial}{\partial \mathbf{x}\_{\boldsymbol{\beta}}} (\mathbf{v}\_{\mathbf{k}\alpha} \mathbf{v}\_{\mathbf{k}\boldsymbol{\beta}}) = -\boldsymbol{\Phi}\_{\mathbf{k}} \frac{1}{\rho\_{\mathbf{k}}} \frac{\partial \mathbf{P}\_{\mathbf{k}}}{\partial \mathbf{x}\_{\alpha}} + \boldsymbol{\Phi}\_{\mathbf{k}} \frac{1}{\rho\_{\mathbf{k}}} \frac{\partial \mathbf{\tau}\_{\mathbf{k}\alpha\boldsymbol{\beta}}}{\partial \mathbf{x}\_{\boldsymbol{\beta}}} + \boldsymbol{\Phi}\_{\mathbf{k}} \mathbf{F}\_{\mathbf{k}\alpha} \tag{50}
$$

Averaging Eqs.(49) and (50) , one obtains time averaged conservation equation of mass and momentum conservation of phase k.

(Time averaged mass conservation)

$$\frac{\partial \overline{\mathbf{v}\_{\mathbf{k}\boldsymbol{\beta}}}}{\partial \mathbf{x}\_{\boldsymbol{\beta}}} = -\frac{1}{\overline{\Phi\_{\mathbf{k}}}} \overline{\mathbf{v}\_{\mathbf{k}\boldsymbol{\beta}i}' \mathbf{n}\_{\mathbf{k}\boldsymbol{\beta}i} \mathbf{a}\_{i}} \tag{51}$$

(Time averaged momentum conservation)

$$\begin{split} \frac{\partial \overline{\text{V}\_{\text{kox}}}}{\partial \mathbf{t}} + \frac{\partial}{\partial \mathbf{x}\_{\text{\beta}}} \overline{\text{(} \overline{\text{v}\_{\text{kox}} \text{v}\_{\text{k\beta}})}} &= -\frac{1}{\rho\_{\text{k}}} \frac{\partial \overline{\text{P}\_{\text{k}}}}{\partial \mathbf{x}\_{\text{\alpha}}} + \frac{1}{\rho\_{\text{k}}} \frac{\partial}{\partial \mathbf{x}\_{\text{\beta}}} \overline{\left( \overline{\text{r}\_{\text{k}\alpha\beta}} - \rho\_{\text{k}} \overline{\text{v'}\_{\text{k\alpha}} \mathbf{v}\_{\text{k}\beta}'} \right)} + \overline{\text{F}\_{\text{k\alpha}}} - \frac{\overline{\text{v}\_{\text{k\alpha}}}}{\overline{\Phi}\_{\text{k}}} \overline{\mathbf{v'}\_{\text{k\beta}} \mathbf{n}\_{\text{k\beta}} \mathbf{a}\_{\text{i}}} \\ &- \frac{1}{\overline{\Phi}\_{\text{k}}} \frac{1}{\rho\_{\text{k}}} \overline{\mathbf{P}'\_{\text{k\alpha}} \mathbf{n}\_{\text{k\alpha}} \mathbf{a}\_{\text{i}}} + \frac{1}{\overline{\Phi}\_{\text{k}}} \frac{1}{\rho\_{\text{k}}} \overline{\mathbf{r}'\_{\text{k\alpha}} \mathbf{n}\_{\text{k\beta}} \mathbf{a}\_{\text{i}}} + \frac{1}{\overline{\Phi}\_{\text{k}}} \overline{\mathbf{v'}\_{\text{k\alpha}} \mathbf{v}'\_{\text{k\beta}} \mathbf{n}\_{\text{k\beta}} \mathbf{a}\_{\text{i}}} \end{split} \tag{52}$$

Subtracting Eqs(51) and (52) from Eqs.(49) and (50), the conservation equations of fluctuating terms are obtained.

(Conservation equation of mass fluctuation)

$$
\Phi\_{\mathbf{k}} \frac{\partial \mathbf{v}\_{\mathbf{k}\boldsymbol{\beta}}^{\prime}}{\partial \mathbf{x}\_{\boldsymbol{\beta}}} = \frac{\Phi\_{\mathbf{k}}}{\overline{\Phi\_{\mathbf{k}}}} \overline{\mathbf{v}\_{\mathbf{k}\boldsymbol{\beta}\mathbf{i}}^{\prime} \mathbf{n}\_{\mathbf{k}\boldsymbol{\beta}\mathbf{i}} \mathbf{a}\_{\mathbf{i}}} \tag{53}
$$

(Conservation equation of momentum fluctuation)

<sup>k</sup> <sup>k</sup> k k kk k k k <sup>k</sup> k k kk k <sup>k</sup> k k k kk k k k k ki ki i ki k i i k i ki i k k ki i k kk k k k v P 1 1 (v v v v v v ) ( vv) t x x x 1 1 F v v n a Pn a n a v v n a (54)

Using Eqs(53) and (54), one can derive conservation equation of

$$\overline{\phi\_{\mathbf{k}}\acute{\mathbf{(\infty)}}\phi\_{\mathbf{k}}(\mathbf{x}+\mathbf{r})\mathbf{v}\_{\mathbf{k}}\acute{\mathbf{(\infty}+\mathbf{r})} + \overline{\phi\_{\mathbf{k}}(\mathbf{x})\phi\_{\mathbf{k}}\acute{\mathbf{(\infty)}}\mathbf{(\infty)}\mathbf{v}\_{\mathbf{k}}\acute{\mathbf{(\infty)}}}}$$

Then, conservation equation of the difference between interfacial velocity and averaged velocity of each phase is derived. The result is given by

k k v

x 

k k k k k kk k k k k

k

x 

v 1

k k <sup>k</sup>

11 11 1 Pn a n a v v n a

v P 1 1 (v v ) <sup>F</sup> t x x x 

Averaging Eqs.(49) and (50) , one obtains time averaged conservation equation of mass and

k

kk k k k

Subtracting Eqs(51) and (52) from Eqs.(49) and (50), the conservation equations of

k k k ki ki i k

<sup>k</sup> <sup>k</sup> k k kk k k k <sup>k</sup> k k kk k <sup>k</sup> k k

1 1 F v v n a Pn a n a v v n a

k kk k k k k ki ki i ki k i i k i ki i k k ki i k kk k k k

> k k k kk k (x) (x r)v (x r) (x) (x r)v (x)

Then, conservation equation of the difference between interfacial velocity and averaged

v P 1 1 (v v v v v v ) ( vv) t x x x

 

v

x 

vP v 1 1 (v v ) vv F v n a

0

k k

ki ki i

k k k kk k k ki ki i k k k

vna

ki k i i k i ki i k k ki i

vna

(50)

(49)

(51)

(53)

(52)

(54)

(Mass conservation)

(Momentum conservation)

momentum conservation of phase k. (Time averaged mass conservation)

(Time averaged momentum conservation)

(Conservation equation of mass fluctuation)

(Conservation equation of momentum fluctuation)

Using Eqs(53) and (54), one can derive conservation equation of

velocity of each phase is derived. The result is given by

fluctuating terms are obtained.

t x x x

$$\begin{aligned} &\frac{1}{\delta}\Delta\left(\overline{V}\_{\mathbf{k}}-\overline{V}\_{\mathbf{k}}\right)\Big|\_{\delta\mathbf{k}} + \frac{\delta}{\delta\mathbf{p}}\Big{(}\overline{V}\_{\mathbf{k}}-\overline{V}\_{\mathbf{k}}\Big{)}\Big{(}\overline{V}\_{\mathbf{k}}\overline{V}\_{\mathbf{k}}\Big{)} \\ &= \frac{1}{\delta\mathbf{p}}\frac{1}{2\pi}\int\_{0}^{\delta\mathbf{p}}\frac{1}{\delta\mathbf{p}}\Big{(}\overline{\phi}\_{\delta\mathbf{k}}\overline{\phi}\_{\delta\mathbf{p}}\,\frac{\partial\overline{\nabla}\mathbf{p}\_{\mathbf{k}}}{\partial\mathbf{p}\_{\mathbf{k}}\partial\overline{\nabla}\mathbf{p}\_{\mathbf{k}}}\,\overline{\phi}\_{\delta\mathbf{p}}\overline{\phi}\_{\delta\mathbf{p}}\,\frac{\partial\overline{\nabla}\mathbf{p}\_{\mathbf{k}}}{\partial\mathbf{p}\_{\mathbf{k}}\partial\overline{\nabla}\mathbf{p}\_{\mathbf{k}}}\right) \\ &= \frac{1}{\delta\mathbf{p}}\frac{1}{2\pi}\int\_{0}^{\delta\mathbf{p}}\frac{1}{\delta\mathbf{p}}\Big{(}\overline{\phi}\_{\delta\mathbf{k}}\overline{\phi}\_{\delta\mathbf{k}}\,\frac{\partial\overline{\nabla}\mathbf{p}\_{\mathbf{k}}}{\partial\mathbf{p}\_{\mathbf{k}}\partial\overline{\nabla}\mathbf{p}\_{\mathbf{k}}}\,\overline{\phi}\_{\delta\mathbf{k}}\overline{\phi}\_{\delta\mathbf{k}}\,\frac{\partial\overline{\nabla}\mathbf{p}\_{\mathbf{k}}}\,\overline{\phi}\_{\delta\mathbf{k}}\Big{)} \\ &+ \Phi\_{\delta\mathbf{k}}\delta\_{\delta\mathbf{k}}\frac{\delta}{\delta\mathbf{p}}\Big{(}\overline{\phi}\_{\delta\mathbf{k}}\delta\_{\delta\mathbf{k}}\mathbf{p}\_{\delta}\,\frac{\partial$$

As shown above, the formulation of governing equation of interfacial velocity is derived. Then, the most strict formulation of transport equations of interfacial area concentration is given by conservation equation of interfacial area concentration (Eq.(32) , Eq(34), or Eq.(36)) and conservation equation of interfacial velocity (Eq.(55)). As shown in Eq.(55), the conservation equation of interfacial velocity consists of various correlation terms of fluctuating terms of velocity and local instant volume fractions. These correlation terms represent the turbulent transport of interfacial area, which reflects the interactions between gas liquid interface and turbulence of gas and liquid phases. Equation (55) represents such turbulence transport terms of interfacial area concentration. Accurate predictions of interfacial area transport can be possible by solving the transport equations derived here. However, Eq.(55)

Transport of Interfacial Area Concentration in Two-Phase Flow 101

turbulent velocity of gas phase. In bubbly flow, it is considered that turbulent mixing of gas liquid interface is proportional to bubble diameter, dB and the turbulent velocity of gas phase is proportional to the turbulent velocity of liquid phase. These assumptions were confirmed by experiment and analysis of turbulent diffusion of bubbles in bubbly flow (Kataoka and Serizawa (1991a)). Therefore, turbulent diffusion coefficient of interfacial area

ai 1 L B 1 L

empirical coefficient. For the case of turbulent diffusion of bubble, experimental data were well predicted assuming K1=1/3 For the case of turbulent diffusion of interfacial area concentration, there are no direct experimental data of turbulent diffusion. However, the diffusion of bubble is closely related to the diffusion of interfacial area (surface area of bubble). Therefore, as first approximation, the value of K1 for bubble diffusion can be

Equations (61) is based on the model of turbulent diffusion of interfacial area concentration. In this model, it is assumed that turbulence is isotropic. However, in the practical two-phase flow in the flow passages turbulence is not isotropic and averaged velocities and turbulent velocity have distribution in the radial direction of flow passage. In such non-isotropic turbulence, the correlation terms of turbulent fluctuation of velocity and interfacial area concentration given by Eq.(57) is largely dependent on anisotropy of turbulence field. Such non-isotropic turbulence is related to the various terms consisting of turbulent stress which appear in the right hand side of Eq.(55). Assuming that turbulent stress of gas phase is proportional to that of liquid phase and turbulence model in single phase flow, turbulent

> t L L LTP L L ij

Here, LTP is the turbulent diffusivity of momentum in gas-liquid two-phase flow. For bubbly flow, this turbulent diffusivity is given by various researchers (Kataoka and

> LTP B L <sup>1</sup> d v' 3

Based on the model of turbulent stress in gas-liquid two-phase flow and Eq.(55), it is assumed that turbulent diffusion of interfacial area concentration due to non-isotropic turbulence is proportional to the velocity gradient of liquid phase. For the diffusion of bubble due to non-isotropic turbulence in bubbly flow in pipe, Kataoka and Serizawa (1991b,1993) proposed the following correlation based on the analysis of radial distributions

<sup>2</sup> v v { v ( v )} k3 (62)

(63)

D K v d 6K v

i

a

Here, L is the length scale of turbulent mixing of gas liquid interface and vG

concentration is assumed by following equation.

Here, is the averaged void fraction and vL

stress is given by

Serizawa (1991b,1993)).

of void fraction and bubble number density.

applied to diffusion of interfacial area concentration in bubbly flow.

D vL ai G (60)

(61)

is the turbulent velocity of liquid phase. K1 is

is the

consists of complicated correlation terms of fluctuating terms of local instant volume fraction, velocity, pressure and shear stress. The detailed knowledge of these correlation terms is not available. Therefore, solving Eq.(55) together with basic equations of two-fluid model is difficult at present. More detailed analytical and experimental works on turbulence transport terms of interfacial area concentration are necessary for solving practically Eq.(55).

#### **4. Constitutive equations of transport equations of interfacial area concentration. Source and sink terms, diffusion term, turbulence transport term**

As shown in the previous section, the rigorous formulation of transport equation of interfacial area concentration are given by conservation equation of interfacial area concentration (Eq.(32), Eq.(34) or Eq(36)) and conservation equation of interfacial velocity (Eq.(55)). However, Eq.(55) consists of complicated correlation terms of fluctuating terms of local instant volume fraction, velocity, pressure and shear stress. The detailed knowledge of these correlation terms is not available. Therefore, solving Eq.(55) together with basic equations of two-fluid model is difficult at present. More detailed analytical and experimental works on turbulence transport terms of interfacial area concentration are necessary for solving practically Eq.(55). From Eqs.(45) and (46), interfacial velocity is related to averaged velocity of phase k (gas phase or liquid phase)by following equation.

$$\overleftarrow{\mathbf{V}\_{\mathbf{i}}\mathbf{a}\_{\mathbf{i}}} = \overleftarrow{\mathbf{v}\_{\mathbf{k}}\mathbf{a}\_{\mathbf{i}}} - \frac{1}{2\pi} \int\_{0}^{2\pi} \int\_{0}^{\pi} \underset{\mathbf{l}\left|\to 0\right[}{\lim} \frac{1}{\left|\mathbf{r}\right|} (\overline{\phi\_{\mathbf{k}}' \mathbf{q}\_{\mathbf{k}\mathbf{r}} \mathbf{v}\_{\mathbf{k}\mathbf{r}}'} + \overline{\phi\_{\mathbf{k}} \phi\_{\mathbf{k}\mathbf{r}}' \mathbf{v}\_{\mathbf{k}}'}) \sin\theta \,\mathbf{d}\theta \,\mathbf{d}\phi \tag{56}$$

When one considers bubbly flow and phase k is gas phase, Eq.(56) can be rewritten by

$$\overline{\mathbf{V}\_{\rm i}\mathbf{a}\_{i}} = \overline{\mathbf{v}\_{\rm G}\mathbf{a}\_{i}} - \frac{1}{2\pi} \int\_{0}^{2\pi} \int\_{0}^{\pi} \lim\_{\|\mathbf{r}\| \to 0} \frac{1}{|\mathbf{r}|} \overline{\langle \boldsymbol{\phi}\_{\rm G}^{\prime} \boldsymbol{\phi}\_{\rm G} \mathbf{v}\_{\rm G}^{\prime} + \overline{\boldsymbol{\phi}\_{\rm G} \boldsymbol{\phi}\_{\rm G}^{\prime} \mathbf{v}\_{\rm G}^{\prime}} \rangle \sin \theta \,\mathbf{d}\theta \,\mathbf{d}\phi \tag{57}$$

From Eqs.(28) and (29) , following relation is derived.

$$\mathbf{a}'\_{\mathbf{i}} = \mathbf{a}\_{\mathbf{i}} - \overline{\mathbf{a}\_{\mathbf{i}}} = -\frac{1}{2\pi} \int\_{0}^{2\pi} \int\_{0}^{\pi} \lim\_{\|\mathbf{r}\| \to 0} \frac{1}{\|\mathbf{r}\|} \Big\{ \phi'\_{\mathbf{k}} \left( 1 - 2\overline{\phi\_{\mathbf{k}\mathbf{r}}} \right) + \phi'\_{\mathbf{k}\mathbf{r}} \left( 1 - 2\overline{\phi\_{\mathbf{k}}} \right) - 2 \langle \phi'\_{\mathbf{k}} \phi'\_{\mathbf{k}\mathbf{r}} - \overline{\phi'\_{\mathbf{k}} \phi'\_{\mathbf{k}\mathbf{r}}} \rangle \Big\} \sin\theta d\theta d\phi \tag{58}$$

In Eq.(57), the terms, <sup>G</sup> /r and Gr / r are related to the fluctuating term of interfacial area concentration. On the other hand, the terms, Gr Gr v and G G v are the fluctuating term of gas phase velocity at the location, **x**+**r.** and **x**. Therefore, the second term of right hand side of Eq.(57) is considered to correspond to turbulent transport term due to the turbulent velocity fluctuation. In analogous to the turbulent transport of momentum, energy (temperature) and mass, the correlation term described above is assumed to be proportional to the gradient of interfacial area concentration which is transported by turbulence (diffusion model). Then, one can assume following relation.

$$-\frac{1}{2\pi} \int\_0^{2\pi} \int\_0^{\pi} \lim\_{\left|\mathbf{r}\right| \to 0} \frac{1}{\left|\mathbf{r}\right|} (\overline{\phi\_{\rm G}' \phi\_{\rm G r} \mathbf{v}\_{\rm G r}'} + \overline{\phi\_{\rm G} \phi\_{\rm G r}' \mathbf{v}\_{\rm G}'}) \sin\theta \mathbf{d}\theta d\phi = -\mathbf{D}\_{\rm ai} \mathbf{grad}\_{\mathbf{i}} \tag{59}$$

Here, the coefficient, Dai is considered to correspond to turbulent diffusion coefficient of interfacial area concentration. In analogy to the turbulent transport of momentum, energy (temperature) and mass, this coefficient is assumed to be given by

consists of complicated correlation terms of fluctuating terms of local instant volume fraction, velocity, pressure and shear stress. The detailed knowledge of these correlation terms is not available. Therefore, solving Eq.(55) together with basic equations of two-fluid model is difficult at present. More detailed analytical and experimental works on turbulence transport

terms of interfacial area concentration are necessary for solving practically Eq.(55).

**4. Constitutive equations of transport equations of interfacial area** 

2

2

<sup>G</sup> /r and

(diffusion model). Then, one can assume following relation.

(temperature) and mass, this coefficient is assumed to be given by

2

From Eqs.(28) and (29) , following relation is derived.

In Eq.(57), the terms,

ii ki k kr kr k kr k 0 0 r 0 1 1 Va v a lim ( v v )sin d d 2 r

When one considers bubbly flow and phase k is gas phase, Eq.(56) can be rewritten by

1 1 Va v a lim ( v v )sin d d 2 r

<sup>2</sup>

(58)

area concentration. On the other hand, the terms, Gr Gr v and G G v are the fluctuating term of gas phase velocity at the location, **x**+**r.** and **x**. Therefore, the second term of right hand side of Eq.(57) is considered to correspond to turbulent transport term due to the turbulent velocity fluctuation. In analogous to the turbulent transport of momentum, energy (temperature) and mass, the correlation term described above is assumed to be proportional to the gradient of interfacial area concentration which is transported by turbulence

> G Gr Gr G Gr G ai i 0 0 r 0 1 1 lim ( v v )sin d d D grada 2 r

Here, the coefficient, Dai is considered to correspond to turbulent diffusion coefficient of interfacial area concentration. In analogy to the turbulent transport of momentum, energy

(59)

(56)

(57)

Gr / r are related to the fluctuating term of interfacial

ii Gi G Gr Gr G Gr G 0 0 r 0

i ii k kr kr k k kr k kr 0 0 r 0 1 1 aaa lim (1 2 ) (1 2 ) 2( ) sin d d 2 r

**term** 

**concentration. Source and sink terms, diffusion term, turbulence transport** 

As shown in the previous section, the rigorous formulation of transport equation of interfacial area concentration are given by conservation equation of interfacial area concentration (Eq.(32), Eq.(34) or Eq(36)) and conservation equation of interfacial velocity (Eq.(55)). However, Eq.(55) consists of complicated correlation terms of fluctuating terms of local instant volume fraction, velocity, pressure and shear stress. The detailed knowledge of these correlation terms is not available. Therefore, solving Eq.(55) together with basic equations of two-fluid model is difficult at present. More detailed analytical and experimental works on turbulence transport terms of interfacial area concentration are necessary for solving practically Eq.(55). From Eqs.(45) and (46), interfacial velocity is related to averaged velocity of phase k (gas phase or liquid phase)by following equation.

$$\mathbf{D}\_{\rm ai} \propto \left| \mathbf{v}\_{\rm G}' \right| \mathbf{L} \tag{60}$$

Here, L is the length scale of turbulent mixing of gas liquid interface and vG is the turbulent velocity of gas phase. In bubbly flow, it is considered that turbulent mixing of gas liquid interface is proportional to bubble diameter, dB and the turbulent velocity of gas phase is proportional to the turbulent velocity of liquid phase. These assumptions were confirmed by experiment and analysis of turbulent diffusion of bubbles in bubbly flow (Kataoka and Serizawa (1991a)). Therefore, turbulent diffusion coefficient of interfacial area concentration is assumed by following equation.

$$\mathbf{^\ast D\_{ai} = K\_1 \left| \mathbf{v\_L'} \right|} \mathbf{d\_B} = 6 \mathbf{K\_1} \frac{\alpha}{\mathbf{a\_i}} \left| \mathbf{v\_L'} \right| \tag{61}$$

Here, is the averaged void fraction and vL is the turbulent velocity of liquid phase. K1 is empirical coefficient. For the case of turbulent diffusion of bubble, experimental data were well predicted assuming K1=1/3 For the case of turbulent diffusion of interfacial area concentration, there are no direct experimental data of turbulent diffusion. However, the diffusion of bubble is closely related to the diffusion of interfacial area (surface area of bubble). Therefore, as first approximation, the value of K1 for bubble diffusion can be applied to diffusion of interfacial area concentration in bubbly flow.

Equations (61) is based on the model of turbulent diffusion of interfacial area concentration. In this model, it is assumed that turbulence is isotropic. However, in the practical two-phase flow in the flow passages turbulence is not isotropic and averaged velocities and turbulent velocity have distribution in the radial direction of flow passage. In such non-isotropic turbulence, the correlation terms of turbulent fluctuation of velocity and interfacial area concentration given by Eq.(57) is largely dependent on anisotropy of turbulence field. Such non-isotropic turbulence is related to the various terms consisting of turbulent stress which appear in the right hand side of Eq.(55). Assuming that turbulent stress of gas phase is proportional to that of liquid phase and turbulence model in single phase flow, turbulent stress is given by

$$
\overline{\mathbf{v\_L'v\_L'}} = -\varepsilon\_{\rm LTP} \{ \overline{\nabla \mathbf{v\_L}} + \mathbf{\overline{(\nabla \mathbf{v\_L})}} \} + \frac{2}{3} \mathbf{k} \delta\_{\overline{\mathbf{u}}} \tag{62}
$$

Here, LTP is the turbulent diffusivity of momentum in gas-liquid two-phase flow. For bubbly flow, this turbulent diffusivity is given by various researchers (Kataoka and Serizawa (1991b,1993)).

$$\varepsilon\_{\rm LTP} = \frac{1}{3} \mathbf{c} \mathbf{d}\_{\rm B} \left| \mathbf{v}'\_{\rm L} \right| \tag{63}$$

Based on the model of turbulent stress in gas-liquid two-phase flow and Eq.(55), it is assumed that turbulent diffusion of interfacial area concentration due to non-isotropic turbulence is proportional to the velocity gradient of liquid phase. For the diffusion of bubble due to non-isotropic turbulence in bubbly flow in pipe, Kataoka and Serizawa (1991b,1993) proposed the following correlation based on the analysis of radial distributions of void fraction and bubble number density.

Transport of Interfacial Area Concentration in Two-Phase Flow 103

Similarly, the transport equation of interfacial area concentration for droplet flow is given by

ai i iL G G

Here D/Dt denotes material derivative following the liquid phase motion andCO and Bk are sink and source terms due to droplet coalescence and break up and L is the mass generation rate of liquid phase per unit volume of two-phase flow due to condensation. Here, turbulent diffusion coefficient of interfacial area concentration is approximated by turbulent diffusion coefficient of droplet (Cousins and Hewitt (1968)) as first approximation. The coefficient C for turbulent diffusion of interfacial area concentration due to nonisotropic turbulence (or lift force term) can be approximated by lift force coefficient of solid

The research and development of source and sink terms in transport equation of interfacial area concentration have been carried out mainly for the bubbly flow based on detailed

Hibiki and Ishii(2000a,2002) developed the transport equation of interfacial area mentioned above and carried out detailed modeling of source and sink terms of interfacial area concentration. They assumed that the sink term of interfacial area concentration is mainly due to the coalescence of bubble. On the other hand, they assumed the source term is mainly contributed by the break up of bubble due to liquid phase turbulence. Based on detailed mechanistic modeling of bubble liquid interactions, they finally obtained the constitutive

The sink term of interfacial area concentration due to the coalescence of bubbles, CO is composed of number of collisions of bubbles per unit volume and the probability of

CO 11/3 C 3

collision. db, and are bubble diameter, turbulent dissipation and surface tension. max is maximum permissible void fraction in bubbly flow and assumed to be 0.52. C and KC are

i b max

a d( )

5 32 6 b L C 3

<sup>2</sup> 2 1/3 5 32 C 6 b L

<sup>d</sup> exp K

represents the number of collisions of bubbles per unit

represents the probability of coalescence at

analysis and experiment of interfacial area concentration which is shown below.

div(D garda ) div{C(1 )a (v v ) rot(v )}

L CO BK

(69)

(70)

i

sphere as first approximation.

coalescence at collision and given by

where the term

volume and the term

i L

<sup>a</sup> div(a v ) <sup>t</sup>

L

equations for sink and source terms of interfacial area transport.

2 1/3 C 11/3 d( ) b max 

empirical constants and following values are given

<sup>d</sup> exp K

i L

2 a <sup>D</sup> (1 ) 3 (1 ) Dt

$$\mathbf{J}\_{\rm B} = \mathbf{K}\_2 \alpha \mathbf{d}\_{\rm B} \mathbf{n}\_{\rm B} \frac{\partial \mathbf{v}\_{\rm L}}{\partial \mathbf{y}} \tag{64}$$

Here, JB is the bubble flux in radial direction and nB is the number density of bubble. y is radial distance from wall of flow passage. K2 is empirical coefficient and experimental data were well predicted assuming K2=10. In analogous to Eq.(64), it is assumed that turbulent diffusion of interfacial area concentration due to non-isotropic turbulence is given by following equation.

$$\mathbf{J}\_{\rm ai} = \mathbf{K}\_2 \mathbf{a} \mathbf{d}\_{\rm B} \mathbf{a}\_{\rm i} \frac{\partial \mathbf{v}\_{\rm L}}{\partial \mathbf{y}} \tag{65}$$

Here, Jai is the flux of interfacial area concentration in radial direction. Equation (64) can be interpreted as equation of bubble flux due to the lift force due to liquid velocity gradient.

As shown above, turbulent diffusion of interfacial area concentration due to non-isotropic turbulence is related to the gradient of averaged velocity of liquid phase and using analogy to the lift force of bubble, Eq.(58) can be rewritten in three dimensional form by

$$\begin{aligned} &-\frac{1}{2\pi} \int\_0^{2\pi} \int\_0^{\pi} \varlimsup\_{\left|\mathbf{r}\right| \to 0} \overline{\frac{1}{|\mathbf{r}|} (\overline{\phi\_{\rm G}' \phi\_{\rm G r} \mathbf{v}'\_{\rm G r}} + \overline{\phi\_{\rm G} \phi\_{\rm G r}' \mathbf{v}'\_{\rm G}})} \sin\theta \mathbf{d}\theta \mathbf{d}\phi = \\ &- \mathbf{D}\_{\rm ai} \overline{\mathbf{grad}\_{\mathbf{i}} + \mathbf{C} \mathbf{ca}\_{\rm i} (\mathbf{v}\_{\rm G} - \mathbf{v}\_{\rm L})} \times \operatorname{rot}(\mathbf{v}\_{\rm L}) \end{aligned} \tag{66}$$

Empirical coefficient C in the right hand side of Eq.(66) should be determined based on the experimental data of spatial distribution of interfacial area concentration and averaged velocity of each phase. However, at present, there are not sufficient experimental data. Therefore, as first approximation, the value of coefficient C can be given by Eq.(67)

$$\mathbf{C} = \mathbf{K}\_2 \mathbf{d}\_\mathcal{B} \;/\; \mathbf{u}\_\mathcal{R} \; \text{ (based on Eq.(65))}\tag{67}$$

Using Eqs(55) and (66), transport equation of interfacial area concentration (Eq.(32),(34) or (36)) can be given by following equation for gas-liquid two-phase flow where gas phase is dispersed in liquid phase for bubbly flow.

$$\begin{aligned} \frac{\partial \overline{\mathbf{a}\_i}}{\partial \mathbf{t}} + \text{div}(\overline{\mathbf{a}\_i \mathbf{v}\_G}) &= \text{div}(\mathbf{D}\_{\text{ail}} \text{grad} \overline{\mathbf{a}\_i}) - \text{div} \{ \text{Caa}\_i (\mathbf{v}\_G - \mathbf{v}\_L) \times \text{rot}(\mathbf{v}\_L) \} + \\ + \frac{2}{3} \frac{\overline{\mathbf{a}\_i}}{\text{a} \rho\_G} \left( \Gamma\_G - \alpha \frac{\text{D} \rho\_G}{\text{D}t} \right) + \phi\_{\text{CO}} + \phi\_{\text{BK}} \end{aligned} \tag{68}$$

Here, D/Dt denotes material derivative following the gas phase motion and turbulent diffusion coefficient of interfacial area concentration, Dai is given by Eq.(61). Coefficient of turbulent diffusion of interfacial area concentration due to non-isotropic turbulence, C is given by Eq.(67). The third term in the right hand side of Eq.(68) is source term of interfacial area concentration due to phase change and density change of gas phase due to pressure change. G is the mass generation rate of gas phase per unit volume of two-phase flow due to evaporation. CO and Bk are sink and source term due to bubble coalescence and break up

Here, JB is the bubble flux in radial direction and nB is the number density of bubble. y is radial distance from wall of flow passage. K2 is empirical coefficient and experimental data were well predicted assuming K2=10. In analogous to Eq.(64), it is assumed that turbulent diffusion of interfacial area concentration due to non-isotropic turbulence is given by

B 2 BB <sup>v</sup> J K dn

ai 2 B i <sup>v</sup> J K da

Here, Jai is the flux of interfacial area concentration in radial direction. Equation (64) can be interpreted as equation of bubble flux due to the lift force due to liquid velocity gradient.

As shown above, turbulent diffusion of interfacial area concentration due to non-isotropic turbulence is related to the gradient of averaged velocity of liquid phase and using analogy

to the lift force of bubble, Eq.(58) can be rewritten in three dimensional form by

ai i i G L L

Therefore, as first approximation, the value of coefficient C can be given by Eq.(67)

D grada C a (v v ) rot(v )

G Gr Gr G Gr G 0 0 r 0

1 1 lim ( v v )sin d d 2 r

Empirical coefficient C in the right hand side of Eq.(66) should be determined based on the experimental data of spatial distribution of interfacial area concentration and averaged velocity of each phase. However, at present, there are not sufficient experimental data.

Using Eqs(55) and (66), transport equation of interfacial area concentration (Eq.(32),(34) or (36)) can be given by following equation for gas-liquid two-phase flow where gas phase is

<sup>a</sup> div(a v ) div(D grada ) div{C a (v v ) rot(v )} <sup>t</sup>

Here, D/Dt denotes material derivative following the gas phase motion and turbulent diffusion coefficient of interfacial area concentration, Dai is given by Eq.(61). Coefficient of turbulent diffusion of interfacial area concentration due to non-isotropic turbulence, C is given by Eq.(67). The third term in the right hand side of Eq.(68) is source term of interfacial area concentration due to phase change and density change of gas phase due to pressure change. G is the mass generation rate of gas phase per unit volume of two-phase flow due to evaporation. CO and Bk are sink and source term due to bubble coalescence and break up

i G ai i iG L L

2

dispersed in liquid phase for bubbly flow.

G

i G

2 a D 3 Dt

G CO BK

i

following equation.

L

(64)

(65)

(68)

y 

L

(66)

C K d /u 2B R (based on Eq.(65)) (67)

y 

Similarly, the transport equation of interfacial area concentration for droplet flow is given by

$$\begin{aligned} \frac{\partial \mathbf{a}\_{i}}{\partial t} + \text{div}(\overline{\mathbf{a}\_{i} \mathbf{v}\_{L}}) &= \\ \text{div}(\mathbf{D}\_{\text{ai}} \text{grad} \overline{\mathbf{a}\_{i}}) - \text{div}\{\mathbf{C}(1-\alpha)\overline{\mathbf{a}\_{i}}(\overline{\mathbf{v}\_{L}} - \overline{\mathbf{v}\_{G}}) \times \text{rot}(\overline{\mathbf{v}\_{G}})\} + \\ + \frac{2}{3} \frac{\overline{\mathbf{a}\_{i}}}{(1-\alpha)\rho\_{L}} \Big(\Gamma\_{L} - (1-\alpha)\frac{\mathbf{D}\rho\_{L}}{\mathbf{D}t}\Big) + \phi\_{\text{CO}} + \phi\_{\text{BK}} \end{aligned} \tag{69}$$

Here D/Dt denotes material derivative following the liquid phase motion andCO and Bk are sink and source terms due to droplet coalescence and break up and L is the mass generation rate of liquid phase per unit volume of two-phase flow due to condensation. Here, turbulent diffusion coefficient of interfacial area concentration is approximated by turbulent diffusion coefficient of droplet (Cousins and Hewitt (1968)) as first approximation. The coefficient C for turbulent diffusion of interfacial area concentration due to nonisotropic turbulence (or lift force term) can be approximated by lift force coefficient of solid sphere as first approximation.

The research and development of source and sink terms in transport equation of interfacial area concentration have been carried out mainly for the bubbly flow based on detailed analysis and experiment of interfacial area concentration which is shown below.

Hibiki and Ishii(2000a,2002) developed the transport equation of interfacial area mentioned above and carried out detailed modeling of source and sink terms of interfacial area concentration. They assumed that the sink term of interfacial area concentration is mainly due to the coalescence of bubble. On the other hand, they assumed the source term is mainly contributed by the break up of bubble due to liquid phase turbulence. Based on detailed mechanistic modeling of bubble liquid interactions, they finally obtained the constitutive equations for sink and source terms of interfacial area transport.

The sink term of interfacial area concentration due to the coalescence of bubbles, CO is composed of number of collisions of bubbles per unit volume and the probability of coalescence at collision and given by

$$\boldsymbol{\phi}\_{\rm CO} = -\left(\frac{\boldsymbol{\alpha}}{\mathbf{a}\_{\rm i}}\right)^{2} \frac{\Gamma\_{\rm C} \boldsymbol{\alpha}^{2} \,\mathrm{c}^{1/3}}{\mathbf{d}\_{\rm b}^{11/3} (\boldsymbol{a}\_{\rm max} - \boldsymbol{\alpha})} \exp\left(-\mathbf{K}\_{\rm C} \sqrt[6]{\frac{\mathbf{d}\_{\rm b}^{5} \boldsymbol{\rho}\_{\rm L}^{3} \,\mathrm{c}^{2}}{\sigma^{3}}}\right) \tag{70}$$

where the term 2 1/3 C 11/3 d( ) b max represents the number of collisions of bubbles per unit

volume and the term 5 32 6 b L C 3 <sup>d</sup> exp K represents the probability of coalescence at

collision. db, and are bubble diameter, turbulent dissipation and surface tension. max is maximum permissible void fraction in bubbly flow and assumed to be 0.52. C and KC are empirical constants and following values are given

Transport of Interfacial Area Concentration in Two-Phase Flow 105

Lb b 2 (d) d We

On the other hand, source term due to break up of bubble by liquid phase turbulence is

(1 ) <sup>1</sup> 1.24 60.3 exp a d 1 0.42(1 ) We /1.24 We

The transport equation of interfacial area concentration and constitutive equations of source terms described above are implemented to CATHARE code which is developed at CEA using three-dimensional two-fluid model and k- turbulence model. Predictions were carried out for thermal hydrodynamic structure of boiling and non-boiling (air-water) twophase flow including of interfacial area concentration. Comparisons were made with experimental data of DEBORA experiment which is boiling experiment using R-12 carried out at CEA and DEDALE experiment which is air-water experiment carried out at EDF, Electricite de France. The predictions reasonably agreed with experimental data of boiling and non-boiling two-phase flow for distribution of void fraction, velocities of gas and liquid phase, turbulent velocity and interfacial area concentration and the validity of transport

The measurements of interfacial area have been carried out earlier in the field of chemical engineering using chemical reaction and/or chemical absorption at gas-liquid interface (Sharma and Danckwerts (1970)). A lot of experimental studies have been reported and reviewed (Ishii et al.,(1982), Kocamustafaogullari and Ishii (1983)). However, in this method, measured quantity is the product of interfacial area concentration and mass transfer coefficient. Light attenuation method and photographic method were also developed and measurement of interfacial area concentration was carried out. However, the measured interfacial area concentration using these methods is volumetric averaged value and measurement of local interfacial area concentration is impossible. In the detailed analysis of multidimensional two-phase flow, measurements of distribution of local interfacial area concentration are indispensable for the validation of interfacial area transport model. Therefore, the establishment of the measurement method of local interfacial area concentration was strongly required. Ishii (1975) and Delhaye (1968) derived following relation among time averaged interfacial area concentration, number of interfaces and velocity of interface. They pointed out local interfacial area concentration can be measured

N

j 1 1 1

**n v ij ij**

Here, T and N are time interval of measurement and number of interfaces passing a measuring point during time interval T. **n**ij and **v**ij are unit normal vector and interfacial velocity of j-th interface. For bubbly flow, assuming that shape of bubble is spherical and sensor of probe passes any part of bubble with equal probability, Eq.(78) can be simplified to

i

a T

<sup>2</sup> 1/3

equation and constitutive equations described above was confirmed.

using two or three sensor probe based on this relation.

**5. Experimental researches on interfacial area concentration** 

BK 11/3 i b

given by

2/3

(76)

(78)

(77)

$$
\Gamma\_{\mathbb{C}} = 0.188, \ K\_{\mathbb{C}} = 1.29 \tag{71}
$$

As for the source term due to the break up of bubble, it is assumed that bubble break up mainly occurs due to the collision between bubble and turbulence eddy of liquid phase. The constitutive equation is given based on the detailed mechanistic modeling of this phenomenon as

$$\boldsymbol{\phi}\_{\rm BK} = \left(\frac{\boldsymbol{\alpha}}{\mathbf{a}\_{\rm i}}\right)^{2} \frac{\Gamma\_{\rm B} \alpha (1 - \alpha) \varepsilon^{1/3}}{\mathbf{d}\_{\rm b}^{11/3} (\alpha\_{\rm max} - \alpha)} \exp\left(-\frac{\mathbf{K}\_{\rm B} \sigma}{\rho\_{\rm L} \mathbf{d}\_{\rm b}^{5/3} \varepsilon^{2/3}}\right) \tag{72}$$

where the term 1/3 B 11/3 b max (1 ) d( ) represents the number of collisions of bubble and

turbulence eddy per unit volume and the term <sup>B</sup> 5/3 2/3 L b K exp <sup>d</sup> represents the

probability of break up at collision. B and KB are empirical constants and following values are given.

$$
\Gamma\_{\rm B} = 0.264, \,\mathrm{K\_{B}} = 1.37 \,\mathrm{} \tag{73}
$$

The validity of transport equation of interfacial area concentration (Eq.(68)) and constitutive equations for sink term due to bubble coalescence (Eq.(70)) and source term due to bubble break up (Eq.(72)) are confirmed by experimental data as will be described later in details.

Hibiki and Ishii (2000b) further modified their model of interfacial area transport and applied to bubbly-to-slug flow transition. In bubbly-to-slug flow transition, bubbles are classified into two groups that are small spherical/distorted bubble (group I) and large cap/slug bubble (group II). They derived transport equations of interfacial area concentration and constitutive equations for sink and source terms for group I and group II bubbles based on the transport equation and constitutive equations for bubbly flow mentioned above.

Yao and More (2004) developed more practical transport equation of interfacial area concentration and constitutive correlations of source terms. They derived these equations based on the basic transport equation developed in CEA and models of source terms developed at Purdue University (Ishii 's group). They also developed sink term due to coalescence of bubbles which is given by

$$\phi\_{\rm CO} = -107.8 \left( \frac{a}{\text{a}\_{\rm i}} \right)^2 \frac{\alpha^2 \varepsilon^{1/3}}{\text{d}\_{\rm b}^{11/3}} \frac{1}{\text{g}(a) + 1.922a \sqrt{\text{We} / 1.24}} \exp \left( -1.017 \sqrt{\frac{\text{We}}{1.24}} \right) \tag{74}$$

where g() and We is given by

$$\mathbf{g}(\alpha) = \frac{(\alpha\_{\text{max}}^{1/3} - \alpha^{1/3})}{\alpha\_{\text{max}}^{1/3}} \qquad \text{ ( $\alpha\_{\text{max}} = 0.52$ )}\tag{75}$$

As for the source term due to the break up of bubble, it is assumed that bubble break up mainly occurs due to the collision between bubble and turbulence eddy of liquid phase. The constitutive equation is given based on the detailed mechanistic modeling of this

a d( ) d

probability of break up at collision. B and KB are empirical constants and following values

The validity of transport equation of interfacial area concentration (Eq.(68)) and constitutive equations for sink term due to bubble coalescence (Eq.(70)) and source term due to bubble break up (Eq.(72)) are confirmed by experimental data as will be described later in details.

Hibiki and Ishii (2000b) further modified their model of interfacial area transport and applied to bubbly-to-slug flow transition. In bubbly-to-slug flow transition, bubbles are classified into two groups that are small spherical/distorted bubble (group I) and large cap/slug bubble (group II). They derived transport equations of interfacial area concentration and constitutive equations for sink and source terms for group I and group II bubbles based on the transport equation and constitutive equations for bubbly flow

Yao and More (2004) developed more practical transport equation of interfacial area concentration and constitutive correlations of source terms. They derived these equations based on the basic transport equation developed in CEA and models of source terms developed at Purdue University (Ishii 's group). They also developed sink term due to

> <sup>1</sup> We 107.8 exp 1.017 a d g( ) 1.922 We /1.24 1.24

> > 1/3 max

( ) g( ) ( 0.52) 

(74)

(75)

1/3 1/3 max

max

B B BK 11/3 5/3 2/3 i b max L b

(1 ) K exp

(72)

represents the number of collisions of bubble and

exp <sup>d</sup>

5/3 2/3

represents the

L b K

B =0.264, KB =1.37 (73)

<sup>2</sup> 1/3

1/3

turbulence eddy per unit volume and the term <sup>B</sup>

B 11/3 b max (1 ) d( ) 

phenomenon as

where the term

mentioned above.

coalescence of bubbles which is given by

where g() and We is given by

CO 11/3

<sup>2</sup> 2 1/3

i b

are given.

C =0.188, KC =1.29 (71)

$$\text{We} = \frac{2\rho\_{\text{L}} (\text{\text{\textdegree\text{d}\_{b}}})^{2/3} \text{d}\_{\text{b}}}{\sigma} \tag{76}$$

On the other hand, source term due to break up of bubble by liquid phase turbulence is given by

$$\phi\_{\rm BK} = 60.3 \left( \frac{a}{\mathbf{a\_i}} \right)^2 \frac{\mathbf{c^{1/3}} a (1 - a)}{\mathbf{d\_b}^{11/3}} \frac{1}{1 + 0.42 (1 - a) \sqrt{\rm We / 1.24}} \exp \left( - \sqrt{\frac{1.24}{\rm We}} \right) \tag{77}$$

The transport equation of interfacial area concentration and constitutive equations of source terms described above are implemented to CATHARE code which is developed at CEA using three-dimensional two-fluid model and k- turbulence model. Predictions were carried out for thermal hydrodynamic structure of boiling and non-boiling (air-water) twophase flow including of interfacial area concentration. Comparisons were made with experimental data of DEBORA experiment which is boiling experiment using R-12 carried out at CEA and DEDALE experiment which is air-water experiment carried out at EDF, Electricite de France. The predictions reasonably agreed with experimental data of boiling and non-boiling two-phase flow for distribution of void fraction, velocities of gas and liquid phase, turbulent velocity and interfacial area concentration and the validity of transport equation and constitutive equations described above was confirmed.

#### **5. Experimental researches on interfacial area concentration**

The measurements of interfacial area have been carried out earlier in the field of chemical engineering using chemical reaction and/or chemical absorption at gas-liquid interface (Sharma and Danckwerts (1970)). A lot of experimental studies have been reported and reviewed (Ishii et al.,(1982), Kocamustafaogullari and Ishii (1983)). However, in this method, measured quantity is the product of interfacial area concentration and mass transfer coefficient. Light attenuation method and photographic method were also developed and measurement of interfacial area concentration was carried out. However, the measured interfacial area concentration using these methods is volumetric averaged value and measurement of local interfacial area concentration is impossible. In the detailed analysis of multidimensional two-phase flow, measurements of distribution of local interfacial area concentration are indispensable for the validation of interfacial area transport model. Therefore, the establishment of the measurement method of local interfacial area concentration was strongly required. Ishii (1975) and Delhaye (1968) derived following relation among time averaged interfacial area concentration, number of interfaces and velocity of interface. They pointed out local interfacial area concentration can be measured using two or three sensor probe based on this relation.

$$\overline{\mathbf{a}\_{\rm i}}^{\prime} = \frac{1}{\mathbf{T}} \sum\_{j=1}^{N} \frac{1}{\left| \mathbf{n}\_{\rm ij} \bullet \mathbf{v}\_{\rm ij} \right|} \tag{78}$$

Here, T and N are time interval of measurement and number of interfaces passing a measuring point during time interval T. **n**ij and **v**ij are unit normal vector and interfacial velocity of j-th interface. For bubbly flow, assuming that shape of bubble is spherical and sensor of probe passes any part of bubble with equal probability, Eq.(78) can be simplified to

Transport of Interfacial Area Concentration in Two-Phase Flow 107

Probe 2

**n**s2

The direction cosines of unit vector of each double sensor probe (**n**sk, as shown in Fig.9) are denoted by cosxk, cosyk, coszk. Then, the inverse of product of interfacial velocity and

Sensor 1

s

2

Sensor 2

1 AAA

0 x2 y2 z2

cos cos cos

A cos cos cos

cos cos cos

<sup>2</sup> i i <sup>0</sup>

n v A

2 2 2 123

x1 y1 z1

x3 y3 z3

(82)

s3

Probe 3

**n**s3

(83)

Fig. 9. Three double sensor probe (four sensor probe)

Probe 1

**n**s1

Here, |A0|, |A1|,|A2| and |A3| are given by

unit normal vector of interface which appears in Eq.(78) is given by

s1

$$\overline{\mathbf{a}\_{i}} = 4 \frac{\mathbf{N}}{\mathbf{T}} \frac{\overline{\mathbf{1}}}{\|\mathbf{v}\_{\text{sz}}\|}\tag{79}$$

Here, vsz is the z directional (flow directional) component of velocity of interface measured by double sensor probe as shown in Fig.8**.** vsz is obtained by

$$\mathbf{v}\_{\rm sx} = \frac{\Delta \mathbf{s}}{\Delta \mathbf{t}}\tag{80}$$

where s is spacing of two sensors (Fig.8) and t is the time interval where interface passes upstream sensor and downstream sensor.

Fig. 8. Double sensor probe and velocity of interface

Later, based on local instant formulation of interfacial area concentration, Kataoka et al. (1986) proposed three double sensor probe method (four sensor probe method) as shown in Fig.9. Using this method, time averaged interfacial area is measured without assuming spherical bubble and statistical behavior of bubbles. The passing velocities measured by each double sensor probe are denoted by vsk which are given by

$$\mathbf{v}\_{\text{ sk}} = \frac{\Delta \mathbf{s}\_{\text{k}}}{\Delta \mathbf{t}\_{\text{k}}} \quad \text{ (k=1,2,3)}\tag{81}$$

N 1 a 4

Here, vsz is the z directional (flow directional) component of velocity of interface measured

where s is spacing of two sensors (Fig.8) and t is the time interval where interface passes

Later, based on local instant formulation of interfacial area concentration, Kataoka et al. (1986) proposed three double sensor probe method (four sensor probe method) as shown in Fig.9. Using this method, time averaged interfacial area is measured without assuming spherical bubble and statistical behavior of bubbles. The passing velocities measured by

k

(k=1,2,3) (81)

k s v t

sk

sz

sz

s v t

T v (79)

(80)

i

by double sensor probe as shown in Fig.8**.** vsz is obtained by

upstream sensor and downstream sensor.

Fig. 8. Double sensor probe and velocity of interface

each double sensor probe are denoted by vsk which are given by

Fig. 9. Three double sensor probe (four sensor probe)

The direction cosines of unit vector of each double sensor probe (**n**sk, as shown in Fig.9) are denoted by cosxk, cosyk, coszk. Then, the inverse of product of interfacial velocity and unit normal vector of interface which appears in Eq.(78) is given by

$$\frac{1}{\left|\mathbf{n}\_{i}\bullet\mathbf{v}\_{i}\right|}=\frac{\sqrt{\left|\mathbf{A}\_{1}\right|^{2}+\left|\mathbf{A}\_{2}\right|^{2}+\left|\mathbf{A}\_{3}\right|^{2}}}{\sqrt{\left|\mathbf{A}\_{0}\right|^{2}}}\tag{82}$$

Here, |A0|, |A1|,|A2| and |A3| are given by

$$\begin{vmatrix} \mathbf{A}\_{0} \\ \end{vmatrix} = \begin{vmatrix} \cos \eta\_{\mathbf{x}1} & \cos \eta\_{\mathbf{y}1} & \cos \eta\_{\mathbf{z}1} \\ \cos \eta\_{\mathbf{x}2} & \cos \eta\_{\mathbf{y}2} & \cos \eta\_{\mathbf{z}2} \\ \cos \eta\_{\mathbf{x}3} & \cos \eta\_{\mathbf{y}3} & \cos \eta\_{\mathbf{z}3} \end{vmatrix} \tag{83}$$

Transport of Interfacial Area Concentration in Two-Phase Flow 109

Hibiki, Hognet and Ishii (1998) carried out more detailed analysis of configuration of gasliquid interface and double sensor probe and proposed more accurate formulation of

> 2 N 3

2 2 <sup>2</sup> <sup>0</sup> <sup>0</sup> z iz 1 ( /v ) <sup>3</sup> sin 2 <sup>1</sup> 2 2 1 3( / v )

Double sensor probe or three double sensor probe (four sensor probe) has finite spacing between sensors. In relation to sensor spacing and size of bubble, some measurement errors are inevitable. In order to evaluate such measurement errors, a numerical simulation method using Monte Carlo approach is proposed (Kataoka et al., (1994), Wu and Ishii (1999)) for sensitivity analysis of measurement errors of double sensor probe or three double sensor probe. Using this method, Wu and Ishii (1999) carried out comprehensive analysis of accuracy of interfacial area measurement using double sensor probe including the probability of missing bubbles. They obtained formulation of interfacial area concentration measurement similar to Eqs.(91) and (92). The method using Eqs.(89) and (90)

For adiabatic two-phase flow, many research groups all over the world, carried out measurements of interfacial area concentration mainly using double sensor or four sensor electrical resistivity probes. Most of experiments were carried out for vertical upward airwater two-phase flow in pipe. Some data were reported in annulus or downward flow. Flow regime covers bubbly flow to bubbly-to-slug transition. Some data are reported for annular flow. The experimental database of interfacial area concentration for non-boiling system

Measurement of interfacial area concentration in boiling two-phase flow is quite important in view of practical application to nuclear reactor technology. However, in boiling two-phase flow, measurement of interfacial area is much more difficult compared with the measurement in non-boiling two-phase flow because of the durability of electrical resistivity and optical probes in high temperature liquid. Therefore, the accumulation of experimental data in boiling system was not sufficient compared with those in non-boiling system. However, recently, based on the establishment of measurement method of interfacial area as described above and improvement of electrical resistivity and optical probes, detailed measurements of interfacial area concentration become possible and experimental works have been carried out by various research groups. Most of experiments are carried out in annulus test section where inner pipe is heated. However, recently, some experimental studies are reported in rod bundle geometry. The experimental database of interfacial area concentration for boiling system

 

j 1 szj 0 0 1 1 a 2 I( ) T v 3( sin ) 

z iz 0

0

<sup>2</sup> <sup>2</sup>

(91)

(92)

interfacial area concentration measurement using double sensor probe. It is given by

i 0

underestimated the interfacial area concentration up to 50%.

described above is summarized in Table 1 (Kataoka (2010)).

described above is summarized in Table 2 (Kataoka (2010)).

Here 0 is given by

$$\begin{vmatrix} \mathbf{A}\_1 \\ \mathbf{A}\_1 \end{vmatrix} = \begin{vmatrix} 1/\,\mathrm{v}\_{s1} & \cos\eta\_{\mathrm{y}1} & \cos\eta\_{\mathrm{z}1} \\ 1/\,\mathrm{v}\_{s2} & \cos\eta\_{\mathrm{y}2} & \cos\eta\_{\mathrm{z}2} \\ 1/\,\mathrm{v}\_{s3} & \cos\eta\_{\mathrm{y}3} & \cos\eta\_{\mathrm{z}3} \end{vmatrix} \tag{84}$$

$$\begin{vmatrix} \mathbf{A}\_{2} \end{vmatrix} = \begin{vmatrix} \cos \mathfrak{n}\_{\mathbf{x}1} & 1/\operatorname{v}\_{\mathbf{s}1} & \cos \mathfrak{n}\_{\mathbf{x}1} \\ \cos \mathfrak{n}\_{\mathbf{x}2} & 1/\operatorname{v}\_{\mathbf{s}2} & \cos \mathfrak{n}\_{\mathbf{x}2} \\ \cos \mathfrak{n}\_{\mathbf{x}3} & 1/\operatorname{v}\_{\mathbf{s}3} & \cos \mathfrak{n}\_{\mathbf{x}3} \end{vmatrix} \tag{85}$$

$$\begin{vmatrix} \mathbf{A}\_{3} \end{vmatrix} = \begin{vmatrix} \cos \mathfrak{n}\_{\mathbf{x}1} & \cos \mathfrak{n}\_{\mathbf{y}1} & 1 \end{vmatrix} \begin{pmatrix} \mathbf{v}\_{s1} \\ \mathbf{v}\_{s1} \\ \cos \mathfrak{n}\_{\mathbf{x}2} & \cos \mathfrak{n}\_{\mathbf{y}2} & 1 \end{pmatrix} \begin{pmatrix} \mathbf{v}\_{s1} \\ \mathbf{v}\_{s2} \\ \mathbf{v}\_{s3} \end{pmatrix} \tag{86}$$

When three double sensor probes as shown in Fig.9 are orthogonal (perpendicular to each other), Eq.(82) is simply given by

$$\frac{1}{\left|\mathbf{n\_i}\bullet\mathbf{v\_i}\right|} = \sqrt{\left(\frac{1}{\mathbf{v\_{s1}}}\right)^2 + \left(\frac{1}{\mathbf{v\_{s2}}}\right)^2 + \left(\frac{1}{\mathbf{v\_{s3}}}\right)^2} \tag{87}$$

Then time averaged interfacial area concentration is given by

$$\overline{\mathbf{a}\_{i}} = \frac{1}{\mathcal{T}} \sum\_{j=1}^{N} \sqrt{\left(\frac{1}{\mathbf{v}\_{s1j}}\right)^{2} + \left(\frac{1}{\mathbf{v}\_{s2j}}\right)^{2} + \left(\frac{1}{\mathbf{v}\_{s3j}}\right)^{2}} \tag{88}$$

Most of recent experimental works of local interfacial area measurement are carried out by double sensor probe or three double sensor probe (four sensor probe) using electrical resistivity probe or optical probe.

For practical application, Kataoka et al.(1986) further proposed a simplified expression of Eq.(88) for double sensor probe which is given by

$$\overline{\mathbf{a}\_{i}} = 4 \frac{1}{T} \sum\_{j=1}^{N} \sqrt{\left(\frac{1}{\mathbf{v}\_{\text{szj}}}\right)^{2}} \frac{1}{1 - \cot \frac{1}{2} a\_{0} \ln(\cos \frac{1}{2} a\_{0}) - \tan \frac{1}{2} a\_{0} \ln(\sin \frac{1}{2} a\_{0})} \tag{89}$$

where α0 is given by

$$\frac{\sin 2\alpha\_0}{2\alpha\_0} = \frac{1 - \left(\sigma\_x^{\,^2} / \left| \overline{\mathbf{v}\_{ix}} \right|^2 \right)}{1 + 3(\sigma\_x^{\,^2} / \left| \overline{\mathbf{v}\_{ix}} \right|^2)}\tag{90}$$

Here, v and iz <sup>z</sup> are the mean value and fluctuation of the z component interfacial velocity.

Hibiki, Hognet and Ishii (1998) carried out more detailed analysis of configuration of gasliquid interface and double sensor probe and proposed more accurate formulation of interfacial area concentration measurement using double sensor probe. It is given by

$$\overline{\mathbf{a}\_{i}} = 2 \frac{1}{T} \sum\_{j=1}^{N} \sqrt{\left(\frac{1}{\mathbf{v}\_{s\bar{x}j}}\right)^{2}} I(\mathbf{a}\_{0}) \frac{\mathbf{a}\_{0}^{\,^{\,3}}}{\Im(\mathbf{a}\_{0} - \sin \alpha\_{0})} \tag{91}$$

Here 0 is given by

108 Nuclear Reactors

1 s2 y2 z2

1 / v cos cos

2 x2 s2 z2

3 x2 y2 s2

A cos cos 1 / v

When three double sensor probes as shown in Fig.9 are orthogonal (perpendicular to each

i i s1 s2 s3 1 111 nv v v v

Then time averaged interfacial area concentration is given by

i

2

a

Eq.(88) for double sensor probe which is given by

N

N

j 1 s1j s2j s3j 11 1 1

Tv v v

Most of recent experimental works of local interfacial area measurement are carried out by double sensor probe or three double sensor probe (four sensor probe) using electrical

For practical application, Kataoka et al.(1986) further proposed a simplified expression of

j 1 szj 0 000 1 1 <sup>1</sup> a 4 T v <sup>1111</sup> 1 cot ln(cos ) tan ln(sin ) 2222

> z iz 0 <sup>2</sup> <sup>2</sup> <sup>0</sup> z iz 1 ( /v ) sin 2 2 1 3( / v )

Here, v and iz <sup>z</sup> are the mean value and fluctuation of the z component interfacial

<sup>2</sup> <sup>2</sup>

(89)

other), Eq.(82) is simply given by

resistivity probe or optical probe.

i

where α0 is given by

velocity.

cos cos 1 / v

cos cos 1 / v

cos 1 / v cos A cos 1 / v cos cos 1 / v cos

A 1 / v cos cos

1 / v cos cos

s1 y1 z1

(84)

(85)

(86)

(87)

(90)

s3 y3 z3

x1 s1 z1

 

x3 s3 z3

x1 y1 s1

x3 y3 s3

 

222

222

(88)

$$\frac{3}{2\alpha\_0^2} \left( 1 - \frac{\sin 2\alpha\_0}{2\alpha\_0} \right) = \frac{1 - \left( \mathbf{\sigma}\_x \,^2 \slash\begin{bmatrix} \mathbf{\overline{v}\_{ix}} \\ \mathbf{\overline{v}\_{ix}} \end{bmatrix} \right)}{1 + 3\left( \mathbf{\overline{\sigma}\_x \,^2 \slash\begin{bmatrix} \mathbf{\overline{v}\_{ix}} \end{bmatrix}} \right)} \tag{92}$$

Double sensor probe or three double sensor probe (four sensor probe) has finite spacing between sensors. In relation to sensor spacing and size of bubble, some measurement errors are inevitable. In order to evaluate such measurement errors, a numerical simulation method using Monte Carlo approach is proposed (Kataoka et al., (1994), Wu and Ishii (1999)) for sensitivity analysis of measurement errors of double sensor probe or three double sensor probe. Using this method, Wu and Ishii (1999) carried out comprehensive analysis of accuracy of interfacial area measurement using double sensor probe including the probability of missing bubbles. They obtained formulation of interfacial area concentration measurement similar to Eqs.(91) and (92). The method using Eqs.(89) and (90) underestimated the interfacial area concentration up to 50%.

For adiabatic two-phase flow, many research groups all over the world, carried out measurements of interfacial area concentration mainly using double sensor or four sensor electrical resistivity probes. Most of experiments were carried out for vertical upward airwater two-phase flow in pipe. Some data were reported in annulus or downward flow. Flow regime covers bubbly flow to bubbly-to-slug transition. Some data are reported for annular flow. The experimental database of interfacial area concentration for non-boiling system described above is summarized in Table 1 (Kataoka (2010)).

Measurement of interfacial area concentration in boiling two-phase flow is quite important in view of practical application to nuclear reactor technology. However, in boiling two-phase flow, measurement of interfacial area is much more difficult compared with the measurement in non-boiling two-phase flow because of the durability of electrical resistivity and optical probes in high temperature liquid. Therefore, the accumulation of experimental data in boiling system was not sufficient compared with those in non-boiling system. However, recently, based on the establishment of measurement method of interfacial area as described above and improvement of electrical resistivity and optical probes, detailed measurements of interfacial area concentration become possible and experimental works have been carried out by various research groups. Most of experiments are carried out in annulus test section where inner pipe is heated. However, recently, some experimental studies are reported in rod bundle geometry. The experimental database of interfacial area concentration for boiling system described above is summarized in Table 2 (Kataoka (2010)).

Transport of Interfacial Area Concentration in Two-Phase Flow 111

Table 2. Summary of Experimental Database of Interfacial Area Concentration for Boiling

In order to confirm the validity of transport equation of interfacial area, comparisons with experimental data were carried out mainly for bubbly flow and churn flow. The transport equation for bubbly flow is given by Eq.(68). This equation includes turbulent diffusion term

**6. Validation of interfacial area transport models by experimental data** 

System


Table 1. Summary of Experimental Database of Interfacial Area Concentration for Non-Boiling System

Table 1. Summary of Experimental Database of Interfacial Area Concentration for Non-

Boiling System


Table 2. Summary of Experimental Database of Interfacial Area Concentration for Boiling System
