**2. Direct numerical simulations of bubbly flows**

## **2.1 Configuration**

Figure 1 shows the spatial configuration considered in the present study. The *x*, *y* and *z* axes are assigned to the streamwise, wall-normal and spanwise directions, respectively. The gravitational force is assumed to point to the negative *x* direction. Here, we consider the system where the flow is heated with a uniform heat flux from both the walls and the mean temperature increases linearly in the streamwise (or vertical) direction. The wall heat flux, *qW* , is determined so that the energy (enthalpy) of the system is conserved. In this situation, the wall heat flux is kept nearly constant for stationary turbulence.

Fig. 1. Configuration.

It is assumed that the fluids are incompressible for both of the liquid (or continuous) and gas (or dispersed) phases in the present study. The buoyancy force originated from the nonuniform spatial distribution of temperature is assumed to be negligibly small. It is also assumed that the fluid properties do not depend on the temperature.

#### **2.2 Basic equations**

In this study, simulations of the bubbly flow are conducted by using VOF (Volume of Fluid) method. In the VOF mothod, the fraction of the bubble gas, *F*, occupying each computational grid cell is advected by using the equation

$$\frac{\partial F}{\partial t} = -\mu\_i \frac{\partial F}{\partial \mathbf{x}\_i} \tag{1}$$

where *x* = (*x*1,*x*2,*x*3) = (*x*,*y*,*z*) denotes the position, and *u* = (*u*1,*u*2,*u*3) = (*u*,*v*,*w*) represents the fluid velocity. In the present study, summation convention is applied to repeated subscripts if not otherwise specified.

The motions of incompressible fluid are governed by the Navier-Stokes (momentum) equations

$$\rho \left( \frac{\partial \mathbf{u}\_i}{\partial t} + \mathbf{u}\_k \frac{\partial \mathbf{u}\_i}{\partial \mathbf{x}\_k} \right) = -\frac{\partial (p - P)}{\partial \mathbf{x}\_i} + \frac{\partial \tau\_{ij}}{\partial \mathbf{x}\_j} + f\_{\sigma i} - \left( \rho - \{\rho\} \right) \mathbf{g} \delta\_{i1} + \rho \delta\_{i1'} \tag{2}$$

supplemented with the continuity equation

$$\frac{\partial \mathbf{u}\_k}{\partial \mathbf{x}\_k} = \mathbf{0},\tag{3}$$

where

120 Computational Simulations and Applications

Figure 1 shows the spatial configuration considered in the present study. The *x*, *y* and *z* axes are assigned to the streamwise, wall-normal and spanwise directions, respectively. The gravitational force is assumed to point to the negative *x* direction. Here, we consider the system where the flow is heated with a uniform heat flux from both the walls and the mean temperature increases linearly in the streamwise (or vertical) direction. The wall heat flux, *qW* , is determined so that the energy (enthalpy) of the system is conserved. In this situation,

It is assumed that the fluids are incompressible for both of the liquid (or continuous) and gas (or dispersed) phases in the present study. The buoyancy force originated from the nonuniform spatial distribution of temperature is assumed to be negligibly small. It is also

In this study, simulations of the bubbly flow are conducted by using VOF (Volume of Fluid) method. In the VOF mothod, the fraction of the bubble gas, *F*, occupying each

> *F F u t x*

where *x* = (*x*1,*x*2,*x*3) = (*x*,*y*,*z*) denotes the position, and *u* = (*u*1,*u*2,*u*3) = (*u*,*v*,*w*) represents the fluid velocity. In the present study, summation convention is applied to repeated subscripts

The motions of incompressible fluid are governed by the Navier-Stokes (momentum)

*k ij*

*u u p P <sup>u</sup> f g tx xx*

 

( ) , *ij i i k i ii*

 

, *<sup>i</sup> i*

(1)

1 1

(2)

**2. Direct numerical simulations of bubbly flows** 

the wall heat flux is kept nearly constant for stationary turbulence.

assumed that the fluid properties do not depend on the temperature.

computational grid cell is advected by using the equation

**2.1 Configuration** 

Fig. 1. Configuration.

**2.2 Basic equations** 

if not otherwise specified.

equations

$$
\sigma \pi\_{ij} = \mu \left( \frac{\partial u\_j}{\partial \mathbf{x}\_i} + \frac{\partial u\_i}{\partial \mathbf{x}\_j} \right) \tag{4}
$$

is the viscous stress. Here, and denote the density and viscosity of the fluid. Their values for each grid cell are respectively given by the simple average,

$$
\rho = F\rho\_d + (1 - F)\rho\_{c'} \tag{5}
$$

and the harmonic average,

$$\frac{1}{\mu} = F \frac{1}{\mu\_d} + (1 - F) \frac{1}{\mu\_c} \tag{6}$$

of the two fluids. The subscripts *d* and *c* denote the dispersed and continuous phases, respectively. *p* is the pressure and *P* denotes the mean pressure linearly decreasing in the vertical direction. **f** represents the volume force associated with the interfacial surface tension of the bubbles, and *g* denotes the gravitational acceleration. represents the spatial average over the whole computational domain. The last term in Eq.(2) represents the total driving force exerted on the fluids flowing through the channel, and

$$
\beta = -\frac{dP}{d\mathbf{x}} - \langle \rho \rangle \mathbf{g}.\tag{7}
$$

For stationary turbulence, the temporal average of equals to *W* , where *<sup>W</sup>* represents the wall shear stress and denotes the channel half width.

Since the mean temperature increases linearly downstream, the temperature *T* is decomposed as *T Gx* , where *G* denotes the mean temperature gradient in the streamwise direction and represents the temperature variance. The governing equation for (or energy equation) is described by

$$\left(\rho \mathbf{C}\_{p}\right)\left(\frac{\partial \Theta}{\partial t} + \mu\_{k}\frac{\partial \Theta}{\partial \mathbf{x}\_{k}} - \mathbf{G}\mu\_{1}\right) = \frac{\partial}{\partial \mathbf{x}\_{j}}\left(k\frac{\partial \Theta}{\partial \mathbf{x}\_{j}}\right)\tag{8}$$

where

$$
\rho \mathbf{C}\_P = \mathbf{F} \rho\_d \mathbf{C}\_{Pd} + (\mathbf{1} - \mathbf{F}) \rho\_c \mathbf{C}\_{Pc} \tag{9}
$$

and

$$\frac{1}{k} = F \frac{1}{k\_d} + (1 - F) \frac{1}{k\_c} \tag{10}$$

represent the volume-averaged heat capacity (per unit volume) and heat conductivity, respectively. Here, *CP* denotes specific heat.

Numerical Study on Flow Structures and

**2.4 Numerical methods 2.4.1 PLIC-VOF algorithm** 

and *z* directions.

where 

**2.4.2 Interfacial tension** 

the interface, and *<sup>S</sup>*

continuum limit. The curvature,

a static drop (Francois et al., 2006).

**2.4.3 Collision between bubbles** 

convolution using the kernel:

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

where *W* is the averaged temperature variance at the walls. The enhancement of heat

<sup>2</sup> 2Re Pr . *<sup>t</sup> <sup>c</sup> c m W*

In the present study, we use a PLIC-VOF (Piecewise Linear Interface Calculation - Volume of Fluid) method to capture the interface of the bubble. The position of the interface is determined by the volume fraction *F* of the bubble. *F* 1 represents a cell filled with the gas (dispersed-phase fluid) of the bubble, while 0 *F* indicates that the cell is filled with the continuous-phase liquid. The cells of 0 1 *F* include the interface. The time evolution of *F* is estimated with a PLIC-VOF algorithm (Rudman; 1998, Scardovelli & Zaleski, 2003). In PLIC-VOF algorithm, we take account of the slope of the interface. Young and Parkar's method is used to estimate the normal vectors of the interface from the values of *F* in adjacent cells (Parker & Youngs, 1992). A mass-conserving scheme which also satisfies the consistency condition ( 0 1 *F* ) is applied for the advection of the VOF function. The EI-LE (Eulerian implicit Lagrangian explicit) scheme (Scardovelli & Zaleski, 2003; Aulisa et al., 2003; Aulisa et al., 2007), which is an advection scheme originally designed for twodimensional incompressible velocity field, is applied for the three-dimensional (3D) incompressible velocity field. The extension is conducted by decomposing the 3D velocity field into three two-component velocity fields by the use of Fourier transformation in the *x*

In the continuum surface force (CSF) method, the interfacial tension force is calculated as

height functions (Lorstad & Fuchs, 2004), which effectively eliminates spurious currents for

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

 

is the coefficient of the interfacial tension,

force in Eq.(2) is calculated by using Eq.(18) and the relation **n**

*F*1 and *F*<sup>2</sup> . We make use of their smoothed functions *F*<sup>1</sup>

,

is a delta function concentrated on the interface. The interfacial tension

*<sup>S</sup>* **f n** (18)

, is calculated with a high degree of accuracy by using

and *F*<sup>2</sup>

is the curvature, **n** is the normal to

*<sup>S</sup> F* , which holds in the

which are obtained by the

 (17)

transfer can be estimated by the Nusselt number which is defined by

*<sup>h</sup> Nu k*

#### **2.3 Control parameters**

The parameters which control the system considered in the present study are the channel Reynolds number Re*<sup>m</sup>* , Eötvös number *Eo* , Morton number *M* , Archimedes number *Ar* , dispersed-phase Prandtl number Pr*d dd a* , continuous-phase Prandtl number Pr*c cc a* , the density ratio *d c* , the viscosity ratio *d c* , and the ratio of specific heat, *C C Pd Pc* , where ( ) *d d d Pd ak C* and ( ) *c c c Pc ak C* are thermal diffusivities in the dispersed and continuous phases, respectively. First four parameters are defined as

$$\mathrm{Re}\_m = \frac{2\,\mathrm{U}\_m \mathcal{S}}{\nu\_c}, \mathrm{E}\sigma = \frac{\Delta\rho \,\mathrm{g}\,\mathrm{d}\_0^2}{\sigma}, \mathrm{M} = \frac{\mathrm{g}\,\mu\_c^4 \Delta\rho}{\rho\_c^2 \sigma^3}, \mathrm{Ar} = \frac{\rho\_c \Delta\rho \,\mathrm{g}\,\mathrm{d}\_0^3}{\mu\_c^2}, \tag{11}$$

where *U u <sup>m</sup>* , , 0 *d* , and denote the mean velocity, kinematic viscosity, bubble diameter, and the surface tension coefficient, respectively. *c d* represents the density difference between the two phases.

The wall units are used to normalize the time, length, velocity and temperature. The friction velocity is defined by *u W c* . The friction time and the friction length are defined by 2 *<sup>c</sup> t u* and *<sup>c</sup> l u* , respectively. The friction Reynolds number is defined as

$$\text{Re}\_{\tau} = \frac{\mu\_{\tau} \delta}{\nu\_{c}} = \frac{\delta}{l\_{\tau}} = \frac{\text{Re}\_{m}}{2\text{LI}\_{m}^{+}} \text{ } \tag{12}$$

where the superscript '+' represents the values normalized by the wall units. The temperature is normalized by the friction temperature defined by

$$\Theta\_{\tau} = \frac{q\_W}{\rho\_c \mathbb{C}\_{Pc} u\_{\tau}} \, ^\prime \tag{13}$$

where *qW* denotes the wall heat flux.

Eötvös number, *Eo* , represents the ratio of the buoyancy and surface-tension forces, while Archimedes number, *Ar* , represents the ratio of the buoyancy and viscous forces. Morton number, *M* , is not independent from Eötvös number *Eo* and Archimedes number *Ar* and is expressed as 2 2 *M Eo Ar* . The 'buoyancy parameter', which is defined by

$$B\mu = \frac{\beta}{g\Delta\rho} = \left(\frac{\text{Re}\_r^2}{Ar}\right) \left(\frac{d\_0}{\delta}\right)^3 \,\text{\,\,\,\,\tag{14}$$

is also an important parameter. The buoyancy effect of bubbles may be large when this parameter is small.

The mean temperature

$$\Theta\_m = \frac{\left\langle (\rho \mathbb{C}\_P) \Theta \,\mu \right\rangle}{\left\langle (\rho \mathbb{C}\_P) \mu \right\rangle},\tag{15}$$

is used as the representative temperature of the fluid in the present study. The heat transfer coefficient can be defined as

$$\eta h\_t = \frac{q\_W}{\Theta\_m - \Theta\_W}'\tag{16}$$

where *W* is the averaged temperature variance at the walls. The enhancement of heat transfer can be estimated by the Nusselt number which is defined by

$$N\mu = \frac{2\,\delta h\_t}{k\_c} = \frac{2\,\text{Re}\_\tau\text{Pr}\_c}{\Theta\_m^+ - \Theta\_W^+}.\tag{17}$$

#### **2.4 Numerical methods 2.4.1 PLIC-VOF algorithm**

122 Computational Simulations and Applications

The parameters which control the system considered in the present study are the channel Reynolds number Re*<sup>m</sup>* , Eötvös number *Eo* , Morton number *M* , Archimedes number *Ar* ,

and ( ) *c c c Pc ak C*

<sup>2</sup> Re , , , , *m c <sup>c</sup> <sup>m</sup> c c c <sup>U</sup> gd g gd Eo M Ar*

The wall units are used to normalize the time, length, velocity and temperature. The friction

Re Re , <sup>2</sup>

where the superscript '+' represents the values normalized by the wall units. The

Eötvös number, *Eo* , represents the ratio of the buoyancy and surface-tension forces, while Archimedes number, *Ar* , represents the ratio of the buoyancy and viscous forces. Morton number, *M* , is not independent from Eötvös number *Eo* and Archimedes number *Ar* and

> <sup>2</sup> <sup>3</sup> Re <sup>0</sup> , *<sup>d</sup> Bu g Ar*

is also an important parameter. The buoyancy effect of bubbles may be large when this

is used as the representative temperature of the fluid in the present study. The heat transfer

, *<sup>W</sup> <sup>t</sup> m W*

*m*

( ) , ( ) *P*

*P C u C u*

 

*u*

 

 

is expressed as 2 2 *M Eo Ar* . The 'buoyancy parameter', which is defined by

*c m*

*l U*

, *<sup>W</sup> c Pc q C u*

 

, respectively. The friction Reynolds number is defined as

*m*

*d c* 

24 3 0 0 2 3 2

(11)

 

denote the mean velocity, kinematic viscosity, bubble

 

 

The friction time and the friction length are defined by

(12)

(13)

(15)

*<sup>q</sup> <sup>h</sup>* (16)

, the viscosity ratio

dispersed and continuous phases, respectively. First four parameters are defined as

 

> 

*a* , continuous-phase Prandtl number

, and the ratio of specific heat,

*c d* represents the

(14)

are thermal diffusivities in the

**2.3 Control parameters** 

*a* , the density ratio

*C C Pd Pc* , where ( ) *d d d Pd ak C*

velocity is defined by *u W c* .

where *qW* denotes the wall heat flux.

 and *<sup>c</sup> l u* 

density difference between the two phases.

Pr*c cc* 

where *U u <sup>m</sup>* ,

2

 

parameter is small. The mean temperature

coefficient can be defined as

*<sup>c</sup> t u* 

dispersed-phase Prandtl number Pr*d dd*

*d c* 

diameter, and the surface tension coefficient, respectively.

temperature is normalized by the friction temperature defined by

, 0 *d* , and

  In the present study, we use a PLIC-VOF (Piecewise Linear Interface Calculation - Volume of Fluid) method to capture the interface of the bubble. The position of the interface is determined by the volume fraction *F* of the bubble. *F* 1 represents a cell filled with the gas (dispersed-phase fluid) of the bubble, while 0 *F* indicates that the cell is filled with the continuous-phase liquid. The cells of 0 1 *F* include the interface. The time evolution of *F* is estimated with a PLIC-VOF algorithm (Rudman; 1998, Scardovelli & Zaleski, 2003). In PLIC-VOF algorithm, we take account of the slope of the interface. Young and Parkar's method is used to estimate the normal vectors of the interface from the values of *F* in adjacent cells (Parker & Youngs, 1992). A mass-conserving scheme which also satisfies the consistency condition ( 0 1 *F* ) is applied for the advection of the VOF function. The EI-LE (Eulerian implicit Lagrangian explicit) scheme (Scardovelli & Zaleski, 2003; Aulisa et al., 2003; Aulisa et al., 2007), which is an advection scheme originally designed for twodimensional incompressible velocity field, is applied for the three-dimensional (3D) incompressible velocity field. The extension is conducted by decomposing the 3D velocity field into three two-component velocity fields by the use of Fourier transformation in the *x* and *z* directions.
