1. Introduction

The study of the multiphase flow of water including air bubbles has been attracting interest of many researchers, not only because of a lot of mechanical and chemical engineering applications, but also out of physical interest. When the diameter of an air bubble is around several tens of micrometers, it is called a microbubble, and a nanobubble if the diameter is around hundreds of nanometers. Experimental studies of microbubbles are now well developed in engineering applications, such as drag reduction of pipe flow or boundary layer flow [1] and cleaning of hydrophobic substances [2] . Regarding nanobubble flow, however, the experimental study remains in the early stages, because optical detection of sub-microscale particles is difficult and the method of nanobubble generation has not been well established.

Theoretical and numerical studies of microbubble and nanobubble flows are on the other hand very retarded because numerical calculation of such flows requires a significant computational load. An appropriate numerical method has not been found in spite of many methods proposed so far. Such methods include the volume of fluid method (VOF) [3] , the discrete element method (DEM) [4] , and the immersed boundary method (IBM) [5] . VOF is based on the fractional volume of the fluid, and has been shown to be more flexible and efficient than the other methods for treating complicated free boundary configurations. Consequently, the calculation cost becomes very large for reducing the numerical error due to an increased number of volume elements. DEM tracks individual particle motion at every instant based on the momentum equation [6] . Although the interaction between a particle and fluid is included in DEM, this method cannot be applied to the liquid fluid including air bubbles because the mass of a bubble is too small. IBM was proposed by Peskin [5] and has been developed by many authors to solve the fluid motion including solid particles. In IBM, the boundary condition on the particle surface is satisfied by applying an appropriate volume force. However, IBM cannot be used to solve particle-liquid flow in which particle density is less than 0.3 times that of the liquid [7] .

The force-coupling method (FCM) originally proposed by Maxey and Patel [8] does not require moving boundaries to satisfy the precise boundary conditions on the particle surfaces involved. Instead of that it introduces the body force from the particle to the surrounding fluid in the Navier-Stokes equations. It is noteworthy that FCM can calculate fluid motion including particles with a relatively small computational load in the case of water including air bubbles [9] . However, the accuracy of the FCM concerning the trajectories of bubbles has not been investigated yet by comparison with experimental results.

Note that there exist a few experimental studies on the bubble motion in the quiescent fluid. Takemura et al. [10] [11] published their experimental results on the motion of bubbles near a wall. They used an experimental device capable of simultaneously measuring the velocity and diameter of a rising bubble in fluid, while a microscope camera moved upward with this bubble near the vertical wall. They measured the motion of a spherical air bubble with a slip boundary in silicone oil, as well as the lift force acting on the bubble in the direction away from the wall. In order to examine the terminal velocity (or drag) and the lift force, they proposed theoretical formulas by modifying the analysis of Vasseur and Cox [12] , and compared them with their experimental results. Then, in the [13] , they measured the lift force acting on the bubble in the direction away from the wall when the bubble surface was like a solid sphere by adding surfactant into water to induce the Marangoni effect [14] . A theoretical formula for the lift force was also proposed in a similar manner to the slip boundary case and compared with their experimental results. Although they found good agreement between the theoretical formulas and the experimental results, the validity was limited in the range of small Reynolds numbers (<1) and the formulas were only applicable to the asymptotic values of the terminal velocity and lift force. To the best of the authors’ knowledge, there is no practical simulation method for precisely calculating the trajectory of a rising bubble near a wall or the fluid motion around such a bubble.

In the present paper, we proposed corrections to the original FCM in order to trace the trajectories of spherical bubbles with no-slip and slip surfaces. The bubble with a no-slip surface or a slip surface will be called the no-slip bubble or the slip bubble, respectively, henceforth in this paper. Then, we compared the results of the numerical simulation with the experimental results of Takemura et al. [10] [11] [13] . We used the same kinematic viscosity of the fluid and bubble radius over a wide range of the Reynolds number (<20) to confirm the validity of our new theory. This study gives a new method to elucidate the mechanism of the lift force generation from the vertical wall, and accurately predicting the motion of a bubble. This work is the first step in the simulation method to calculate fluid motion including many bubbles.

2. Original Force-Coupling Method

In the original version of the FCM, the incompressible Navier-Stokes equations with a force term representing the interaction with bubble,

$\frac{\partial u_i}{\partial t}=-u_j\frac{\partial u_i}{\partial x_j}-\frac{1}{{\rho }_f}\frac{\partial p}{\partial x_i}+\nu \frac{{\partial }^2{u}_i}{\partial x_j{}^2}+g_i+f_{b,i},\left(i=1,2,3\right)$ (1)

and the equation of continuity,

$\frac{\partial u_i}{\partial x_i}=0,$ (2)

are used as the basic equations, where the Einstein’s summation convention is used. $u_i$ is the velocity vector of the fluid, t the time, ${\rho }_f$ the density of the fluid, p the pressure, $\nu$ the kinematic viscosity, $g_i$ the gravitational acceleration vector. $f_{b,i}$ represents the body force acting on the fluid from N bubbles given by

$f_{b,i}=\underset{n=1}{\overset{N}{{\displaystyle \sum }}}\left\{F_i^{\left(n\right)}{\Delta }_M\left(x_i-Y_i^{\left(n\right)},{\sigma }_M\right)+G_{ij}^{\left(n\right)}\frac{\partial }{\partial x_j}{\Delta }_D\left(x_i-Y_i^{\left(n\right)},{\sigma }_D\right)\right\},$ (3)

where $i=1,2,3$ indicate the $x,y,z$ -direction components, respectively. $f_{b,i}$ 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, $x_i$ is the spatial position vector in the fluid, and $Y_i^{\left(n\right)}$ is the position vector of the center of the n-th bubble. ${\Delta }_M$ and ${\Delta }_D$ are the force ranges of FMT and FDT, respectively, in terms of the Gaussian distribution function given by

${\Delta }_M\left(x_i,{\sigma }_M\right)=\frac{1}{{\left[2\text{π}{\left({\sigma }_M\right)}^2\right]}^{3/2}}\exp \left[-\frac{x_i^2}{2{\left({\sigma }_M\right)}^2}\right],$ (4)

${\Delta }_D\left(x_i,{\sigma }_D\right)=\frac{1}{{\left[2\text{π}{\left({\sigma }_D\right)}^2\right]}^{3/2}}\exp \left[-\frac{x_i^2}{2{\left({\sigma }_D\right)}^2}\right],$ (5)

according to [8] [15] [16] . The force ranges, ${\Delta }_M$ and ${\Delta }_D$ , have a close relationship with the bubble radius, to which they are proportional. In the present study, the bubble radius is assumed to be the same for all bubbles. Here, ${\sigma }_M$ and ${\sigma }_D$ are the widths of the Gaussian distribution of FMT and FDT, respectively, which are related to the bubble radius R by ${\sigma }_M=R/\sqrt{\text{π}}$ and ${\sigma }_D=R/{\left(6\sqrt{\text{π}}\right)}^{1/3}$ . These values of ${\sigma }_M$ and ${\sigma }_D$ were set by Maxey et al. [8] [15] [16] considering the theoretical solution of the Stokes approximation of a solid particle. Therefore, they cannot be applied to a slip bubble, however, their values ${\sigma }_M$ and ${\sigma }_D$ were used in this paper because the experimental results were successfully reproduced, as shown later.

FMT represents the reaction to the drag force acting on the bubble due to the influence of its translational motion. $F_i^{\left(n\right)}$ is the force acting on the fluid from the n-th bubble, which in the original FCM [9] was given by

$F_i^{\left(n\right)}=\frac{4}{3}\text{π}{R}^3\left({\rho }_f-{\rho }_b\right)\left(g_i-\frac{\text{d}{U}_i^{\left(n\right)}}{\text{d}t}\right),$ (6)

where ${\rho }_b$ is the density of the bubble. The velocity of the n-th bubble, $U_i^{\left(n\right)}$ , is given by the following equation in terms of the Gaussian distribution function.

$U_i^{\left(n\right)}={\displaystyle \iiint u_i\left(x_i,t\right){\Delta }_M\left(x_i-Y_i^{\left(n\right)},{\sigma }_M\right)\text{d}{x}_1\text{d}{x}_2\text{d}{x}_3},$ (7)

which means that the bubble moves with the fluid velocity averaged nearly over the bubble region.

FDT represents the influence of the velocity gradient around the bubble, where $G_{ij}^{\left(n\right)}$ is the sum of two parts, an antisymmetric part, $T_{ij}^{\left(n\right)}$ , and a symmetric part, $S_{ij}^{\left(n\right)}$ . $T_{ij}^{\left(n\right)}$ represents the torque on the fluid from the bubbles and is given by the equation,

$T_{ij}^{\left(n\right)}=\frac{1}{2}{ϵ}_{ijk}{T}_k^{\left(n\right)\left(ext\right)}$ (8)

where ${ϵ}_{ijk}$ is the Eddington’s alternating sign tensor, and $T_k^{\left(n\right)\left(ext\right)}$ is an external torque on the n-th bubble. The symmetric part, $S_{ij}^{\left(n\right)}$ , represents a stresslet acting on the fluid and is given by the following equations for a no-slip bubble,

$S_{ij}^{\left(n\right)}=\frac{20}{3}\text{π}\nu {\rho }_f{R}^3{E}_{ij}^{\left(n\right)},$ (9a)

and for a slip bubble,

$S_{ij}^{\left(n\right)}=\frac{8}{3}\text{π}\nu {\rho }_f{R}^3{E}_{ij}^{\left(n\right)},$ (9b)

respectively. $E_{ij}^{\left(n\right)}$ is the rate of strain, which determined to satisfy the constraint,

$E_{ij}^{\left(n\right)}=\frac{1}{2}{{\displaystyle \int }}^{\text{}}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right){\Delta }_D\left(x_i-Y_i^{\left(n\right)},{\sigma }_D\right)\text{d}{x}_1\text{d}{x}_2\text{d}{x}_3=0,$ (10)

for each bubble using an iteration method [15] [16] [17] . Equations (1)-(10) give the complete set of the governing equations for the system of the fluid and bubbles [8] [9] .

The calculation procedure of the original FCM is summarized as follows:

1) Give the initial conditions.

2) Get $F_i^{\left(n\right)}$ from the initial conditions (Equation (6)).

3) Obtain the FMT from ${\Delta }_M\left(x_i,{\sigma }_M\right)$ , and $F_i^{\left(n\right)}$ (Equations (4), and (6)).

4) Get $G_{ij}^{\left(n\right)}$ from ${\Delta }_D\left(x_i,{\sigma }_D\right)$ , $T_{ij}^{\left(n\right)}$ , $S_{ij}^{\left(n\right)}$ , and $E_{ij}^{\left(n\right)}$ (Equations (5), (8), (9), and (10)).

5) Obtain the FDT from ${\Delta }_D\left(x_i,{\sigma }_D\right)$ , and $G_{ij}^{\left(n\right)}$ (Equations (5), (8), and (9)).

6) Acquire $f_{b,i}$ from the FMT and the FDT (Equation (3)).

7) Solve the Navier-Stokes equations including $f_{b,i}$ and obtain the velocity field $u_i\left(x_i,t\right)$ (Equations (1), and (2)).

8) Obtain the bubble velocity $U_i^{\left(n\right)}$ from $u_i\left(x_i,t\right)$ and ${\Delta }_M\left(x_i,{\sigma }_M\right)$ (Equations (1), (2), (4), and (7)).

9) Solve the next-time step bubble position from $U_i^{\left(n\right)}$ and $\Delta t$ .

10) Return to Step 2.

In the present paper, N was put to be 1, because we studied only a single bubble case, and the affix (n) will be deleted henceforth.

3. Correction of the FCM

FCM was originally developed for the fluid including solid particles. Therefore, the calculation results using the original version of FCM were first compared with the results by Takemura et al. [13] , where the bubble surface was like a solid sphere. The numerical results proved to agree well for small Reynolds numbers, but not very well for comparatively large Reynolds numbers. If the Reynolds number is about 10, the numerical value was about twice the experimental value. Then, the calculation results using the original version of FCM were compared with those obtained by Takemura et al. [10] [11] , where the bubble surface was slip. Note that the applicability of the original FCM was not clear in the calculation by Xu et al. [9] . The results were found to be in poor agreement with the experimental results, even for small Reynolds numbers. Even if the Reynolds number was about 1 or less, the difference was about 50%. Therefore, we proposed several corrections to the original FCM.

3.1. Correction of FCM for a No-Slip Bubble

First, the equation of $U_i$ in FCM (Equation (7)) was changed for a no-slip bubble considering the drag coefficient of a sphere over a wide range of the Reynolds number. It is known that the drag coefficient of a spherical solid particle, $D_n$ , is very accurately approximated by the following equation in the range of the bubble Reynolds number, $R{e}_{\infty }$ , from 0.01 to 20 [18] .

$D_n=6\text{π}\nu {\rho }_{f}R{U}_{\infty }\left(1+0.1315R{e}_{\infty }^{0.82-0.05{\log }_{10}R{e}_{\infty }}\right),$ (11)

where $R{e}_{\infty }=2R{U}_{\infty }/\nu$ , and $U_{\infty }$ is the terminal velocity of the particle in the quiescent fluid. Then, the equation for determining the velocity of a no-slip bubble is changed to

$U_i={\displaystyle \iiint \left[u_i\left(x_i,t\right){\Delta }_M\left(x_i-Y_i,{\sigma }_M\right)\cdot \frac{1}{1+0.1315R{e}^{0.82-0.05{\log }_{10}Re}}\right]\text{d}{x}_1\text{d}{x}_2\text{d}{x}_3},$ (12)

where $Re=2R\sqrt{{\left({\tilde{u}}_i-U_i\right)}^2}/\nu$ , and ${\tilde{u}}_i$ is the fluid velocity though which the bubble travels. The system of Equations (1), (2), (3), (6), and (12) instead of Equation (7) will be called the renormalized FCM (RFCM) henceforth in this paper.

3.2. Correction of FCM in FMT

Next, the reaction of a no-slip bubble to the fluid was considered in terms of the drag force instead of the acceleration term, $\text{d}{U}_i/\text{d}t$ , for FMT in Equation (6). The following equation is introduced,

$F_i=-\frac{4}{3}\text{π}{R}^3\left({\rho }_f-{\rho }_b\right)g_i+6\text{π}\nu {\rho }_{f}R\left(U_i-{\tilde{u}}_i\right)\left(1+0.1315R{e}_{\infty }^{0.82-0.05{\log }_{10}R{e}_{\infty }}\right),$ (13)

taking into consideration of the drag force formula (11). The terminal velocity, $U_{\infty }$ , in the infinite quiescent fluid is given by

$U_{\infty }=\frac{2}{9}\frac{g_i{R}^2}{\nu }\cdot \frac{1}{1+0.1315R{e}_{\infty }^{0.82-0.05{\log }_{10}R{e}_{\infty }}}\text{.}$ (14)

Therefore, $R{e}_{\infty }$ is obtained by solving the following equation.

$R{e}_{\infty }=\frac{4}{9}\frac{g_i{R}^3}{{\nu }^2}\cdot \frac{1}{1+0.1315R{e}_{\infty }^{0.82-0.05{\log }_{10}R{e}_{\infty }}}\text{.}$ (15)

Note that $R{e}_{\infty }$ is uniquely determined if R is specified.

The drag coefficient of a spherical slip bubble, $D_s$ , is known to be very accurately approximated by the following equation in the Reynolds number range from 0 to 100 [19] ,

$D_s=4\text{π}\nu {\rho }_{f}R{U}_{\infty }\left\{1+{\left[\frac{8}{R{e}_{\infty }}+0.5\left(1+\frac{3.315}{R{e}_{\infty }^{0.5}}\right)\right]}^{-1}\right\}\text{.}$ (16)

Then, the reaction from a slip bubble to the fluid was considered in terms of the drag force instead of the acceleration term $\text{d}{U}_i/\text{d}t$ for FMT in Equation (6) as

$F_i=-\frac{4}{3}\text{π}{R}^3\left({\rho }_f-{\rho }_b\right)g_i+4\text{π}\nu {\rho }_{f}R\left(U_i-{\tilde{u}}_i\right)\left\{1+{\left[\frac{8}{R{e}_{\infty }}+0.5\left(1+\frac{3.315}{R{e}_{\infty }^{0.5}}\right)\right]}^{-1}\right\}\text{.}$ (17)

Using the drag force formula (16), $U_{\infty }$ is given by

$U_{\infty }=\frac{1}{3}\frac{g_i{R}^2}{\nu }\cdot \frac{1}{1+{\left[\frac{8}{R{e}_{\infty }}+0.5\left(1+\frac{3.315}{R{e}_{\infty }^{0.5}}\right)\right]}^{-1}},$ (18)

where $R{e}_{\infty }$ is obtained by solving the following equation.

$R{e}_{\infty }=\frac{2}{3}\frac{g_i{R}^3}{{\nu }^2}\cdot \frac{1}{1+{\left[\frac{8}{R{e}_{\infty }}+0.5\left(1+\frac{3.315}{R{e}_{\infty }^{0.5}}\right)\right]}^{-1}}.$ (19)

Now, $R{e}_{\infty }$ is uniquely determined if R is specified. The system of Equations (1), (2), (3), (7), and (13) instead of Equation (6) will be called the modified FCM-n (MFCM-n) for a no-slip bubble henceforth in this paper, while the system of Equations (1), (2), (3), (7), and (17) instead of Equation (6) will be called the modified FCM-s (MFCM-s) for a slip bubble.

4. Numerical Scheme

The calculation domain of the present study was a cube, the length of each side of which was 0.06 m, as shown in Figure 1. Note that the gravitational acceleration is in the negative z-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 $n_x=256$ . The finite-difference method was used with equally spaced staggered grids, where the number of grids $n_y=256$ in the y-direction and $n_z=256$ in the z-direction. Note that the number of grids within the bubble radius in this study was about 3 at most.

The test section used in the experiment by Takemura et al. [10] [13] had dimensions of 0.5 m in the x-direction and 0.06 m in both the y- and z-directions. The numerical calculation of the present study, on the other hand, has been performed at a smaller length of x in order to reduce the calculation load. In spite of this, it was possible to compare the experimental data and the results of

numerical calculations, because the measurement in the experiment was conducted in a relatively small region of x where the rising bubble reached a nearly asymptotic state.

Figure 2 shows comparisons of the test sections in the y-z plane and of the position of the bubble in the experiment performed by Takemura et al. [10] [13] and in the present numerical study. In this study, the gravitational acceleration was 9.81 m/s² and the kinematic viscosity, $\nu$ , and the bubble radius, R, were taken to be equal to be the experimental conditions. L shows the distance from the bubble center to the nearest vertical wall.

Because of the periodic boundary condition in the x-direction of the present study, a pressure gradient opposite to the gravitational acceleration was needed to keep the flow in the calculation domain quiescent as a whole. A bubble moving in the calculation region was accelerated from the initial stationary state by the buoyancy force. The initial position of the bubble center was put at the point $\left(0,0,0.03-L_0\right)$ , where $L_0$ , being the initial value of L, was put to be 2R.

The Navier-Stokes equations were solved by the simplified marker and cell method (SMAC). The fourth-order central difference scheme was used for spatial differentiation, the Adams-Bashforth and the Crank-Nicolson method were used for the advection and viscous terms, respectively, in the time advancement. In order to solve the pressure equation, the bi-conjugate gradient stabilization method (Bi-CGSTAB) was used.

5. Results of the Numerical Calculation and Comparison with the Experimental Results

In order to study the applicability of the proposed formulas, RFCM, MFCM-n, and MFCM-s, comparison with the experimental data for the trajectory of a rising bubble by Takemura et al. [10] [11] [13] was performed.

5.1. Experimental Results for Rising No-Slip and Slip Bubbles

Here, an outline of the experimental study by Takemura et al. [10] [11] [13] is given. Takemura et al. obtained the lift force by observing the motion of a bubble, but the horizontal (z-direction) velocity data was not given in the paper. In this study, thus, the lift force was calculated from the motion of the bubble

obtained by the present simulations and compared with the experimental results.

Takemura et al. introduced the wall distance Reynolds number, $R{e}_L$ , defined in terms of L as

$R{e}_L=\frac{U_{\infty }L}{\nu },$ (20)

and analyzed the force acting upon bubbles moving parallel to the wall under the assumption of the Oseen approximation when a wall exists in the infinite fluid. They modified the formula for the lift force obtained by Vasseur and Cox [12] considering the motion in the oblique direction to the wall. They obtained $I_{D-n}$ and $I_{D-s}$ , the inverse Fourier transform of the horizontal (z-direction) component of the drag force on no-slip and slip bubbles, as in Equations (21) and (22), respectively:

$I_{D-n}\cdot \frac{U_3}{U_1}=\frac{3U_3}{4U_1\text{π}R{e}_L^2}{\displaystyle {\int }_0^{\infty }{\displaystyle {\int }_0^{\text{2π}}\frac{{\xi }^2}{i\cos \phi }\left[\left(\frac{{\text{e}}^{-2\xi }}{\xi }-\frac{{\text{e}}^{-2\chi }}{\chi }\right)+\frac{2{\left({\text{e}}^{-\xi }-{\text{e}}^{-\chi }\right)}^2}{\chi -\xi }\right]\text{d}\xi \text{d}\phi }},$ (21)

$I_{D-s}\cdot \frac{U_3}{U_1}=\frac{U_3}{2U_1\text{π}R{e}_L^2}{\displaystyle {\int }_0^{\infty }{\displaystyle {\int }_0^{\text{2π}}\frac{{\xi }^2}{i\cos \phi }\left[\left(\frac{{\text{e}}^{-2\xi }}{\xi }-\frac{{\text{e}}^{-2\chi }}{\chi }\right)+\frac{2{\left({\text{e}}^{-\xi }-{\text{e}}^{-\chi }\right)}^2}{\chi -\xi }\right]\text{d}\xi \text{d}\phi }},$ (22)

where $U_1$ is the bubble velocity in the vertical direction (x-direction) and $U_3$ is that in the horizontal direction (z-direction), $\chi$ is defined by the following equations.

$\chi =R{e}_L\eta ,{\eta }^2=k_1^2+k_2^2+i{k}_1,k_1=\frac{\xi \cos \phi }{R{e}_L},k_2=\frac{\xi \sin \phi }{R{e}_L},$

where i denotes an imaginary unit. Takemura et al. [10] [11] [13] assumed the following approximate horizontal force balance, using the drag force equation of [18] for a no-slip bubble and that for a slip bubble [19] .

$\begin{array}{l}3\text{π}\nu {\rho }_{f}R{U}_{1}R{e}_{\infty }{I}_{L-n}\\ =6\text{π}\nu {\rho }_{f}R{U}_3\left\{1+0.1315R{e}_{\infty }^{0.82-0.05{\log }_{10}R{e}_{\infty }}\right\}+3\text{π}\nu {\rho }_{f}R{U}_{3}R{e}_{\infty }{I}_{D-n},\end{array}$ (23)

$\begin{array}{l}2\text{π}\nu {\rho }_{f}R{U}_{1}R{e}_{\infty }{I}_{L-s}\\ =4\text{π}\nu {\rho }_{f}R{U}_3\left\{1+{\left[\frac{8}{R{e}_{\infty }}+0.5\left(1+\frac{3.315}{R{e}_{\infty }^{0.5}}\right)\right]}^{-1}\right\}+2\text{π}\nu {\rho }_{f}R{U}_{3}R{e}_{\infty }{I}_{D-s},\end{array}$ (24)

where $I_{L-n}$ and $I_{L-s}$ are the inverse Fourier transform of the horizontal (z-direction) components of the lift force acting on a no-slip bubble and a slip bubble, respectively. Finally, $I_{L-n}$ and $I_{L-s}$ are given by

$I_{L-n}=\frac{2}{R{e}_{\infty }}\cdot \frac{U_3}{U_1}\left(1+0.1315R{e}_{\infty }^{0.82-0.05{\log }_{10}R{e}_{\infty }}+\frac{R}{L}R{e}_L{I}_{D-n}\right),$ (25)

and

$I_{L-s}=\frac{2}{R{e}_{\infty }}\cdot \frac{U_3}{U_1}\left\{1+{\left[\frac{8}{R{e}_{\infty }}+\frac{1}{2}\left(1+\frac{3.315}{R{e}_{\infty }^{0.5}}\right)\right]}^{-1}+\frac{R}{L}R{e}_L{I}_{D-s}\right\}.$ (26)

Because the data for $I_{L-n}$ or $I_{L-s}$ were given by Takemura et al. [10] [11] [13] , the values of $I_{L-n}$ and $I_{L-s}$ were calculated from $U_1$ and $U_3$ in the present calculation and compared with their experimental results. Note that the drag components $I_{D-n}$ and $I_{D-s}$ can be neglected when $R{e}_L\ge 10$ [10] [11] [13] .

5.2. Results of RFCM and Comparison with Experimental Data for a No-Slip Bubble

5.2.1. Trajectory of a Rising No-Slip Bubble near the Wall

Figure 3 shows the calculated trajectory of a no-slip bubble rising near the wall, starting from the position (0, 0, 0.03 − 2R) at zero initial velocity for $R{e}_{\infty }=0.1$ , 0.5, 1, and 5. In Figure 3, the bubble for $R{e}_{\infty }=0.1$ is shown to rise almost vertically, but for $R{e}_{\infty }\ge 0.5$ , it rises in a direction that gradually moves away from the wall. This indicates that, for larger $R{e}_{\infty }$ , a relatively stronger lift force acts on the bubble in the direction away from the wall due to the influence of the wall. Although it is not shown here, it was found that a bubble finally rises in the vertical direction for all values of $R{e}_{\infty }$ if the simulation is performed over a sufficiently long time. The lift force due to the influence of the wall decreases as the bubble moves farther away.

5.2.2. Comparison with the Experimental Data of a No-Slip Bubble for Re_∞ > 0.1

Figure 4 shows the value of $I_{L-n}$ as a function of $R{e}_L\left(=L{U}_{\infty }/\nu \right)$ for $R{e}_{\infty }=0.9$ , 2.5, and 5. The black symbols indicate the experimental data of Takemura et al. [13] , and the color lines are $I_{L-n}$ values obtained by Equation

(25) based on the calculation results using RFCM in the present study. This figure also shows that the results by RFCM are in good agreement with the experimental data for $R{e}_{\infty }=0.9$ , 2.5, and 5.

Figure 5 shows $I_{L-n}$ versus $R{e}_L$ for $R{e}_{\infty }=7.7$ , 14.2, and 17.7. They show good agreement between the present numerical results and the experimental data, although the discrepancy is about 17% at most for $R{e}_{\infty }=7.7$ and about 15% at most for $R{e}_{\infty }=14.2$ . In the experiment of Takemura et al. [13] , the bubble surface becomes like a solid sphere by adding surfactant. Therefore, it is considered that incomplete formation of a solid-like surface may cause the discrepancy between the experiment and the present theory. Note that the difference between the experiment and the present theory is less than 5% for other Reynolds numbers, and it is concluded that RFCM is applicable for calculating multiphase flow with water and a no-slip bubble.

5.3. Results by MFCM-N and Comparison with the Results by RFCM and the Experimental Data for a No-Slip Bubble at Re_∞ = 17.7

In order to study the applicability of MFCM-n to the no-slip bubble motion, the calculation results by MFCM-n were compared with those by RFCM and the experimental data of Takemura et al. [13] for the case of $R{e}_{\infty }=17.7$ . Figure 6 shows the trajectories of a no-slip bubble rising near the wall starting from the position (0, 0, 0.03 − 2R) at zero initial velocity for $R{e}_{\infty }=17.7$ , as obtained by RFCM (a purple line) and MFCM-n (a red broken line). Both results are found to be very close. Figure 7 shows $I_{L-n}$ , the lift force acting on the bubble, as obtained by RFCM, MFCM-n, and the experiment for $R{e}_{\infty }=17.7$ . It is found that

MFCM-n is applicable for $R{e}_{\infty }=17.7$ , and probably for $R{e}_{\infty }\le 17.7$ as RFCM.

5.4. Results by MFCM-S and Comparison with the Experimental Data for a Slip Bubble

Although RFCM was applied to predicting the trajectory of a no-slip bubble, it proved to be inapplicable to a slip bubble. Even if the Reynolds number was about 1 or less, the difference was about by 50%. Because the boundary conditions for no-slip and slip bubbles are different, the drag force is larger for a

no-slip bubble than that for a slip bubble and the terminal velocity of a slip bubble is at least 1.5 times larger than that for a no-slip bubble. Therefore, the lift force for a slip bubble is larger than that for a no-slip bubble if the Reynolds number is same (Equations (15), (19), (25), and (26)). The MFCM-s was applied to trace a slip bubble trajectory. Figure 8 shows the values of $I_{L-s}$ as a function of $R{e}_L$ for $R{e}_{\infty }=2.1$ and 4.2, and Figure 9 shows that for $R{e}_{\infty }=5.5$ , 7.3, and 15, where the black symbols are the experimental data of Takemura et al. [10] [11] . The results by MFCM-s are found to be in good agreement with the experimental data for all $R{e}_{\infty }$ except for $R{e}_{\infty }=15$ . The details of MFCM-s were already discussed in the preceding paper with the mechanism of the production of the lift force [20] .

6. Conclusions

Three new methods for simulating multiphase flow consisting of water and air bubbles were proposed based on FCM in order to trace the trajectories of a spherical bubble with a slip or no-slip surface rising near the wall in the quiescent fluid. The calculation results were compared with the experimental data by Takemura et al. [10] [11] [13] over the Reynolds number from 0.9 to 17.7.

The first method is called RFCM, where the induced bubble velocity was renormalized according to the formula of the drag force for a no-slip bubble according to [18] . This method was successfully applied to trace the trajectory of a bubble and gave the values of the lift force on a no-slip bubble for $0.9\le R{e}_{\infty }\le 17.7$ , where $R{e}_{\infty }$ is the bubble Reynolds number determined by its radius.

The second method is called MFCM-n, where FMT was modified according to the drag force equation given in [18] . This method was successfully applied to obtain the values of the lift force on a no-slip bubble for $R{e}_{\infty }=17.7$ . The results agreed with those found by RFCM.

The third method is called MFCM-s, where FMT was modified according to the drag force equation given in [19] for a slip bubble. This method was successfully applied to obtain the values of the lift force on a slip bubble for $2.1\le R{e}_{\infty }<15$ . However, calculation results using RFCM for a slip bubble could not reproduce the experimental data.

In the present paper, it is found that the three newly proposed methods can correctly trace the trajectory of a bubble near a vertical wall for relatively large Reynolds numbers, regardless of whether the bubble has a slip or no-slip surface. In future study, we investigate the effects of the value of the Gaussian distribution’s width, $\sigma$ , in FMT and the modification of $G_{ij}^{\left(n\right)}$ in FDT and develop simulation methods to deal with the fluid which includes many bubbles.

Acknowledgements

S. Yanase would like to give his cordial thanks to the Ministry of Education, Culture, Sports, Science and Technology for the financial support through the Grant-in-Aids for Scientific Research No. 15K05814 and 15H0419919.