**2. Numerical simulation method**

The fluid mixing phenomena dominates water and steam distribution in the fuel bundles and is one of the most important phenomena that affect boiling transition (BT) of BWRs. Surface tension at gas-liquid interface (interface), shapes of interfaces and pressure difference fluctuation between subchannels play important roles in these phenomena. Then an advanced simulation method developed in this study must evaluate surface tension, interface shape and pressure fluctuation accurately. To satisfy these requirements, we developed an advanced simulation code: TPFIT and an advanced interface tracking method. Detail of the TPFIT and the advanced interface tracking method are explained in this section.

#### **2.1 Governing equation and numerical simulation**

To calculate surface tension and interface shapes by realistic computational resources, the TPFIT adopt interface tracking method. In two-phase flow in rod bundles, compressibility of steam has important influence for pressure difference fluctuation between subchannels. In TPFIT, considering the time-dependent Navier-Stokes equation for two phase compressible flow, conservative equations of mass of liquid, mass of gas, momentum and energy are described as follows:

Mass of liquid:

$$\frac{\left\|\left(\rho f\right)\right\|\_{\mathbb{L}}}{\left\|t\right\|} + \frac{\left\|u\_{\mathbb{L}}\left(\rho f\right)\right\|\_{\mathbb{L}}}{\left\|\mathbf{x}\_{\mathbb{k}}\right\|} = 0 \tag{1}$$

Mass of gas:

$$\frac{\left\|\left(\rho f\right)\_{\boldsymbol{s}}\right\|\_{\boldsymbol{s}} + \frac{\left\|u\_{\boldsymbol{s}}\left(\rho f\right)\_{\boldsymbol{s}}\right\|\_{\boldsymbol{s}}}{\left\|\mathbf{x}\_{\boldsymbol{s}}\right\|} = 0 \tag{2}$$

Momentum:

$$\frac{\partial \mathbf{u}\_i}{\partial t} + \boldsymbol{\mu}\_k \frac{\partial \mathbf{u}\_i}{\partial \mathbf{x}\_k} = -\frac{\partial p}{\partial \mathbf{x}\_i} + \frac{1}{\rho} \frac{\partial \boldsymbol{\tau}\_{ik}}{\partial \mathbf{x}\_k} + \mathbf{g}\_i + \boldsymbol{\sigma}\_i \tag{3}$$

Energy:

288 Computational Simulations and Applications

entail great cost, development of a method that enables the thermal-hydraulic design of

It is expected that large scale numerical simulations (numerical simulations by large scale computer) would replace certain large-scale tests and thermal-hydraulic information, some of which is currently difficult or impossible to measure experimentally, would be obtained.

For this reason we developed an advanced thermal-hydraulic design method for BWRs using innovative two-phase flow simulation technology. For this, the following are required: (1) an advanced simulation method with high accuracy prediction, (2) verification of simulation method; and, (3) analytical method to confirm or modify the correlations by detailed numerical simulation data. In this study, we are developing the method to create or modify the two-phase flow correlations for the fluid mixing phenomena in BWR based on

The fluid mixing phenomena dominates water and steam distribution in the fuel bundles and is one of the most important phenomena that affect boiling transition (BT) of BWRs. Surface tension at gas-liquid interface (interface), shapes of interfaces and pressure difference fluctuation between subchannels play important roles in these phenomena. Then an advanced simulation method developed in this study must evaluate surface tension, interface shape and pressure fluctuation accurately. To satisfy these requirements, we developed an advanced simulation code: TPFIT and an advanced interface tracking method. Detail of the TPFIT and the advanced interface tracking method are explained in this

To calculate surface tension and interface shapes by realistic computational resources, the TPFIT adopt interface tracking method. In two-phase flow in rod bundles, compressibility of steam has important influence for pressure difference fluctuation between subchannels. In TPFIT, considering the time-dependent Navier-Stokes equation for two phase compressible flow, conservative equations of mass of liquid, mass of gas, momentum and energy are

<sup>0</sup> *l l <sup>k</sup>*

*f uf t x*

1 *i i ik*

*f uf t x*

 

 

*k*

0 *g g <sup>k</sup> k*

*k i i ki k u u <sup>p</sup> u g t xx x*

 

(1)

(2)

(3)

 

And new design method will be developed based on these numerical simulations.

BWRs without these actual size tests is desired.

advanced numerical simulation technic.

**2. Numerical simulation method** 

**2.1 Governing equation and numerical simulation** 

section.

described as follows: Mass of liquid:

Mass of gas:

Momentum:

$$\frac{\partial e}{\partial t} + u\_k \frac{\partial e}{\partial \mathbf{x}\_k} = -\frac{p}{\rho} \frac{\partial u\_k}{\partial \mathbf{x}\_k} + \frac{1}{\rho} \frac{\partial}{\partial \mathbf{x}\_k} \left(\mathcal{A} \frac{\partial T}{\partial \mathbf{x}\_k}\right) + q \tag{4}$$

where *u*, *p*, *e*, *T* are velocity, static pressure, internal energy and temperature, and *g* and in the momentum equation are gravity and surface tension force, respectively. In above equations, subscripts *g* and *l* are used to represent gas and liquid phases. The mass of liquid and gas phases in two phase flow are defined as following equation.

$$\begin{aligned} \left(\rho f\right)\_l &= \rho\_l f\_l\\ \left(\rho f\right)\_s &= \rho\_s f\_s \end{aligned} \tag{5}$$

In Eq. (5), *f* is volume fraction of fluid, and Volume fraction of gas phase is evaluated by use of volume fraction of liquid.

$$f\_x = \mathbf{1} \cdot f\_l \tag{6}$$

Two-phase fluid density is calculated as sum of mass of both phases:

$$
\rho = \left(\rho f\right)\_! + \left(\rho f\right)\_! \tag{7}
$$

Eqs. (1) and (2) are solved in the advanced interface tracking method described in sec.2.2. The momentum equation (Eq. (3)) is solved by the CIP (Cubic Interpolated Pseudo-particle) method (Yabe and Aoki, 1991). The energy equation (Eq. (4)) is used to obtain the Poison equation for static pressure. Temperature is estimated by means of a fluid property routine based on the static pressure and local density of both phases. The ILUCGS method is used to solve the Poison equation for static pressure. In the TPFIT code, a Cartesian coordinate system and staggered grid are used. The surface tension force in the momentum equation is estimated using the CSF model (Blackbill et al., 1992). In the CSF model, volume fraction of liquid *fl* is required and evaluated in the advanced interface tracking method, too. The local viscosities and thermal conductivities of liquid and gas were evaluated using solved static pressure and temperature fields based on the fluid property routine.

#### **2.2 Outline of advanced interface tracking method**

The fundamental concept of the advanced interface tracking method is quite simple. That is, a transported volume of liquid and gas between neighboring calculation control volumes during each time step is calculated through the movement of approximated gas-liquid interfaces, as estimated in the Lagrangian system. Schematic drawings of the major three operational steps within each time step in the two-dimensional case are shown in Figure 1. In the first step (reconstruction step), as shown in Figure 1 (a), a gas-liquid interface in each control area for calculations is reconstructed, taking account of the liquid fraction represented by it and its surroundings. In the reconstruction step, the gas-liquid interface is approximated by a linear function: *F*(**x**) as same as PLIC (Gueyffier, et al., 1999) method (see Fig.2 for 3-dimensional case). The function *F*(**x**) is expressed as follows:

$$F(\mathbf{x}) = \sum\_{i=1}^{n\_{\mathbf{x}}} a\_i \mathbf{x}\_i + b\_i \; \mathbf{x} = \left(\mathbf{x}\_1, \mathbf{x}\_2, \mathbf{x}\_3\right) \tag{8}$$

Development of Two-Phase Flow Correlation

**u** *t*

Control volume 1

interface tracking method (2-dimensional case).

*x*1

Fig. 2. Approximated fluid segment in control volume (3-dimensional case).

*-b/a*<sup>1</sup>

summing these transferred mass.

(CV1)

Liquid

for Fluid Mixing Phenomena in Boiling Water Reactor 291

Control volume 2

(b) Transportation step

Volume fraction

(c) Redistribution step

In transportation step, transfer volume of both phases is calculated. In the last step (redistribution step), as shown in Fig.1(c), transferred mass of both phases can be evaluated by multiplying density by transferred volume. Mass at new time step can be calculated by

*x*<sup>2</sup>

*-b/a*<sup>3</sup>

*x*3

*ii bxa*

0

*i*

3 1

Fig. 1. Schematic drawings of volume fraction transport calculation using the advanced

*x*3

*<sup>x</sup>*1 Transported volume *V* from CV1 to CV2

*x*2

*-b/a*<sup>2</sup>

Gas

(CV2)

In the equation, *nm* represents dimension of simulation, and is 2 or 3. *xi* is the coordinate position of the definition point of scalar quantities. A unit normal vector to the interface **a**= (*a*1,*a*2,*a*3) and segment *b* in Eq.(8) must be estimated.

To evaluate the function *f*(**x**), we assumed that the interface exists in the position where volume fraction of liquid *fl* equals to 0.5. By least-squares method (choose eight nearest neighbors for 2-dimensional case, 26 nearest neighbors for 3-dimensional case), linear function is obtained.

$$F(\mathbf{x}) = f\_l - \mathbf{0.5} = \sum\_{i=1}^{n\_\*} a\_i \mathbf{x}\_i + b\_0 \tag{9}$$

In the equation, *b*0 is segment of linear function evaluated by least square method. Control volume is divided into two polyhedrons by the linear functions. If gray polyhedron in Fig.2 corresponds to liquid region in control volume, volume of this polyhedron: *Vl* must be equals to proper volume of liquid. Thus,

$$V\_l = f\_l \Delta V \tag{10}$$

where *V* is a volume of computational cell given by following equation for 3-dimensional cases.

$$
\Delta V = \Delta \mathbf{x}\_1 \cdot \Delta \mathbf{x}\_2 \cdot \Delta \mathbf{x}\_3 \tag{11}
$$

To satisfy Eq. (10), the segment *b* is adjusted as fraction *Vl*/*V* agrees with volumetric fraction of fluid. Between the volume of fluid *Vl* and segment *b*, there is the relation that is expressed as the following equation for three dimensional cases.

$$V\_m = \frac{1}{6a\_1a\_2a\_3} \left| b^3 - \sum\_{i=1}^3 \left[ \max\left( b - a\_i \Delta x\_i, 0 \right) \right]^3 - \sum\_{i=1}^3 \left[ \max\left( b - b\_m + a\_i \Delta x\_i, 0 \right) \right]^3 \right| \tag{12}$$

In this equation, *bm* is the maximum value that *b* can take:

$$b\_m = \sum\_{i=1}^{3} a\_i \Delta \mathbf{x}\_i \tag{13}$$

In this study, the Newton's method is used to estimate the segment *b* that satisfies Eq. (12). In the next step (transportation step), polyhedrons within each control volume are transported along a surrounding velocity field (Fig.1 (b)). Some parts of the volume of this polyhedron are transported to surrounding control volume. For example, in Fig.1(b), *V* is transferred from control volume 1 to control volume 2. Each transferred volume of liquid is calculated in this step. And each transferred volume of gas is also calculated by a procedure same as the volume of liquid.

Fig. 1. (a) Reconstruction step

In the equation, *nm* represents dimension of simulation, and is 2 or 3. *xi* is the coordinate position of the definition point of scalar quantities. A unit normal vector to the interface **a**=

To evaluate the function *f*(**x**), we assumed that the interface exists in the position where volume fraction of liquid *fl* equals to 0.5. By least-squares method (choose eight nearest neighbors for 2-dimensional case, 26 nearest neighbors for 3-dimensional case), linear

*<sup>l</sup> F f ax b*

In the equation, *b*0 is segment of linear function evaluated by least square method. Control volume is divided into two polyhedrons by the linear functions. If gray polyhedron in Fig.2 corresponds to liquid region in control volume, volume of this polyhedron: *Vl* must be

where *V* is a volume of computational cell given by following equation for 3-dimensional

To satisfy Eq. (10), the segment *b* is adjusted as fraction *Vl*/*V* agrees with volumetric fraction of fluid. Between the volume of fluid *Vl* and segment *b*, there is the relation that is

> 6 *m i <sup>i</sup> <sup>m</sup> <sup>i</sup> <sup>i</sup> i i V b bax bb ax*

> > 3

1 *m ii i b ax* 

transferred from control volume 1 to control volume 2. Each transferred volume of liquid is calculated in this step. And each transferred volume of gas is also calculated by a procedure

Fig. 1. (a) Reconstruction step

In this study, the Newton's method is used to estimate the segment *b* that satisfies Eq. (12). In the next step (transportation step), polyhedrons within each control volume are transported along a surrounding velocity field (Fig.1 (b)). Some parts of the volume of this polyhedron are transported to surrounding control volume. For example, in Fig.1(b),

3 3 3 3 <sup>3</sup>

max , 0 max , 0

(12)

1 <sup>0</sup> 0.5 *mn i i*

**<sup>x</sup>** (9)

*V fV l l* (10)

<sup>123</sup> *V xxx* · · (11)

. (13)

*V* is

*i*

(*a*1,*a*2,*a*3) and segment *b* in Eq.(8) must be estimated.

equals to proper volume of liquid. Thus,

1

same as the volume of liquid.

expressed as the following equation for three dimensional cases.

In this equation, *bm* is the maximum value that *b* can take:

1 1 123

*aaa*

function is obtained.

cases.

(c) Redistribution step Fig. 1. Schematic drawings of volume fraction transport calculation using the advanced

interface tracking method (2-dimensional case).

In transportation step, transfer volume of both phases is calculated. In the last step (redistribution step), as shown in Fig.1(c), transferred mass of both phases can be evaluated by multiplying density by transferred volume. Mass at new time step can be calculated by summing these transferred mass.

Fig. 2. Approximated fluid segment in control volume (3-dimensional case).

Development of Two-Phase Flow Correlation

(Rayleigh, 1879).

*l*=2×10-5 m2/s, density,

results. In the table 1,

the liquid.

for Fluid Mixing Phenomena in Boiling Water Reactor 293

were used, one is a fine grid (41×41×41, case A) and another is a course grid (9×9×9, case B). Time changes of diameters of the droplet are shown in Fig.4. The vibration cycle is about 0.022 seconds and agrees with the theoretical value that expressed as following equation

3

(14)

*ave\_cal*

*<sup>N</sup>* represents Nusselt's mean film thickness (Nusselt, 1916), and is

1/3

*ave\_exp* =2.06×10-2 N/m. And air

(15)

*N*

2 *l r*

The results of the case A and those of the case B are almost the same, and the effect of the different grid number is small. The diffusion at the gas-liquid interface was not observed.

In the fluid mixing phenomena, bubble, bubble/slug, slug and slug/churn flow are important. However, in slug flow, liquid film is observed between slug bubble and wall.

The TPFIT code was applied to numerical simulation of liquid film falling down on inclined flat plate. Simulations were performed with the same conditions as the experiment by Moran et al. (2003) (see Fig.5). Physical properties of the liquid were as follows: kinematic viscosity,

*l*=960 kg/m3, and surface tension,

properties at 300K and atmospheric pressure were used as gas properties. On all walls, nonslip boundary condition was assigned, and inlet pressure was fixed at atmospheric pressure. The flow conditions were summarized in Table 1. The analysis conditions were set up to compare the probability density function (PDF) of local film thickness with the experimental

3 *l*

In this equation, *gz* is flow direction acceleration by gravity force, and *J* is mass flow rate of

1 0.333 13 0.91 0.85 0.84 2 2.55 106 1.73 1.67 1.66 3 5.45 220 2.31 2.15 2.14

*z*

*J g* 

*N*

Film Reynolds number

 

Then we confirmed the effectiveness of the advanced interface tracking method.

**3.1.2 Liquid film falling down on inclined flat plate** 

evaluated by the following equation:

Case Inlet flow rate *J* (l/min)

Fig. 5. Analytical geometry of a liquid film.

Table 1. Numerical conditions.

Then verification of the TPFIT must be performed for film flow.
