Numerical Study on Flow Structures and Heat Transfer Characteristics of Turbulent Bubbly Upflow in a Vertical Channel

Turbulent bubbly flows have attracted a lot of attention because of their importance for many practical applications such as flows in chemical plants and nuclear power plants. Enhancement of heat transfer by bubble-induced turbulence also attracts a lot of attention from the view point of energy saving. Many studies have been conducted for the motion of bubbles and the characteristics of heat-transfer in turbulent bubbly flows. It is expected that the research on turbulent bubbly flows is accelerated by fully resolved simulations of bubble-turbulence interaction (Tryggvason et al., 2011). The characteristics of bubbly upflow strongly depend on the motions of bubbles and resulting void fraction distribution in the flow. Serizawa et al. (1975) found that the local void fraction is high near the walls and is lower in the core region of upflow in a pipe. Liu (1993) also found in the experiments of turbulent bubbly upflow in a vertical channel that the void fraction has peaks near the walls for the bubbles smaller than 5-6mm, while it has a peak in the core of the channel for the bubbles larger than 5-6mm. Lu & Tryggvason (2008) also showed in their direct numerical simulations of turbulent bubbly upflow in a vertical channel that nearly spherical bubbles tend to concentrate on the near-wall regions, while strongly deformable bubbles tend to be expelled from the near-wall regions. They also showed that the turbulence structures are changed by the motions of bubbles. The detailed mechanism of turbulence modulation due to the bubbles, however, has not been fully clarified yet. Some experimental studies have been conducted for heat-transfer enhancement by the injection of bubbles. Tamari & Nishikawa (1976) showed in their experiments of laminar natural convection heat transfer in water from a vertical plate that the heat transfer is enhanced by the injection of air bubbles. The enhancement of heat transfer by bubble injection was studied further in detail by Tokuhiro & Lykoudis (1994) and Kitagawa et al. (2008, 2010). However, the mechanism of the heat-transfer enhancement has not yet been fully clarified especially in turbulent flows. In the present study, direct numerical simulations have been conducted for turbulent bubbly upflow between two parallel heating walls in order to clarify its heat transfer characteristics. The mechanism of the heat-transfer enhancement is examined by performing simulations with different values of control parameters. The performance of the heat-transfer enhancement is discussed based on the numerical results.


Introduction
Turbulent bubbly flows have attracted a lot of attention because of their importance for many practical applications such as flows in chemical plants and nuclear power plants. Enhancement of heat transfer by bubble-induced turbulence also attracts a lot of attention from the view point of energy saving. Many studies have been conducted for the motion of bubbles and the characteristics of heat-transfer in turbulent bubbly flows. It is expected that the research on turbulent bubbly flows is accelerated by fully resolved simulations of bubble-turbulence interaction (Tryggvason et al., 2011). The characteristics of bubbly upflow strongly depend on the motions of bubbles and resulting void fraction distribution in the flow. Serizawa et al. (1975) found that the local void fraction is high near the walls and is lower in the core region of upflow in a pipe. Liu (1993) also found in the experiments of turbulent bubbly upflow in a vertical channel that the void fraction has peaks near the walls for the bubbles smaller than 5-6mm, while it has a peak in the core of the channel for the bubbles larger than 5-6mm. Lu & Tryggvason (2008) also showed in their direct numerical simulations of turbulent bubbly upflow in a vertical channel that nearly spherical bubbles tend to concentrate on the near-wall regions, while strongly deformable bubbles tend to be expelled from the near-wall regions. They also showed that the turbulence structures are changed by the motions of bubbles. The detailed mechanism of turbulence modulation due to the bubbles, however, has not been fully clarified yet. Some experimental studies have been conducted for heat-transfer enhancement by the injection of bubbles. Tamari & Nishikawa (1976) showed in their experiments of laminar natural convection heat transfer in water from a vertical plate that the heat transfer is enhanced by the injection of air bubbles. The enhancement of heat transfer by bubble injection was studied further in detail by Tokuhiro & Lykoudis (1994) and Kitagawa et al. (2008Kitagawa et al. ( , 2010. However, the mechanism of the heat-transfer enhancement has not yet been fully clarified especially in turbulent flows. In the present study, direct numerical simulations have been conducted for turbulent bubbly upflow between two parallel heating walls in order to clarify its heat transfer characteristics. The mechanism of the heat-transfer enhancement is examined by performing simulations with different values of control parameters. The performance of the heat-transfer enhancement is discussed based on the numerical results. 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, W q , 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. 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.

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 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.
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 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. 1 F  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.

Interfacial tension
In the continuum surface force (CSF) method, the interfacial tension force is calculated as where  is the coefficient of the interfacial tension,  is the curvature, n is the normal to the interface, and S  is a delta function concentrated on the interface. The interfacial tension force in Eq.
(2) is calculated by using Eq.(18) and the relation S F    n , which holds in the continuum limit. The curvature,  , is calculated with a high degree of accuracy by using height functions (Lorstad & Fuchs, 2004), which effectively eliminates spurious currents for a static drop (Francois et al., 2006).

Collision between bubbles
In a general VOF method, two interfaces inside the same grid cell cannot be distinguished.
is exerted on Bubble 2. The integral of 12 RR  ff 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 34 xx    , where x  is the grid spacing.

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 n tt  is used to advect the VOF function.

Validation of numerical scheme
We compared the rise velocity of a bubble in otherwise quiescent fluid with that in Bunner & Tryggvason (2002 . 38 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. 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.

Computational conditions
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 simulations are performed on the domain of Lu & Tryggvason (2008). In the simulation of Lu & Tryggvason (2008), the constant pressure gradient to drive the flow was set so that the friction Reynolds number 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 www.intechopen.com constant. Notice that the friction velocity (and the friction Reynolds number) is generally changed by the effects of the bubbles. Non-slip boundary conditions are imposed in the wall-normal direction for the velocity components. Periodic boundary conditions are imposed in the x and z directions for the velocity, the pressure variance, p P  , and the temperature variance, .  As mentioned above, we assume a constant temperature gradient in the vertical (streamwise) direction. We impose a uniform heat flux from both walls. In the present study, the energy (enthalpy) of the system is kept constant, so that the instantaneous wall heat flux W q is given by Here,  represents the spatial average in the x and z directions.

Bubbly flow
The simulations are performed with 256 256 128   rectangular grid cells. We set the fluid density inside the bubbles (density of the dispersed-phase fluid) to be one-tenth of that of , and we set the viscosities to be equal to reduce the computational cost as in Lu & Tryggvason (2008). Air bubbles with a diameter of 1-2mm in water are considered in the present study. Eötvös number, Morton number, and Archimedes number of the bubbles are 0.36, 10 2.91 10   , and 12700, respectively. These parameters correspond to a 1.64mm air bubble in the fluid whose viscosity is 1.84 times higher than that of the water at room temperature. Twelve bubbles with a diameter of 0.4 are introduced randomly into the turbulent singlephase flow in the channel. Computational conditions are summarized in Table 2. Although most of the parameters employed here are quite close to those in Lu & Tryggvason (2008), the buoyancy parameter, Bu, is considerably higher than their value of 0.018. This indicates that the buoyancy effects are less important in our simulation. In order to assess the importance of the buoyancy effects, we conduct a simulation for neutrally buoyant droplets, where Bu is infinite , as will be explained in 2.6.2. The thermal properties employed in the present simulation are summarized in Table. 3. The Prandtl number for the liquid (continuous-phase fluid) is set at a low value of 2.0 ( Pr 2.0 c  ) to maintain high numerical accuracy. The Prandtl number for the gas (dispersed-phase fluid) is also set at 2.0 ( Pr 2.0 d  ). In this case, the ratio of heat capacities per unit volume, dP d cP c CC   , is 0.1 and the ratio of thermal conductivity, dc kk , is 1.0, respectively. Hereafter, we call this run Case B1. For comparison, three cases (Cases B2-B4) with different thermal porperties are simulated. In Case B2, we change the Prandtl number of the continuous-phase fluid to 1.0 to examine Prandtl number dependence of the heat transfer characteristics of the turbulent bubbly flow. In Cases B3 and B4, we change the thermal properties inside the bubbles. In Case B3, the thermal conductivity inside the bubbles is set at 1/10 of that in the surrounding liquid in order to examine the heat insulating effect due to the bubbles. In Case B4, the specific heat of the gas inside the bubble is set to be 10 times larger ( 1.0 dP d cP c CC   ) to clarify the effect of lower heat capacity inside the bubbles. In order to check the accuracy of the simulation, we have conducted a simulation under the same physical conditions with a lower grid resolution. The parameters for this simulation are summarized in Table 4 Table 4. Computational conditions for the simulation with lower grid resolution.

Droplet flow
As mentioned above, we conduct a simulation for neutrally buoyant droplets in order to assess the importance of the buoyancy effects. The density ratio of the dispersed-phase fluid is changed to 1.0 from 0.1 in the bubbly flow. Computational conditions are summarized in www.intechopen.com In order to estimate the effects of the numerical diffusion caused by the QUICK scheme, a simulation was conducted by using the centered 2nd-order scheme (consistent scheme) for the convection terms in a staggered grid system (Kawamura et al., 1998). It is found that the amplitude of the streamwise component of vorticity is slightly (about 1.5%) lower in the www.intechopen.com case with the QUICK scheme than that with the centered 2nd-order scheme. The Nusselt number for higher Prandtl number of 2 Pr  c is also slightly (about 1%) lower. It can be concluded that the effects of the numerical diffusion due to the QUICK scheme are small.

Results and discussion
The turbulent bubbly (or droplet) flow and the temperature field reached a fully developed state about 1000   t after the injection of the bubbles (or the droplets). After the turbulence reached the fully developed state, the simulation has further been conducted for the period of 700 t  (or 1400 t  ) for the bubbly (or droplet) turbulent flow to obtain statistical quantities. The longer simulation for 2900 t  has been performed for the bubbly flow with the lower grid resolution.

Flow structures of turbulent bubbly flow 3.1.1 Wall shear stress and friction Reynolds number
In  Fig. 3(a) shows a typical snapshot of the bubbles and the vortical structures visualized by the second invariant of velocity gradient tensor,

Bubble motions and vortices
It is clearly seen that the bubbles tend to collect on the near-wall regions of the channel. The bubbles are slightly deformed form the spherical shape. As is shown in Fig. 3(b), the droplets are distributed rather uniformly throughout the channel though some droplets are located close to the walls as in the case of the bubbly flow.

www.intechopen.com
It is seen in the case of the bubbly flow that the vortices are locally activated by the bubbles in the near-wall regions. In the case of the droplet flow, strong vortices are observed even in the central region of the channel. Notice that the normalized values of Q are used for visualization in Figs. 3 and that the vortices are considerably strengthened due to the injection of the bubbles. Figs. 4(a) and 4(b) show the time evolution of the locations of the bubbles and the droplets, respectively. The bubbles rise along the walls almost all the time. As is shown in Fig. 4(c), the void fraction has sharp peaks near the walls as a result of bubble accumulation toward the walls. Horizontal motions are much more noticeable for the droplets. Some droplets rise along the walls as the bubbles, however. This is also confirmed by the profile of the volume fraction of the droplet in Fig. 4(c).  . The vorticities in the bubbly and droplet flows have high values compared with the single-phase flow, indicating that the generation of the streamwise vortices is enhanced by the injection of the bubbles or droplets. There are peaks at 0.07 y   and 1.93 near the walls in the profile of the bubbly flow. The corresponding peaks for the droplet flow are more distantly positioned from the walls, and the peaks in the single-phase flow are further away from the walls. In general, these peaks approach the walls as the friction Reynolds number is increased. We have conducted a simulation of the single-phase turbulent flow at a higher Re ( 160)   , which is comparable with that of the turbulent bubbly flow, to find that the peaks in the vorticity profile are located at 0.19 y   and 1.81. This indicates that the bubbles (or droplets) rising along the walls enhance the generation of quasi-streamwise vortices in the regions very close to the walls. This is verified by the visualization of the vortical structures around the bubbles (see Figs. 3). In the profile for the bubbly flow, there are small bumps at 0.3 y   and 1.7 in addition to the two peaks near the walls. As is shown in Fig. 6, these peaks and bumps correspond to www.intechopen.com the trailing vortices around the bubbles. Relatively small trailing vortices are seen on the wall side of the bubble (e.g. a vortex pair marked by A), while relatively large vortices are also seen on the side of the channel center (e.g. vortices marked by B). In Fig. 5(b), the profiles of the streamwise vorticity squared which is normalized by the friction time for each flow, 2 () x   , are shown. The vortices is relatively weak in the central region of the channel in the bubbly flow. Fig. 7(a) shows the mean velocity profiles for the bubbly flow. The black line represents the liquid velocity in the single-phase flow, and the red line and blue circles represent the liquid and gas velocities in the bubbly flow, respectively. Note that the gas velocity is not the rise velocity of the bubble centroids, but just the velocity of the gas inside the bubbles. Note also that the bubble diameter is 0.4 .  Since the wall shear stress is increased with the flow rate fixed, the mean velocity normalized by the friction velocity is reduced compared with that for the single-phase flow. The gas velocity is higher than the liquid velocity near the walls, while it is lower in the regions around 0.25 y   and 1.75. The bubbles are exposed to high shear near the walls. The balance between this shear stress and the interfacial surface tension leads to the higher gas velocity near the walls and the lower gas velocity on the central side of the channel. In fact, the streamwise velocity is homogenized by the circulating flow inside the bubble. Fig. 7(b) shows the mean velocity profiles for the droplet flow. The velocity for the dispersed-phase fluid is remarkably lower than that for the continuous-phase fluid in the regions around 0.25 y   and 1.75, indicating that the droplets are moving more slowly than the surrounding fluid. www.intechopen.com

Turbulence intensities
Figs. 8 show the profiles of velocity fluctuations of the liquid (continuous-phase fluid) for the single-phase, bubbly and droplet flows. In the case of the droplet flow, the streamwise component of the turbulence intensities decreases near the walls due to the droplets. The wall-normal and spanwise components, on the other hand, increase due to the presence of the droplets in the near-wall regions. The fluid motions normal to the interfaces of droplets are directly suppressed or induced by the presence of droplets. The droplets, on the other hand, indirectly affect the fluid motions. The turbulence is augmented by the droplets, and the redistribution mechanism of Reynolds stresses is enhanced. This may be one of the reasons for the decrease in the streamwise component and the increase in the wall-normal and spanwise components near the walls. The profiles of turbulence intensities in the bubbly flow resemble those in the droplet flow when they all are normalized by the friction velocity of the single-phase flow (figures not shown). When normalized by each value of the friction velocity, all components are considerably low in the turbulent bubbly flow compared with those in the droplet flow as is shown in Figs. 8. This is because the increase in the wall shear stress (and the friction velocity) is brought about by the factors other than turbulence augmentation, as is discussed in 3.1.5.

Shear-stress profiles
The wall-friction drag is increased due to the injection of the bubbles (Note that the volume flow rate is kept constant in the computation). Now, we examine the mechanism for the increase of the drag by considering the balance of forces in the channel. Taking the average of Eq.(2) over time and the x and z directions and integrating the averaged equation with respect to y, we obtain the relation  The profiles of these four terms are drawn in Fig. 9. They are normalized by the wall shear stress. The sum of the four terms and the straight line of (1 ) y   are also plotted in the figure. They agree well with each other, which indicates that the overall balance of forces is satisfied. The viscous shear stress is dominant in the near-wall regions as in the case of single-phase flows. The surface-tension term has large values in the regions of high void fractions (see Fig. 4). The bubbles are deformed by the mean shear in the near-wall regions and a restoring force due to the interfacial tension is acting to the liquid fluid. Since this term has large values near the walls, it makes a major contribution to the increase in the friction drag. The buoyancy term has relatively large values in the core region of the channel, which reduces the relative magnitude of the Reynolds shear stress there. For the droplet flow, the buoyancy term is obviously zero. The surface-tension term has relatively large values near the walls since some droplets are located there. Its magnitude is smaller than that in the bubbly flow, however. Instead, the relative magnitude of the Reynolds shear stress is large in the droplet flow.  Fig. 10(a). It is clearly seen that the magnitude of Reynolds shear stress is increased in the near-wall regions due to the presence of the bubbles (or droplets). This is because the momentum exchange is enhanced by the vortices generated around the bubbles (or droplets) near the walls. This increase of the Reynolds shear stress near the walls contributes to the increase in the wall friction. When normalized by W  , the amplitude of the Reynolds stress is substantially reduced in the bubbly flow except in the close vicinity of the walls. This indicates that the relative role of the turbulence in the shear stress becomes diminished in the bubbly flow.

Heat transfer characteristics of turbulent bubbly flow 3.2.1 Friction temperature and Nusselt number
The friction temperature is also altered by the injection of the bubbles or droplets. The relative magnitude of the friction temperature is shown in Table 9. The time-averaged values of the Nusselt number in the single-phase flow are 15.8 and 20.6 for Pr 1 c  and Pr 2 c  , respectively (see Table 10). In the bubbly flow, the time-averaged Nusselt number for Pr 1 c  (Case B2) and Pr 2 c  (Case B1) are 19.5 and 27.3, which are 1.23 and 1.33 times higher than the corresponding values in the single-phase flow, respectively. In the droplet flow, the Nusselt numbers for Pr 1 c  (Case D2) and Pr 2 c  (Case D1) are respectively 19.8 and 27.1, which are very close to the corresponding ones in the bubbly flow in spite of the difference in the wall shear stress. By comparing Case B1 and Case B3, it is found that the reduction in the Nusselt number due to the insulating effect of the bubbles is very small. By comparing Case B1 and Case B4, we notice that the low heat capacity of the gas inside the bubbles leads to some amount of reduction in the Nusselt number.  Table 10. Time-averaged Nusselt numbers. The values in the parentheses represent those for the bubbly flow with the lower grid resolution.

Mean temperature profiles
The profiles of the mean temperature variance, W     , are drawn in Figs. 11. The temperature variance is decreased in the whole region of the channel for the droplet flow. In the case of the bubbly flow, the temperature difference is decreased except in the core region of the channel. This increase in the core region indicates that the enhancement of fluid mixing due to the bubbles is rather confined to the near wall regions. The difference between the mean fluid temperature and the wall temperature is smaller in the multiphase flows than in the single-phase flow, which means that the increase of the Nusselt number exceeds that of the friction Reynolds number (see Eq. (17)). www.intechopen.com

Heat flux profiles
Now, we examine the mechanism for the enhancement of heat transfer by considering the energy balance in the channel. The averaging of the energy equation, Eq.(8), over time and the x and z directions, and the integration of the averaged equation with respect to y yield The left-hand side of Eq.(25) represents the total heat flux in the wall-normal direction, and the first and the second terms on the right-hand side represent the molecular heat flux and the turbulent heat flux, respectively. In the figures below, each term is normalized by the wall heat flux of each case. The profiles of these three terms are drawn in Figs. 12 for Case B1 and Case D1. The sum of the molecular and turbulent fluxes is also plotted in the figures. In both cases, the sum agrees well with the left-hand side of Eq.(25), which indicates that the overall balance of heat transfer is satisfied. Fig. 13(a) shows the profiles of the left-hand side of Eq.(25) for the single-phase flow, Case B1 and Case B2. Clear differences are not seen among three cases, indicating that the change in the profile of heat-capacity flow rate due to the bubbles or droplets has an insignificant effect on the enhancement of heat transfer. The profiles of the molecular heat flux are shown in Fig.13(b). The molecular heat flux in the bubbly (or droplet) flow is reduced near the walls compared with the case of the single-phase flow. This indicates that the effective heat conductivity is increased due to the bubbles (or droplets). Fig. 13(c) shows the profiles of turbulent heat flux. The turbulent heat flux is increased by the effects of bubbles (or droplets) in the regions near the walls. As is shown below, this is caused by the vortices whose generation is activated by the bubbles (or droplets). This increase in the turbulent heat flux in the near-wall regions is responsible for the increase in the effective heat conductivity of the fluid near the wall, which enhances the heat transfer. In summary, the enhancement of heat transfer in the bubbly or droplet flow is caused by the increase of the turbulent heat flux near the walls.

Effects of thermal properties of gas inside bubbles
Figs. 14 show the distribution of the temperature variance, , in the x-y plane for the three cases (Case B1, B3, B4) with different thermal properties for the gas inside the bubbles. Red and blue represent the regions of , respectively. Contour lines represent the cross-sections of the bubbles. The temperature is high near the walls and low in the center of the channel. The flow is going upward. It is found in Figs. 14 that the temperature field is almost uniform inside the droplet. This is due to the circulating flow inside the bubble (figures not shown). The change of the temperature distribution is small if the thermal conductivity of the gas is reduced. This is because the effect of convection dominates that of conduction inside the bubbles. In the single-phase turbulent flow, the vortices are located only in the low-temperature regions away from the walls. In the bubbly and droplet turbulent flows, on the other hand, they are also located in the near-wall regions where the temperature gradient is relatively high. It is obvious that the vortices near the walls play a major role in the heat-transfer enhancement.
It is interesting to know how the heat-transfer enhancement depends on the continuousphase Prandtl number. As is shown in Table 11, the heat-transfer enhancement is more noticeable at higher Prandtl numbers. Since the thermal boundary layer is thinner at a higher Prandtl number, the vortices near the walls, which are generated by the bubbles or droplets, more effectively enhance the heat transfer. The ratio (Pr 2) (Pr 1) Nu Nu   is weakly dependent on the friction Reynolds number in single-phase turbulent flows. We performed a simulation for the single-phase turbulence at Re 160   and obtained the value of 1.33, which is considerably lower than 1.40 in the bubbly flow at Re 160   .

 
1.31 1.40 1.37 Table 11. The ratio of the Nusselt number at Pr 2 c  to that at Pr 1 c  . The ratio of the Nusselt number in Case B1 to that in Case B2 is shown for the bubbly flow.

Performance of heat transfer enhancement
As was shown above, the Nusselt number is increased by the injection of the bubbles or droplets in the present simulations. This heat transfer enhancement is accompanied by the increase of the wall-friction, however. Reynolds analogy provides a useful concept for the evaluating the performance of heat transfer enhancement. Colburn (1933)  As is shown in Table 12, the injection of the bubbles or droplets leads to the reduction of 2 f j c . The forces resulting from the interfacial surface tension (and the buoyancy) significantly contribute to the increase of the wall shear stress in addition to the convection in turbulence. Heat transfer enhancement, on the other hand, is mainly caused by the increase in turbulent heat flux. Since the effects of the surface tension are more significant in the bubbly flow, the reduction is more noticeable in the bubbly flow than in the droplet flow.
The value of 2 f j c is larger for higher Prandtl numbers for all cases. The reduction of 2 f j c due to the injection of the bubbles or droplets is less significant for higher Prandtl numbers where the convection term plays more important roles. The above results indicate that the performance of heat transfer is not so good in the bubbly and droplet turbulent flows. In the case of bubbly flows, however, the buoyancy force exerted on bubbles, g 

Future research
We still have many issues to resolve in order to clarify the heat transfer characteristics of turbulent bubby upflows. Firstly, the size of the computational domain is possibly too small to obtain the correct statistics of the turbulent bubbly flow since we utilized a minimal channel. In addition, statistical errors may be large since the number of bubbles is small and the computational time is relatively short to obtain the turbulence statistics. Larger and longer simulations are required to resolve these problems. Secondly, we set the density ratio at 0.1 and the viscosity ratio at 1.0, which are much higher than the corresponding values in an air-water system. The bubble motions and the generation of vortices due to the bubbles should be examined by changing the values of these two ratios. The continuous-phase Prandtl number of 2.0 employed in the present study is low compared with that of about 7 at room temperature. Simulations at higher Prandtl numbers are desirable. Thirdly, only two flow patterns (the bubbly and the droplet flows) have been simulated in the present study. As was shown in Lu & Tryggvason (2008), the bubble distribution consists of a core where the flow is essentially homogeneous, and a wall layer with a larger number of bubbles sliding along the wall in the case where the buoyancy effect is dominant. This interesting situation should be examined in the future study. The heat transfer characteristics of bubbly drag-reducing flows is also an interesting topic to explore.

Conclusion
Direct numerical simulations have been conducted for turbulent bubbly upflow between two parallel heating walls at a constant volume flow rate in order to clarify its heat transfer characteristics. For comparison, simulations for neutrally buoyant droplets have also been performed. We have obtained the following results.  The bubbles accumulate in the vicinity of the wall and slide along the wall in the turbulent channel upflow.  The droplets are distributed rather uniformly throughout the channel though some tendency of accumulation in the vicinity of the walls is observed.  The turbulence production is enhanced by the bubbles or the droplets in the near-wall regions.  The wall friction is increased by the injection of bubbles. This is mainly caused by the interfacial surface tension resulting from the deformation of the bubbles due to high shear near the walls.  The heat transfer is enhanced by the injection of bubbles (or droplets). This is because the turbulent heat flux is augmented by the generation of the vortices due to the bubbles (or droplets).  The reduction in the Nusselt number due to the insulating effect of the bubbles is very small, while the low heat capacity of the gas inside the bubbles causes some amount of reduction.  The heat-transfer enhancement is more noticeable at higher Prandtl numbers.  The performance of heat transfer enhancement is not good in the bubbly and droplet turbulent flows. However, the performance is improved in the bubbly flow if the buoyancy force exerted on bubbles is available as a driving force of the upflow.

References
Aulisa The purpose of this book is to introduce researchers and graduate students to a broad range of applications of computational simulations, with a particular emphasis on those involving computational fluid dynamics (CFD) simulations. The book is divided into three parts: Part I covers some basic research topics and development in numerical algorithms for CFD simulations, including Reynolds stress transport modeling, central difference schemes for convection-diffusion equations, and flow simulations involving simple geometries such as a flat plate or a vertical channel. Part II covers a variety of important applications in which CFD simulations play a crucial role, including combustion process and automobile engine design, fluid heat exchange, airborne contaminant dispersion over buildings and atmospheric flow around a re-entry capsule, gas-solid two phase flow in long pipes, free surface flow around a ship hull, and hydrodynamic analysis of electrochemical cells. Part III covers applications of non-CFD based computational simulations, including atmospheric optical communications, climate system simulations, porous media flow, combustion, solidification, and sound field simulations for optimal acoustic effects.

How to reference
In order to correctly reference this scholarly work, feel free to copy and paste the following: Mitsuru Tanaka (2011). Numerical Study on Flow Structures and Heat Transfer Characteristics of Turbulent Bubbly Upflow in a Vertical Channel, Computational Simulations and Applications, Dr. Jianping Zhu (Ed.), ISBN: 978-953-307-430-6, InTech, Available from: http://www.intechopen.com/books/computationalsimulations-and-applications/numerical-study-on-flow-structures-and-heat-transfer-characteristics-of-turbulentbubbly-upflow-in-a