**2.4.3 Collision between bubbles**

In a general VOF method, two interfaces inside the same grid cell cannot be distinguished. Thus, coalescence occurs when two bubbles are very close to each other. To avoid this type of coalescence, separate VOF functions are assigned to the bubbles. A repulsive force is applied when two bubbles come very close to each other to avoid the overlap of the bubbles. Suppose Bubble 1 and Bubble 2 approach each other and their VOF functions are given by *F*1 and *F*<sup>2</sup> . We make use of their smoothed functions *F*<sup>1</sup> and *F*<sup>2</sup> which are obtained by the convolution using the kernel:

Numerical Study on Flow Structures and

Table 1. Computational conditions for a rising bubble.

Fig. 2. Time evolution of rise velocities for a single bubble.

spurious current is as small as theirs.

simulations are performed on the domain of

flow was set so that the friction Reynolds number Re

**2.6 Computational conditions** 

Heat Transfer Characteristics of Turbulent Bubbly Upflow in a Vertical Channel 125

Eötvös number 1.0 Archimedes number 900 Density ratio 0.02 Viscosity ratio 0.02 Void fraction 0.06

Fig. 2 shows the time evolution of the bubble rise velocity for the three cases. The horizontal straight line of 74.31Re *<sup>d</sup>* represents the terminal velocity in the case of 9.38 *xd*0 in Bunner & Tryggvason (2002). The tendency of convergence is clearly seen as the grid

It was found by an oscillation test that more than 24 grid points per bubble diameter are needed to simulate oscillation motions of bubbles. A static drop test was also conducted under the same condition as in Francois et al. (2006) to find that the magnitude of the

Since many grid points are needed to resolve the fluid motions around and inside the bubbles, simulations are limited to low Reynolds numbers with small computational domains. We utilize so-called 'minimal turbulent channel' in the present study. The

> )2/(2

In the simulation of Lu & Tryggvason (2008), the constant pressure gradient to drive the

computational domain was 2004.254400 in wall units, and the domain was sufficiently large to sustain turbulence, as was shown in Jimenez & Moin (1991). The resulting channel Reynolds number was 3786 in their simulation of the single-phase turbulent flow. In their flow laden with nearly spherical bubbles, the channel Reynolds number was reduced to less than 2000, which may be too low to examine turbulence statistics of a bubbly flow. In the present study, the channel Reynolds number is set at 3786, and the volume flow rate is kept

as in Lu & Tryggvason (2008).

was 127.2. Therefore, the size of the

resolution is increased. This indicates the validity of our numerical scheme.

$$K(r,\gamma) = \begin{cases} 1 - 6\left(r/\gamma\right)^2 + 6\left(r/\gamma\right)^3 & \text{(r/\gamma < 1/2)}\\ 2\left(r/\gamma\right)^3 & \text{(1/2 \le r/\gamma < 1)}\\ 0 & \text{(1 \le r/\gamma)} \end{cases} \tag{19}$$

The repulsive force

$$\mathbf{f}\_{\mathbb{R}1} = -CF\_1 \nabla \bar{F}\_{2'} \tag{20}$$

is exerted on Bubble 1, where *C* is a positive constant. Similarly, the force

$$\mathbf{f}\_{R2} = -CF\_2\nabla \tilde{F}\_{1'} \tag{21}$$

is exerted on Bubble 2. The integral of **f f** *R R* 1 2 over the whole computational domain is exactly zero. Therefore, this repulsive force does not affect the total momentum in the system. In the present study, we set 3 4 *x x* , where *x* is the grid spacing. Smoothed VOF functions are also utilized for the estimation of the fluid and thermal properties in Eqs.(5), (6), (9) and (10) with 2*x* (Lorstad & Fuchs, 2004).

#### **2.4.4 Time integration and spatial discretization**

The momentum and energy equations are solved on a collocated grid by using the finite difference schemes. All of the velocity components and the pressure are defined at the center of the grid cell. The time-integration is based on a fractional-step method where a pseudo-pressure is used to correct the velocity field so that the continuity equation is satisfied. A balanced force algorithm developed by Francois et al. (2006) is used to suppress unphysical flows (or spurious currents) resulting from the unbalance between the interfacial tension and the pressure difference across the interface. Poisson's equation for the pseudopressure is solved by using a multigrid method.

The QUICK method (Leonard, 1979) is applied in the finite differencing of the convection terms of momentum and energy equations. The second-order central difference scheme is applied for the finite differencing of the viscous terms of the momentum equations and the diffusion term of the energy equation. The 2nd-order improved Euler method is used for the time-integration of the convective and viscous (or diffusion) terms (Rudman, 1998). The velocity field at /2 *<sup>n</sup> t t* is used to advect the VOF function.

#### **2.5 Validation of numerical scheme**

We compared the rise velocity of a bubble in otherwise quiescent fluid with that in Bunner & Tryggvason (2002). The bubble rise velocity can be described by the bubble Reynolds number, Re*db c u d*<sup>0</sup> , where *ub* and 0 *d* denote the rise velocity and the bubble diameter, respectively. Numerical parameters were the same as that in Bunner & Tryggvason (2002). The computational conditions are summarized in Table 1. There cases with different grid resolutions were examined with 32 32 32 , 64 64 64 , and 128 128 128 grid points. For the three cases, the bubble diameter equals to 15.5*x* , 31.1*x* , and 62.2*x* , respectively. Here, *x* denotes the grid spacing. Periodic boundary conditions were imposed in all directions. Domain size was set to be 2 2 2 .


Table 1. Computational conditions for a rising bubble.

124 Computational Simulations and Applications

2 3 3 16 6 1 2

is exerted on Bubble 2. The integral of **f f** *R R* 1 2 over the whole computational domain is exactly zero. Therefore, this repulsive force does not affect the total momentum in the

Smoothed VOF functions are also utilized for the estimation of the fluid and thermal

The momentum and energy equations are solved on a collocated grid by using the finite difference schemes. All of the velocity components and the pressure are defined at the center of the grid cell. The time-integration is based on a fractional-step method where a pseudo-pressure is used to correct the velocity field so that the continuity equation is satisfied. A balanced force algorithm developed by Francois et al. (2006) is used to suppress unphysical flows (or spurious currents) resulting from the unbalance between the interfacial tension and the pressure difference across the interface. Poisson's equation for the pseudo-

The QUICK method (Leonard, 1979) is applied in the finite differencing of the convection terms of momentum and energy equations. The second-order central difference scheme is applied for the finite differencing of the viscous terms of the momentum equations and the diffusion term of the energy equation. The 2nd-order improved Euler method is used for the time-integration of the convective and viscous (or diffusion) terms (Rudman, 1998). The

We compared the rise velocity of a bubble in otherwise quiescent fluid with that in Bunner & Tryggvason (2002). The bubble rise velocity can be described by the bubble Reynolds

respectively. Numerical parameters were the same as that in Bunner & Tryggvason (2002). The computational conditions are summarized in Table 1. There cases with different grid resolutions were examined with 32 32 32 , 64 64 64 , and 128 128 128 grid points. For the three cases, the bubble diameter equals to 15.5*x* , 31.1*x* , and 62.2*x* , respectively. Here, *x* denotes the grid spacing. Periodic boundary conditions were

, where *ub* and 0 *d* denote the rise velocity and the bubble diameter,

 .

*t t* is used to advect the VOF function.

 

*rrr*

(, ) 2 12 1

*K r r r*

is exerted on Bubble 1, where *C* is a positive constant. Similarly, the force

The repulsive force

system. In the present study, we set

properties in Eqs.(5), (6), (9) and (10) with

**2.4.4 Time integration and spatial discretization** 

pressure is solved by using a multigrid method.

velocity field at /2 *<sup>n</sup>*

number, Re*db c u d*<sup>0</sup>

**2.5 Validation of numerical scheme** 

imposed in all directions. Domain size was set to be 2 2 2

0 1.

2*x* (Lorstad & Fuchs, 2004).

*r*

 

1 12 , **<sup>f</sup>***<sup>R</sup> CF F* (20)

2 21 , **<sup>f</sup>***<sup>R</sup> CF F* (21)

3 4 *x x* , where *x* is the grid spacing.

(19)

Fig. 2 shows the time evolution of the bubble rise velocity for the three cases. The horizontal straight line of 74.31Re *<sup>d</sup>* represents the terminal velocity in the case of 9.38 *xd*0 in Bunner & Tryggvason (2002). The tendency of convergence is clearly seen as the grid resolution is increased. This indicates the validity of our numerical scheme.

Fig. 2. Time evolution of rise velocities for a single bubble.

It was found by an oscillation test that more than 24 grid points per bubble diameter are needed to simulate oscillation motions of bubbles. A static drop test was also conducted under the same condition as in Francois et al. (2006) to find that the magnitude of the spurious current is as small as theirs.
