Numerical Study of the Lift Force, Velocities and Pressure Distribution of a Single Air Bubble and Two Interacting Air Bubbles Rising in Quiescent Liquid

The study of multiphase flow consisting of liquid and air bubbles has been attracting the interest of many researchers. Numerical methods for such a system are, however, facing difficulty in numerical accuracy and a heavy computational load. In this paper, we made corrections to the modified force-coupling method in our previous papers and applied it to the numerical studies of a single air bubble rising near a vertical wall and two interacting air bubbles rising in line in quiescent liquid. Corrections were made to the effective ranges of the force-coupling method. The calculation results showed that the lift force acting on an air bubble obtained by the experimental data was more accurately reproduced than those by our previous method. We accurately calculated the time evolution of the velocities of interacting two air bubbles rising in line obtained in the previous experiments and resolved the physical mechanism of the relative movement of two bubbles. We also found the present method is much quicker and needs much smaller memory capacity than other methods, such as the volume of fluid method.


Introduction
The study of multiphase flow of liquid (water) and air bubbles has been attracting the interest of many researchers. Experimental investigations of the flow in-cluding microbubbles are now well developed regarding engineering applications, such as drag reduction of pipe flow [1], cleaning of polymer ink [2] and separation of carbon fibers [3] 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) [4] [5], the discrete element method (DEM) [6] [7], the immersed boundary method (IBM) [8], the moving particle semi-implicit method (MPS) [9], and the theoretical model equations taking into account of inertial and added-mass body acceleration forces [10] [11].
Although VOF can be applied to many complex configurations, it needs much CPU time and heavy memory capacity. IBM proposed by Peskin [8] 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 [12]. MPS also has a problem in a heavy calculation load. Takemura et al. performed experimental studies of a rising air bubble [13] and a solid-like air bubble [14] near a vertical wall. There are some works which studied multiphase flow including plural bubbles. At low Reynolds numbers, Katz and Meneveau [15] 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 [16] and Kusuno et al. [17] performed experiments of two rising bubbles in a quiescent silicon oil. Kusuno et al. [17] observed that a trailing bubble collided with a leading bubble, but two bubbles did not coalesce. Gumulya et al. [5] numerically studied two rising bubbles which collided and coalesced using VOF.
The force-coupling method (FCM) originally proposed by Maxey and Patel [18], 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 [19] [20] [21] [22]. OFCM has been developed to calculate multiphase flow consisting of liquid and solid particles. However, the applicability of OFCM to air bubbles was notclear [19].
Very recently, Guan et al. [23] [24] 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. [13] [14]. However, inconsistency remained in the definition of the effective ranges for the force terms of MFCM in the previous papers [23] [24], 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 [23] [24]. As a result, we reproduced more accurately the results of the experimental data of Takemura et al. [13] than those in the previous study [23]

Formulation of the Modified FCM with Corrections
In OFCM and MFCM, the incompressible Navier-Stokes equations with a body-force term, , b i f , representing the interaction of the liquid with bubbles, 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. i u is the velocity vector of the liquid, t the time, f ρ the density of the liquid, p the pressure, ν the kinematic viscosity of the liquid, and i g the gravitational acceleration vector.
, b i f denotes the body-force acting on the liquid from N bubbles, which is given by , b i f 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) [20] as  [24]. In this paper, we assumed that the radius of each bubble, ( ) n R , 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 ρ is the density of the bubble, ( ) n i U the velocity of the n-th bubble.
in OFCM [15]. 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. [23], and derived the following equation for ( ) n i U according to Mei et al. [26]: (8) in order to apply MFCM over a wide range of the bubble Reynolds number. Here, ( ) n Re ∞ is the bubble Reynolds number determined by the terminal velocity of a bubblerising in quiescent liquid, ( ) n U ∞ , given by Note that in Guan et al. [23], the equation corresponding to Equation (8) (Equation (12) in Guan et al. [23]) 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 is an an- is related to the torque to the fluid from the bubble and is given by the equation, where ijk  is the Eddington's alternating sign tensor, and ( )( ) represents a stresslet acting on the fluid and has been already given by the following equation for an air bubble [23].
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 [13]. The gravitational acceleration, , (11), and (12)).

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 [23], comparison was conducted between the results of the numerical calculation using the present method and the experimental data by Takemura et al. [13] for the motion of a single air bubble rising near a vertical wall. Since the number of a concerning bubble was unity, 1 N = was set. Figure 2 shows the test sections in the y-z plane and the position of the bubble in the experiment [13] and the present numerical calculation, where ( ) L t 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 ( ) Takemura et al. [13] 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,  Re is also a function of t. Takemura et al. [13] also evaluated the inverse Fourier transform of the lift force on a bubble, L I , using the following equation: where D I is the inverse Fourier transform of the horizontal component of the drag force. We evaluated L I in the same manner of Takemura et al. [13] and compared it with the experimental ones. Note D I was theoretically evaluated under the assumption of the Oseen approximation [13] [27]. previous study [23]. Note that the maximum discrepancy between the experimental data [13] and the present numerical results for ( ) 1

2.1
Re ∞ = was 5%and that for ( ) 1 15 Re ∞ = was 28%. The present numerical predictions were much better than the previous ones [23]. Thus, the present MFCM proved to be applicable to the calculation of multiphase flow consisting of water and an air bubble.

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 [15]. For this case, 2 N = was set. Figure 4 shows the initial positions of two bubbles on the ver-    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. produce the experimental data [15] 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 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 [15] 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 0 y z = = and the vertical velocity distribution for  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

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 [28] 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 [4] [5], MPS [9] and the use of the theoretical model equations [10] [11], 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.

Conclusion
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. [23] 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.

Appendix
In this section, we describe the derivation of ( ) Inverse Fourier transform of the above equation gives the following equation,