The study of multiphase flow of liquid (water) and air bubbles has been attracting the interest of many researchers. Experimental investigations of the flow including microbubbles are now well developed regarding engineering applications, such as drag reduction of pipe flow , cleaning of polymer ink  and separation of carbon fibers  using microbubbles. However, numerical investigations of the liquid flow containing air bubbles are facing difficulty in numerical accuracy and a heavy computational load. Multiphase flow has been computationally studied by the volume of fluid method (VOF)  , the discrete element method (DEM)  , the immersed boundary method (IBM) , the moving particle semi-implicit method (MPS) , and the theoretical model equations taking into account of inertial and added-mass body acceleration forces  . Although VOF can be applied to many complex configurations, it needs much CPU time and heavy memory capacity. IBM proposed by Peskin  has been largely applied to solve the liquid motion including solid particles. But it cannot be used to solve particle-liquid flow if the particle density is less than 0.3 of the liquid density . MPS also has a problem in a heavy calculation load.
Takemura et al. performed experimental studies of a rising air bubble  and a solid-like air bubble  near a vertical wall. There are some works which studied multiphase flow including plural bubbles. At low Reynolds numbers, Katz and Meneveau  conducted experiments of the interacting air bubbles rising in a quiescent distilled water, and observed the velocities of the two bubbles rising in line. Watanabe and Sanada  and Kusuno et al.  performed experiments of two rising bubbles in a quiescent silicon oil. Kusuno et al.  observed that a trailing bubble collided with a leading bubble, but two bubbles did not coalesce. Gumulya et al.  numerically studied two rising bubbles which collided and coalesced using VOF.
The force-coupling method (FCM) originally proposed by Maxey and Patel , which we call the original FCM (OFCM) in this paper, requires no moving boundaries of particles, but represents the center position of a particle by the point-body force acting from the particle to the liquid    . OFCM has been developed to calculate multiphase flow consisting of liquid and solid particles. However, the applicability of OFCM to air bubbles was notclear . Very recently, Guan et al.   proposed the modification of OFCM in order to simulate the multiphase flow composed of liquid and air bubbles, where they called it the modified FCM (MFCM), and obtained fairly good agreement between the numerical results and the experimental data by Takemura et al.  . However, inconsistency remained in the definition of the effective ranges for the force terms of MFCM in the previous papers  , where the values valid for flow containing solid particles were used for that containing air bubbles. In the present paper, we proposed the calculation method with a more appropriate definition of the effective ranges. Specifically, the effective range of the force monopole term and that of the force dipole term of OFCM were changed to be applied to air bubbles, which had not been made in the previous papers  . As a result, we reproduced more accurately the results of the experimental data of Takemura et al.  than those in the previous study  for a single air bubble rising near a vertical wall, and the calculated velocities of rising interacting two air bubbles in line were very close to those in the experimental study by Katz and Meneveau .
2. Formulation of the Modified FCM with Corrections
In OFCM and MFCM, the incompressible Navier-Stokes equations with a body-force term, , representing the interaction of the liquid with bubbles,
where indicate the x, y, z-direction components, respectively, and the equation of continuity,
are the basic equations, where the Einstein’s summation convention is used. is the velocity vector of the liquid, t the time, the density of the liquid, p the pressure, the kinematic viscosity of the liquid, and the gravitational acceleration vector. denotes the body-force acting on the liquid from N bubbles, which is given by
consists of two parts, the force monopole term (FMT, the first term on the right-hand side of Equation (3)) and the force dipole term (FDT, the second term on the right-hand side of Equation (3)). Here, is the position vector in the liquid, and is the position vector of the center of the n-th bubble. and are the effective ranges of FMT and FDT around the center of the particle, respectively, of the n-th bubble, and and are the force distribution of FMT and FDT, respectively, given by Maxey and Patel  and Lomholt et al.  as
and are related to the radius of the n-th particle as and in OFCM     . However, because these values of and were obtained by considering the theoretical solution of the Stokes approximation of the flow around a solid particle, they are not guaranteed to be applied to an air bubble.
Therefore, we proposed in this paper that and , where the former was obtained in consideration of the theoretical solution of the Stokes approximation of the flow around an air bubble, while the latter of the theoretical solution of dissipation rate for an air bubble by Batchelor  (§4.11). Please refer Appendix for the evaluation of and .These are the corrections newly introduced in the present paper to OFCM      and MFCM  . In this paper, we assumed that the radius of each bubble, , was the same, so that we put .
In Equation (3), FMT represents the reaction to the drag force acting on the bubble for its translational motion. is the force acting on the liquid from the n-th bubble, which is given by
where is the density of the bubble, the velocity of the n-th bubble. was given by the following equation in terms of the Gaussian distribution function.
in OFCM . Equation (7) means that the particle (bubble) moves with the liquid velocity averaged over the region, where the body force of FMT is effective. In the present study, we used the similar method as Guan et al. , and derived the following equation for according to Mei et al.  :
in order to apply MFCM over a wide range of the bubble Reynolds number. Here, is the bubble Reynolds number determined by the terminal velocity of a bubblerising in quiescent liquid, , given by
Note that in Guan et al. , the equation corresponding to Equation (8) (Equation (12) in Guan et al.  ) was developed for the flow including a solid particle. This is a new correction introduced in the present paper for the air bubble case.
In Equation (3), FDT represents the influence of the velocity gradient around the bubble. in Equation (3) is , where is an anti-symmetric and is symmetric parts, respectively. is related to the torque to the fluid from the bubble and is given by the equation,
where is the Eddington’s alternating sign tensor, and is calculated by use of , , and . represents a stresslet acting on the fluid and has been already given by the following equation for an air bubble .
where is the rate of strain, and determined to satisfy the constraint, for each particle using an iteration method in order to keep the bubble shape a sphere.
The calculation domain of the present study was a cube with 0.05-m-length of each side, as shown in Figure 1. The size of the computational domain was the same to the experimental ones . The gravitational acceleration, , is directed to the negative x-direction. In the numerical calculation, the Fourier transform was used in the x-direction with the collocation method, where the number of Fourier modes was . The finite-difference method was used with equally spaced staggered grids, where the number of grids in the y-direction and in the z-direction. The kinematic viscosity, , and the bubble radius, R, were taken to be equal to those of the experimental conditions.
The calculation procedure of the MFCM is summarized as follows:
1) Give the initial conditions.
2) Get from the initial conditions (Equation (6)).
3) Obtain the FMT from , and (Equations (4), and (6)).
Figure 1. Computational domain and coordinate system.
4) Get from , , , and (Equations (5), (10), (11), and (12)).
5) Obtain the FDT from , and (Equations (5) and Step 4).
6) Acquire from the FMT and the FDT (Equation (3)).
7) Solve the Navier-Stokes equations including and obtain the velocity field (Equations (1), and (2)).
8) Obtain the bubble velocity from and (Equations (8), (9) and Step 3).
9) Solve the next-time step bubble position from and .
10) Return to Step 2.
3. Comparison of the Numerical Results with the Experimental Data for a Single Air Bubble Rising near a Vertical Wall
To show better accuracy of the proposed method applied to the bubble motion in quiescent liquid than our previous method , comparison was conducted between the results of the numerical calculation using the present method and the experimental data by Takemura et al.  for the motion of a single air bubble rising near a vertical wall. Since the number of a concerning bubble was unity, was set. Figure 2 shows the test sections in the y-z plane and the position of the bubble in the experiment  and the present numerical calculation, where denotes the distance from the bubble center to the nearest vertical wall. Note that L is a function of t. The initial position of the bubble center was put at the point , where is .
Takemura et al.  presented the values related to the lift force by observing the motion of the bubble in their paper. We also calculated those values from the motion of the bubble obtained by our simulations for comparison. They introduced the wall distance Reynolds number, , defined in terms of as
Figure 2. Experimental test section and computational domain.
Note that is also a function of t. Takemura et al.  also evaluated the inverse Fourier transform of the lift force on a bubble, , using the following equation:
where is the inverse Fourier transform of the horizontal component of the drag force. We evaluated in the same manner of Takemura et al.  and compared it with the experimental ones. Note was theoretically evaluated under the assumption of the Oseen approximation  .
Figure 3(a) and Figure 3(b) show the values of as a function of for and 15, respectively. The black symbols indicate the experimental data of Takemura et al. . The black dashed lines show by Guan et al.  in which the numerical simulations of this case were conducted without the correction of the effective ranges, while the black solid lines show of the present numerical simulations using MFCM with the correction of the effective ranges. It is found that the present results are more accurate than those in the
Figure 3. The values of as a function of for and 15. (a) ; (b) .
previous study . Note that the maximum discrepancy between the experimental data  and the present numerical results for was 5%and that for was 28%. The present numerical predictions were much better than the previous ones . Thus, the present MFCM proved to be applicable to the calculation of multiphase flow consisting of water and an air bubble.
4. Comparison of the Numerical Results with the Experimental Data for Two Interacting Air Bubbles Rising in Quiescent Liquid
Then, we applied the present MFCM with the corrections of the effective ranges to the case of two interacting bubbles rising in quiescent liquid and compared the results with the experimental data of Katz and Meneveau . For this case, was set. Figure 4 shows the initial positions of two bubbles on the vertical line, , in the computational domain of the present study. The initial position of the lower bubble center was set at the point , and that of the upper bubble center at the point , where was the initial value of S, the distance between two bubbles, and was taken to be equal to the experimental value of Katz and Meneveau . Note that two bubbles moved only vertically.
Figures 5(a)-(c) show the dimensionless bubble-rising velocities, and , as a function of a dimensionless distance between two bubbles, S/R, for , 19 and 35, respectively. Solid and blank circle symbols indicate the dimensionless bubble velocities of the experimental value , and red long-dashed and dash-dotted lines are those obtained by the present calculations.
Figure 5(a) shows that the results by the present MFCM are in good agreement with the experimental values for . It is interesting that the velocities of the two bubbles became faster than the terminal velocity of the one bubble case when two bubbles were present. In the experiment  two bubbles
Figure 4. Initial position of two bubbles in the computational domain.
Figure 5. Change of the velocities of two bubbles as a function of the nondimensional distance between two bubbles for , 19, and 35. (a) ; (b) ; (c) .
finally coalesced, while in the present numerical calculations they were stuck and moved together. In Figure 6, the change of the distance between two bubbles is plotted as a function of time.
For and 35, the present numerical simulations could reproduce the experimental data  very well if two bubbles were not very close as shown in Figure 5(b) and Figure 5(c). However, when two bubbles were very close, discrepancy between the present calculation and the experimental data
Figure 6. Change of the distance between two bubbles as a function of time for .
happened for these bubble Reynolds numbers, because if two bubbles approached, deformation of the bubble surface and coalescence of bubbles might occur, which was not considered in the present calculation. Note also that the experiments of Katz and Meneveau  were on many (more than two) rising air bubbles aligned like a chain, which might cause a little difference between the present simulations and their experiments.
To understand the physical mechanism of the acceleration of the bubble velocity and the approach of two bubbles, we plot the pressure contours (the unit is Pascal [Pa]), the pressure distribution along the line and the vertical velocity distribution for in Figure 7. When two bubbles were apart from each other at t = 0.05 s as shown in Figure 7(a1) and Figure 7(a2), the pressure distribution and vertical velocity distribution were similar for two bubbles. If two bubbles came closer at t = 2.3 s as shown in Figure 7(b1) and Figure 7(b2), the high pressure region of the upper side of the lower bubble shrunk because the lower bubble entered the wake region of the upper bubble, which caused the acceleration of the lower bubble velocity. When two bubbles were very close at t = 2.9 s as shown in Figure 7(c1) and Figure 7(c2), pressure variation between two bubbles became smaller, which implies that two bubbles tended to move together. Remark that the maximum rising velocity is , which is realized if two bubbles coalesce to form a single spherical bubble whose volume is twice of each bubble. In Figure 7(a1) and Figure 7(a2) to Figure 7(c1) and Figure 7(c2), we plotted two circles on the abscissa which show the radius, R, of each bubble. It is reasonable that the pressure peak occurred at the head surface of each bubble and the valley at the tail. Therefore, we found that the present MFCM could approximately exhibit the boundary of the bubble because the effective ranges of FMT and FDT, both of which are proportional to R, were properly set in the present method.
(a1) (a2) (b1) (b2) (c1) (c2)
Figure 7. Pressure contours and velocity distribution in the plane , and pressure distribution along the line for . Circles plotted on the abscissa show the surface of each bubbles. (a1) Pressure contours and vertical velocity distribution around two bubbles in the plane at and for . (a2) Pressure distribution along the line , at and for . (b1) Pressure contours and vertical velocity distribution around two bubbles in the plane at for . (b2) Pressure distribution along the line , at and for . (c1) Pressure contours and vertical velocity distribution around two bubbles in the plane at for . (c2) Pressure distribution along the line , at and for .
5. Comparison of CPU Time and Memory Capacity between MFCM and VOF
We conducted a calculation of mutually interacting two air bubbles rising in linetreated in § 4 by using MFCM as well as VOF  to show the quickness of the present MFCM for the simulation of multiphase flow. For both VOF and MFCM calculations, we used the same computer as shown in Table 1. In the present case, we did not use the parallel computing system.
Figure 8 shows the typical results for both calculations. For this case, the bubble Reynolds number of both the bubbles was set to be 3. Figure 8 shows the distance between the two bubbles once increased and gradually decreased. Both calculations agreed with each other. We found that VOF required 28 grids (using regular grid) in one dimension within the diameter of a bubble for an accurate calculation, while MFCM required only 6 grids to capture the accurate motion of the bubble. An increase in the number of the grid for VOF required much larger computer memory and CPU time for calculations.
Required CPU time and memory for both calculations are also listed in Table 2. We found that VOF required over 50 times CPU time of MFCM and over 100 times memory of MFCM. Table 2 clearly shows the advantage of MFCM to study the multiphase flow including many air bubbles numerically with keeping the ability of the accurate prediction of bubble orbits. Remark that the study of rising two bubbles in line were dealt with numerically by VOF  , MPS  and the use of the theoretical model equations  , but we applied the present MFCM to this case, and found that the present method needed a very small calculation load, and resulted in very accurate solutions and a good physical insight into the interaction of two bubbles.
Table 1. The specifications of the computer used in this study.
Table 2. Required CPU time and memory for two bubbles calculation using MFCM and VOF.
Figure 8. Change of the distance between two bubbles as a function of time for.
In the present paper, we changed the effective ranges of the force monopole and force dipole terms in the modified force-coupling method (MFCM) proposed by Guan et al.  to calculate multiphase flow composed of liquid (water) and air bubbles more accurately. First, we performed numerical simulations of a single air bubble rising near a vertical wall. We found that the present method reproduced the experimental data more accurately than our previous method. Second, we conducted numerical simulations of mutually interacting two air bubbles rising in line. We showed that the present method accurately reproduced the experimental data over a wide range of the bubble Reynolds number. The good agreement proved that our new method is proper. Third, we compared the CPU time and memory needed for the present method with those for VOF. We found that the present method was much quicker and needed much smaller memory capacity. This indicates that our new method has the potential to provide accurate physical insights into the interaction of many bubbles with a much smaller computational load.
This study was supported by a grant-in-aid for Scientific Research for Fundamental Sciences (No. 18K12069 to K. M.) from The Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan.
In this section, we describe the derivation of and for air bubble. that can be applied to air bubble is derived from the following equations. The Stokes equation of the incompressible flow (without considering the inertia term of the Navier-Stokes equation), which acts only on the force of the FMT, is given as the following equation when considering the steady flow,
Fourier transform of the above equation gave the following equation,
Substituting the Oseen operator and assuming holds, the following equation was obtained,
Inverse Fourier transform of the above equation gives the following equation,
where, and are given as follows ,
where erf represents the error function. Considering the coordinate origin, , the following equation was obtained,
In the case of Maxey and Patel , the velocity of movement of the solid spherical particle under the action of only gravity is considered to be equivalent to the velocity of the spherical particle predicted by Stokes law, and substituting the Oseen operator into the Stokes law equation of air bubble gives,
at the coordinate origin, assuming, the above equation is as follows,
Substituting Equation (A4) into Equation (A10) gives,
Using Equation (A7) and Equation (A8), the above equation is changed to,
Solving the above equation gave the following equation,
where, assuming that, the above equation is given as follows,
From the above equation, the effective ranges of FMT for air bubble is, but considering other forces (e.g. the added mass force), we found that which was more approximate.
that can be applied to air bubble is derived from the following equations. The Stokes equation of the incompressible flow (without considering the inertia term of the Navier-Stokes equation) acting on the force of the FMT and the FDT is given as the following equation when considering the steady flow.
Fourier transform of the above equation gave the following equation,
The flow field is as follows,
where,gave the following equation,
and are given as follows ,
The local volume-averaged velocity gradient from the center of a single air bubble is given by
Fourier transform of the above equation gave the following equation,
here, without affecting the FMT by, assuming that the above relation is isotropic for, using constants is given as follows.
By the symmetry of i, k and j, m in Equation (A23), and was obtained. As a result, the local volume-averaged velocity gradient is given by,
By summing the Equations (A23), (A24) and (A25), gave the following equation,
Assuming that a single air bubble receiving torque rotates at an angular velocity, the following relationship can be obtained ,
Considering the local volume-averaged velocity gradient, the angular velocity is given as follows,
Combined with Equation (A28), using the symmetric component of the FDT gives,
Substituting Equation (A27) into Equation (A30) gives,
Then, substituting Equation (A28) into Equation (A31) gives,
Solving Equation (A32) gives,
a body-force term of FCM
velocity vector of the liquid
gravitational acceleration vector
density of the liquid
kinematic viscosity of the liquid
N total number of bubbles
n nth bubble
position vector in the liquid
position vector of the center
effective ranges of FMT
effective ranges of FDT
force distribution of FMT
force distribution of FDT
radius of the bubble
force from the bubble
density of the bubble
velocity of the bubble
the bubble Reynolds number
terminal velocity of a bubble
velocity gradient around the bubble
rate of strain
a stress let acting on the fluid
torque to a fluid from nth bubble
bubble-nearest vertical wall distance
initial value of
S distance between two bubbles
initial value of S
wall distance Reynolds number
 Shatat, M.M.E., Yanase, S., Takami, T. and Hyakutake, T. (2009) Drag Reduction Effects of Micro-Bubbles in Straight and Helical Pipes. Journal of Fluid Science and Technology, 4, 156-167.
 Matsuura, K., Ogawa, S., Kasaki, S., Koyama, K., Kodama, M. and Yanase S. (2015) Cleaning Polymer Ink from a Glass Substrate Using Microbubbles Generated by Hydrogen Bubble Method. Separation and Purification Technology, 142, 242-250.
 Matsuura, K., Uchida, T., Guan, C. and Yanase, S. (2018) Separation of Carbon Fibers in Water Using Microbubbles Generated by Hydrogen Bubble Method. Separation and Purification Technology, 190, 190-194.
 Gumulya, M., Utikar, R.P., Evans, G.M., Joshi, J.B. and Pareek, V. (2017) Interaction of Bubbles Rising Inline in Quiescent Liquid. Chemical Engineering Science, 166, 1-10.
 Fuchimoto, T., Yanase, S., Mizushima, J. and Senda, J. (2009) Dynamics of Vortex Rings in the Spray from a Swirl Injector. Fluid Dynamics Research, 41, Article ID: 045503.
 Peskin, C.S. (1972) Flow Patterns around Heart Valves: A Digital Computer Method for Solving the Equations of Motion. PhD Thesis, Albert Einstein College of Medicine, University Microfilms, 72-30.
 Chen, R.H., Tian, W.X., Su, G.H., Qiu, S.Z., Ishiwatari, Y. and Oka, Y. (2011) Numerical Investigation on Coalescence of Bubble Pairs Rising in a Stagnant Liquid. Chemical Engineering Science, 66, 5055-5063.
 Ramírez-Munoz, J., Baz-Rodríguez, S., Salinas-Rodríguez, E., Castellanos-Sahagún, E. and Puebla, H. (2013) Forces on Aligned Rising Spherical Bubbles at Low-to-Moderate Reynolds Number. Physics of Fluids, 25, Article ID: 093303.
 Kempe, T. and Frohlich, J. (2012) An Improved Immersed Boundary Method with Direct Forcing for the Simulation of Particle Laden Flows. Journal of Computational Physics, 231, 3663-3684.
 Takemura, F., Takagi, S., Magnaudet, J. and Matsumoto, Y. (2002) Drag and Lift Forces on a Bubble Rising Near a Vertical Wall in a Viscous Liquid. Journal of Fluid Mechanics, 461, 277-300.
 Takemura, F., Takagi, S. and Matsumoto, Y. (2002b) Lift Force on a Bubble Rising Near a Wall in Water below Reynolds Number of 20. Transactions of the Japan Society of Mechanical Engineers, Series B, 68, 1684-1690 (In Japanese).
 Maxey, M.R. and Patel, B.K. (2001) Localized Force Representations for Particles Sedimenting in Stokes Flow. International Journal of Multiphase Flow, 27, 1603-1626.
 Xu, J., Maxey, M.R. and Karniadakis, G.E. (2002) Numerical Simulation of Turbulent Drag Reduction Using Micro-Bubbles. Journal of Fluid Mechanics, 468, 271-281.
 Lomholt, S., Stenum, B. and Maxey, M.R. (2002) Experimental Verification of the Force Coupling Method for Particulate Flows. International Journal of Multiphase Flow, 28, 225-246.
 Guan, C., Yanase, S., Matsuura, K., Kouchi, T. and Nagata, Y. (2017) Application of the Modified Force-Coupling Method of Tracing the Trajectories of Spherical Bubbles with Solid-Like and Slip Surfaces. Open Journal of Fluid Dynamics, 7, 657-672.
 Guan, C., Yanase, S., Matsuura, K., Kouchi, T. and Nagata, Y. (2018) Numerical Study on Motion of a Single Spherical Bubble Near Vertical Wall in the Quiescent Fluid Using the Modified Force-Coupling Method. Nagare, 37, 281-289 (In Japanese).
 Vasseur, P. and Cox, R.G. (1977) The Lateral Migration of Spherical Particle Sedimenting in Stagnant Bounded Fluid. Journal of Fluid Mechanics, 80, 561-591.
 OpenFOAM. The Open Source CFD Toolbox.