12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19-22, 2022

A DIRECT NUMERICAL SIMULATION STUDY OF TWO-PHASE TURBULENT FLOWWITH LARGE-SIZED MULTIPLE BUBBLES IN A HORIZONTAL FLAT CHANNEL

Sang Won Kim
Complex Phenomena Unified Simulation Research
Team, RIKEN Center for Computational Science
7 Chome-1-26 Minatojima Minamimachi,Chuo Ward, Kobe, Hyogo, Japan
sangwon.kim@riken.jp

Nobuyuki Oshima
Division of Mechanical and Aerospace Engineering
Hokkaido University
Nishi 8 Chome Kita 13 Jo, Kita Ward,Sapporo, Hokkaido, Japanoshima@eng.hokudai.ac.jp

Yuichi Murai
Division of Mechanical and Aerospace Engineering
Hokkaido University
Nishi 8 Chome Kita 13 Jo, Kita Ward,Sapporo, Hokkaido, Japanmurai@eng.hokudai.ac.jp

Hyun Jin Park
Division of Mechanical and Aerospace Engineering
Hokkaido University
Nishi 8 Chome Kita 13 Jo, Kita Ward,Sapporo, Hokkaido, Japan
park@eng.hokudai.ac.jp

ABSTRACTThe horizontal turbulent channel flows with injected
multiple large-sized bubbles are performed through direct
numerical simulation. We study the drag reduction mechanismof large-sized bubbles by analyzing four different volumes ofair, while the frictional Reynolds number is kept constant. We
apply a VOF model to compute the two immiscible fluids
water and air, the sharp interface between fluids captured by
the isoAdvector method. We established our simulation setup
using experimental data to compare with it in further studies.
Simulation results show that the drag reduction ratio isindependent of the void fractions ratio. From the instantaneous
bubble shapes and skin friction contours, the higher dragreduction occurs at the liquid film regime of sufficiently
longer bubbles in the streamwise direction. Thus, an additionalparameter such as streamwise bubble length is shown to
influence the drag reduction mechanism.

reduction from the downstream owing to that medium size ofbubbles (2–3 mm) contained in microbubbles frequently
coalesce in the shear layer and become large-sized bubbles (>5 mm) in the downstream region. Therefore, the bubble size atdownstream of microbubble methods is important indetermining the total drag reduction, when considering thepersistency of drag reduction along the streamwise direction.
However, the critical bubble size has not been determined.
In this situation, the drag reduction performance of large-sizedbubbles was found to be dependent on the bubble length(Murai et al., 2007). Oishi and Murai (2014) rigorously
investigated and described influential characteristics such as
the velocity gradient and u'v' contours related to the dragreduction of a single large-sized bubble. Park et al. (2019)
experimentally investigated the film thickness of large-sizedbubbles and BDR. Kim et al. (2020) numerically realized asingle large-sized bubble and investigated the flow fieldaround a large-sized bubble. Recently, Kim et al. (2021, 2022)
reported bubble shape and skin friction characteristics asseveral series of large-sized bubble volumes.
In this study, a direct numerical simulation of large-sizedbubbly flow in turbulent channel flow was performed toinvestigate the drag reduction mechanism of large-sizedbubbles. we present four different cases as volumes of air
increase whereas the frictional Reynolds number is kept
constant for investigating drag reduction performance oflarge-sized bubbles independent with the effect of turbulent
shear flows.

INTRODUCTION
Bubble drag reduction has been researched for over 40
years and applied to liquid transport in pipelines and shipsurfaces in water. This method does not damage theenvironment compared to other techniques and does not
require appendages or modifications on the wall surface. A
mostly investigated situation of drag reduction by the air
bubble injection involves microbubbles smaller than 1 mm.
Since the first report by McCormick and Bhattacharyya (1973),
many researchers have conducted experiments on microbubble drag reduction to study the parameter dependency andskin friction reduction process (Madavan et al. (1985); Kato etal. (1999)).
Meanwhile, Moriguchi and Kato (2002), Shen et al. (2006)
reported that the bubble size dependency on micro-bubbles
drag reduction is no significant effect in the case of the highshear rate. This implies that bubble size is not required to bevery small. Recently, the discussion on the drag reduction
mechanism of the bubble was shifted to its deformability, by
Moriguchi and Kato (2002), Lu et al. (2005). Althoughshowing significant drag reduction from the experimental
scale, this method has a problem with the persistence of drag

NUMERICAL METHOD
In this study, the flow is considered incompressible, andthe fluids are considered Newtonian. The governing equationsin this study for continuity and momentum equation can be
expressed by Eq. (1) and (2).

1u = 0(1) u+   (  uu) = −  p +   T +  f + ft(2)

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19-22, 2022
In a two-phase flow, the integral form of the continuity
equation can be expressed in Eq. (3) by change density term,
ρ(x, t), because the density field has a discontinuity at theinterface. The indicator field in Eq. (4) is defined to remove
these insignificant parameters, and the integral form of thecontinuity equation is reformed as Eq. (5).

d
 (x, t )dV = −   (x, t )u(x, t )  dSv
dt v
 (x, t ) −  B
H (x, t )  A − BdH (x, t )dV = −  Sij  H (x, t )u(x, t )  dSFjdt cijBilines in the surface cell, then split the cell and redistributevolume fraction so that sub-cells occupy α = 1 or 0. (5) Use
the velocity field data to estimate interface motion during atime interval. By this procedure, suppressing numerical
diffusion and precise advection can be achieved.
The numerical domain was modeled according to theexperimental channel set-up of Oishi and Murai (2014). The
length and width of the computational domain can beshortened by the periodic conditions in the streamwise andspanwise directions to save computational resources. Basedon the experimental conditions, the domain dimensions wereas follows: height (H = 2h) = 10 mm, where h is the half-heightof the channel, length (L) = 100 mm, and width (W) = 70 mm.
Table 1 lists the experimental and numerical conditions of thepresent study, and a geometrical description of the numerical
domain is shown in Fig. 1. The present study generated aturbulent Couette flow in the channel to compare theexperimental results and realize similar flow conditions withship surface.

(3)
(4)
(5)

The left-hand side of Eq. (5) can be expressed by thevolume fraction of fluid A in cell i as Eq. (6). On the otherhand, the Right-hand side of Eq. (5) can be updated to total
fraction during the time step, ΔVj (t, Δt), as Eq. (7). From theseequations, Eq. (5) can be expressed as Eq. (8).

 i (t ) =1
H (x, t )dV
Vi iV j (t , t )  t +tt

 j (t )SjFi1
 Sij V j (t, t )Vi jBi(7)100 × 10 × 70

998.7

998.7

1.0(8)

1.11Uwall [m/s]

−

2.16

Fr

3.19

3.54

ReReτ

9260

10278
260(9)dSAj ( )   H (x, )dS = F j A ( )Fj j (t ) t +tSj

6000 × 10 × 100

Umean [m/s]

The remaining surface integral is simply the instantaneous
area of face j submerged in fluid A, which is called Aj (τ) in Eq.
(10). By using the above definition, Eq. (7) is now written as
Eq. (11).

V j (t , t ) Present study

L×H×W[mm]ρ [kg/m3]

To simplify ΔVj (t, Δt) the term, some assumptions areneeded in the isoAdvector method. The first assumption is u(x,τ) ≈ u(x, t), which is regarding the velocity field as constant τduring the interval of time step [t, t+∆t]. Another assumption
is that u on the face Fj dotted with the differential face normal
vector, dS, can be approximated in terms of the constantvolumetric face flux, ϕj, as follows:u (x, t )  dS Oishi and Murai (2014)(6)

 H (x, )u(x, )  dSd (t + t ) = i (t ) −Table 1. Experimental & numerical conditions.

tA j ( )ddSFigure 1. Schematic of the numerical conditions.

(10)

The velocity of the upper wall in the channel (Uwall) for
the Couette flow was adjusted to obtain the same friction
velocity (uτ) of the experimental study to compare the dragreduction ratio. All dimensionless numbers were obtained
from single-phase flow conditions. Numerical simulation wasperformed using the open-source CFD package OpenFOAM®
v2006. The PimpleFoam solver was used to generate a fully
developed turbulent flow, whereas the InterIsoFoam solver
implemented the IsoAdvector method for the interface
sharpening method was used for two-phase flow. Furtherdetails, such as the verification work on InterIsoFoam, arepresented in Kim et al. (2020).
All schemes in this study provide second-order accuracy.
The global time step was set to 2.0 × 10−6 s to ensure that theCourant–Friedrichs–Lewy (CFL) number remained below 0.5;
the interface CFL number was less than 0.3, in each simulation(11)

The isoAdvector algorithm applied introduces a sub-gridscale model for estimating Aj (τ) and intracellular distribution
from the given volume fraction data to capture this local
distribution information. The procedure of constructing theinterface can be summarized as follows. (1) Find surface cells
that contain volume fraction, 0 < α < 1. (2) Interpolate volume
fraction of a vertex of surface cells from all surrounding cells.
(3) Calculate face-interface intersection points between twovertexes. (4) Construct interface from intersection points and2

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19-22, 2022case, because the interface CFL number frequently shot up
above the unit owing to the breakup and coalescence of thebubble. Further, both solvers use the Pimple algorithm forpressure-velocity coupling in a segregated manner, which
computes a velocity field to satisfy the continuity andmomentum equations using an iterative procedure.

Figure 2 shows the mean velocity of the bubbly flow as
the volume fraction of air. Mean velocity profiles were almost
constant with showing small fluctuations and decreasing near
the top wall owing to the existence of large-sized bubbles.
With these trends, the velocities near the bottom wall slightly
increased.
Figure 3 shows the relationship between the volume
fractions and projected void fractions. The projected void
fractions uniformly increased as volume fractions, and they
were approximately twice the volume fractions. Figure 4
shows the relationship between the volume fractions and mean
skin friction coefficients, where all plots indicate timeaveraged values in 25ms. Consequently, similar trends areshown with each period, T1-T4, and the values are gradually
decreased and converged skin friction as time passed.
Typically, the drag reduction ratio is proportional as theprojected void fraction increases. However, present results
have not followed these trends. Thus, we investigated theinstantaneous snapshot of bubble shapes and the contours ofskin friction coefficient. Figures 5 and 6 show four typicalsamples of them, where the case numbers are assigned inchronological order.
As can be observed, the bubbles are coalesced and becomelarger as the volume fraction increases, and finally the flow
state becomes to resemble a churn flow. From the contour ofthe skin friction coefficient, the liquid films of some largesized bubbles, those are relatively small or aligned in a
spanwise direction, show high skin friction similar to theregion without the bubbles. Meanwhile, the liquid film of
relatively large bubbles provides lower skin friction than thatin other regions and the turbulent shear flow is disappeared in
the regions. We expected that turbulent coherent structures are
too large to intrude flows in liquid films and cause therelaminarization on that region, and such large bubbles areshown in C24 cases.
Moreover, by measuring the streamwise bubble length as
spacing 1mm in spanwise direction, C24 clearly shows longer
length distribution over 40mm while other cases are shown
similar trends with each other, as shown in Figure 7. Thesetrends are corresponding with those of skin friction
coefficients.

NUMERICAL RESULTS
The numerical simulations are comprised of four cases as
the volume fraction of air, αvol, as described in Table 2. All
cases injected bubbles after the boundary layer is fully
developed and additional transient calculations wereperformed to achieve a steady state after the bubble injectionin the channel.
Table 2. Simulation cases for large-sized bubbly flows.
Flow condition
Turbulent Couetteflow

Case

αvol.

C9C140.0957 (9.57 %)0.1436 (14.36 %)C19C240.1908 (19.08 %)0.2469 (24.69 %)Figure 2. Mean velocity profile of bubble flow as volume
fractions.

Figure 4. Mean skin friction coefficients as averaged period.
Figure 3. Relationship between volume fraction of air andProjected void fractions.

312th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19-22, 2022
CONCLUSION
In this study, a direct numerical simulation of the largesized bubbly flow in turbulent channel flows was performed.
The numerical results show reasonable even conducting
simulations for longer periods. Interestingly, near-zero skinfriction contours on the liquid film of the large-sized bubble
are successfully confirmed from sufficiently large-sizedbubbles, especially the C24 case, and this effect is maximized
when the bubble reaches sufficiently longer in a streamwise
direction.

ACKNOWLEDGEMENTS
This Research and development work was supported byhigh-performance computing projects (hp 200223, 210160).

REFERENCES
McCormick, M., and Bhattacharyya, R., 1973. Drag
reduction of a submersible hull by electrolysis. Nav. Eng. J.
85(2), 11–16.
Madavan, N.K., Deutsch, S., and Merkle, C.L., 1984.
Reduction of turbulent skin friction by microbubbles. Phys.
Fluids 27(2), 356-363.
Kato, H., Iwashina, T., Miyanaga, M., and Yamaguchi, H.,
1999. Effect of microbubbles on the structure of turbulence ina turbulent boundary layer. J. Marine Sci. Tech. 4, 115-162.
Moriguchi, Y., and Kato, H., 2002. Influence ofmicrobubble diameter and distribution on frictional resistance
reduction. J. Mar. Sci. Technol. 7(2), 79–85.
Shen, X., Perlin, M., and Ceccio S.L., 2006. Influence ofbubble size on microbubble drag reduction. Exp. Fluids 41,415-24.
Lu, J., Fernandez, A., and Tryggvason, G., 2005. The
effect of bubbles on the wall drag in a turbulent channel flow.
Phys. Fluids 17(9), 095102.
Murai, Y., Fukuda, H., Oishi, Y., Kodama, Y., andYamamoto, F., 2007. Skin friction reduction by large airbubbles in a horizontal channel flow. Int. J. Multiphase Flow
33(2), 147–163.
Oishi. Y., and Murai. Y., 2014. Horizontal turbulent
channel flow interacted by a single large bubble. Exp. Therm.
Fluid Sci. 55, 128–139.
Park, H.J., Tasaka, Y., and Murai, Y., 2019. Bubbly dragreduction investigated by time-resolved ultrasonic pulse
echography for liquid films creeping inside a turbulent
boundary layer. Exp. Therm. Fluid Sci. 103, 66–77.
Kim, S.W., Oshima, N., Murai, Y., and Park, H.J., 2020.
Numerical investigation of a single intermediate-sized bubble
in horizontal turbulent channel flow. J. Fluid Sci. Technol.
15(3).
Kim, S.W., Oshima, N., Murai, Y., and Park, H.J., 2021.
Direct numerical simulation of frictional drag modulation inhorizontal channel flow subjected to single large-sized bubble
injection. Int. J. Multiphase Flow, 145, 103838.
Kim, S.W., Oshima, N., and Park, H.J., 2022. Numericalflow visualization of a single large-sized bubble in turbulent
Couette flow using OpenFOAM. J. Vis., (Accepted).

Figure 5. Instantaneous snapshot on bubble shapes.

412th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19-22, 2022
40C9Percent [%]302010040 01020304050607080

90

bubble length, dst [mm]C142080

Percent [%]302010040 010304050607090

bubble length, dst [mm]

C192080

Percent [%]302010040 010304050607090

bubble length, dst [mm]C242080

Percent [%]3020100010304050607090

bubble length, dst [mm]Figure 6. Instantaneous snapshot on contours of skin friction
coefficients.

Figure 7. Instantaneous distributions of streamwise bubble
length.

5