Abstract
The motion of multiple bubbles under gravity in two dimensions is numerically studied through the lattice Boltzmann method for the Eotvos number ranging from 1 to 12. Two kinds of initial arrangement are taken into account: vertical and horizontal arrangement. In both cases the effects of Eotvos number on the bubble coalescence and rising velocity are investigated. For the vertical arrangement, it has been found that the coalescence pattern is similar. The first coalescence always takes place between the two uppermost bubbles. And the last coalescence always takes place between the coalesced bubble and the bottommost bubble. For four bubbles in a horizontal arrangement, the outermost bubbles travel into the wake of the middle bubbles in all cases, which allows the bubbles to coalesce. The coalescence pattern is more complex for the case of eight bubbles, which strongly depends on the Eotvos number.
1. Introduction
Bubble motion is one of the most common gas-liquid flow phenomena and plays an important role in many industrial applications, such as cavitation in fluid machinery, nucleate boiling in reactors, and condenser/evaporator. In medical ultrasound imaging, small encapsulated bubbles called contrast agent are used to enhance the contrast. In thermal inkjet printing, vapor bubbles are used as actuators. They are occasionally used in the microfluidics applications as actuators [1]. The motion of multiple bubbles under gravity is complex due to bubble deformation, coalescence, and breakup. Understanding the dynamic interaction between bubbles is an important aspect of the design and operation of many industrial applications.
A number of investigations of the bubble motion in liquid have been conducted in the past [2–5] due to its scientific and engineering importance. Many numerical methods have been established to simulate the motion of gas bubbles in liquid, such as VOF (volume of fluid) [6] and level set method [7]. In recent decades, the lattice Boltzmann method (LBM) has proved to be a powerful numerical scheme for the simulation of multiphase flow which is based on mesoscopic particle dynamics. Its kinetic nature can provide many of the advantages of molecular dynamics, such as the introduction of repulsive interaction between particles without any boundary conditions for the interface, which makes it more useful than other more conventional methods for numerical modeling of multiphase fluid systems [8]. Several kinds of lattice Boltzmann model for simulating multiphase fluid have been established and applied to the simulations of gas bubbles under gravity successfully. These models include potential method proposed by Shan and Chen [9], color method proposed by Rothman and Keller [10], and free energy method proposed by Swift et al. [11] and improved by Zheng et al. [12]. Recently, Huang et al. [13] evaluated these models’ performance in modeling two-dimensional immiscible two-phase flow in porous media on the pore scale. The free energy method [12] was proved to be a good tool for the study of two-phase flows with high viscosity ratios and high density ratio. Takada et al. [8] used the LBM which was based on the free energy model [11] to simulate single bubble and two bubbles rising under gravity. Their results [8] of single bubble are consistent with those determined by VOF method. They [8] also simulated two bubbles rising in a circular tube full of stationary fluid. It has been found [8] that the two bubbles approach due to the wake formation and then coalesce into a single bubble eventually. Sankaranarayanan et al. [14] simulated the rise behavior of a single bubble in a periodic box by using an implicit formulation of LBM under the conditions of 1 < Eo (Eotvos number) < 10 and (Reynolds number) > 100. The effect of volume fraction on the rise characteristics was analyzed by changing the size of the box relative to that of the bubble. Recently, Amaya-Bower and Lee [15] presented a comprehensive study of the dynamics of a single bubble rising under gravitational force by using lattice Boltzmann method based on the Cahn-Hilliard diffuse interface approach. They [15] analyzed the dependence of terminal bubble shape and velocity on density ratio, viscosity ratio, surface tension, interface thickness, and domain size. In particular, their results [15] showed that the influence of viscosity ratio on the terminal velocity and shape of the bubble is significant, while the density ratio has little effect. Gupta and Kumar [16] used the lattice Boltzmann method based on potential method [9] to study the motion of multiple bubbles under gravity. Two kinds of initial arrangement were considered which are in-line arrangement and staggered arrangement in a vertical orientation, respectively. Their results [16] showed that as the Eotvos number increases, the uppermost bubble deforms the most and the bubbles coalesce eventually due to the wake behind the leading bubble. They [16] also showed the process of bubble coalescence for staggered bubbles in channels, in which lift forces come into play due to the presence of walls. More recently, Yu et al. [17] presented an adaptive lattice Boltzmann method to simulate the interaction between a pair of bubbles with spherical or ellipsoidal shapes under gravity. They observed both attractive and repulsive interactions between bubbles which depend on the relative position and the Reynolds number. They [17] also simulated a group of 14 bubbles and analyzed the effects of the bubble shape and Reynolds number on the spatial distribution of the bubbles.
As shown by Gupta and Kumar [16] and Yu et al. [17], in comparison with the case of single bubble, the motion of multiple bubbles under gravity could be much more complex due to hydrodynamic interactions between bubbles. Besides, the bubble coalescence and break-up could take place occasionally which have great effect on the motion of multiple bubbles. Furthermore, as far as we know, a detailed numerical study of the influence of bubble collision and coalescence on the rising velocity has not been undertaken. Thus, it is important to focus on the fundamental understanding of the bubble collision and coalescence when multiple bubbles are rising under gravity. To fulfil this task, the LBM proposed by Zheng et al. [12] is adopted in this work to study the rise behavior of multiple bubbles which are initially placed in in-line arrangement. The effects of Eotvos number which ranges from 1 to 12 are taken into account. This study will evaluate the coalescence pattern and rising velocity of the bubbles which not only depend on the computational parameters but also depend on the initial arrangements. In addition, the terminal velocity is compared with the result of single bubble under the same conditions to illustrate the influence of multiple bubbles motion.
2. Numerical Method
The discrete lattice Boltzmann equations under external forces for the continuity and momentum equations are given by where is the density distribution function at the th microscopic velocity , is the equilibrium distribution function, is the external force (which is gravity in this work), is the time step, is the relaxation time, is the speed of sound, and are the weights related to the D2Q9 model. According to Zheng et al. [12], is the order parameter that is responsible for the gas-liquid interface. is the chemical potential, which is defined in the following. The macroscopic variables and are determined by the distribution function as follows:The variables and are defined aswhere and are the densities of liquid and gas phase, respectively. The lattice model of D2Q9 is used here and the equilibrium distribution functions are defined asAccording to Zheng et al. [12], the coefficients are defined as The discrete lattice Boltzmann equations for the interface capture equation are given by [12]where is a constant coefficient, which is determined byThe macroscopic variable, that is, the order parameter , is given byAccording to Zheng et al. [12], the chemical potential is computed usingwhere and and are the densities of the liquid phase and the gas phase, respectively. and are parameters related to the thickness of the interface layer and the surface tension coefficient , which are expressed asAccording to Zheng et al. [12] the lattice model of D2Q5 is used for and the equilibrium distribution functions are [12]For simplicity the D2Q5 is used here and the coefficients are taken asThe parameter is used to control the mobility. By Chapman-Enskog analysis, the Navier-Stokes equations and an interface evolution equation can be obtained by (1) and (6)where is the mobility, given by . The viscosity is .
3. Validation
Three cases are carried out to validate the present computational code. The first is the Laplace law, which is given (for the two-dimensional case) bywhere is the pressure jump across the interface and is the bubble radius. In this work, the Laplace law is validated by calculating the pressure jump while varying the bubble radius from to . Other parameters are fixed at , , , and . The interface layer thickness is set to be for and for . The relaxation times are and . The results are shown in Figure 1(a), which shows good agreement between the numerical results and the analytical solution of (14).
(a)
(b)
The second validation case is the profile along the normal direction of the interface, which can be written as [12]where and are the coordinates of the center of the bubble. Figure 1(b) shows the order parameter along the -axis from the center point of the gas bubble. Three density ratio levels are taken into account; that is, , 500, and 100. The mobility coefficient and the surface tension coefficient are taken to be and for and 500, and and for . The bubble radius is fixed at and the interface layer thickness is . The other parameters are the same as in the first case above. The numerical results agree well with the corresponding analytical solution of (15), as shown in Figure 1(b). In all cases, the order parameter captures the interface between the two phases accurately.
Finally, a single bubble rising under gravity is simulated to compare with previous results. The bubble velocity is computed throughIn fact the result obtained by (16) is an averaged velocity of all bubbles. The parameters are taken the same as those used by Takada et al. [8] and Zheng et al. [12] which are , , and . Besides, the Eotvos number is introduced Comparisons of the terminal rising velocity for different Eotvos numbers are shown in Table 1. It is observed that the present results agree with the results obtained by Takada et al. [8] (VOF method) and by Zheng et al. [12] (LBM method). The shapes of the rising bubble under gravity for different Eotvos numbers are also shown in Figure 2. It is observed that as the Eotvos number increases, the deformation of the bubble is becoming more and more remarkable. This is due to the fact that the influence of the surface tension force which is to resist the deformation of the bubble is insignificant when the Eotvos number is large. In terms of the shapes of bubble, the present results agree with the results obtained by Takada et al. [8] (VOF method).
(a) Eo = 5
(b) Eo = 10
(c) Eo = 20
(d) Eo = 40
4. Results and Discussion
In this section multiple bubbles motion for both horizontal arrangement and vertical arrangement is investigated in detail. The parameters are fixed at , , , , , and . The computational domain is 250 × 1500 and stationary walls are applied for all boundaries.
Figure 3 shows the time evolution of the rising velocity of two bubbles for . The two bubbles are initially placed in a vertical orientation. The initial bubble distance is 2.5. As shown in the figure, after an initial period of evolution the velocity becomes constant which indicates that the merged bubble is rising at a constant velocity due to the balance between buoyancy and drag forces. For better understanding the contours of order parameter which visualize the process of bubble collision and coalescence are also shown in Figure 3. It is observed that an oscillation occurs in the bubble velocity before reaching the steady state. Obviously, this oscillation corresponds to the bubble coalescence process. Furthermore, as shown in Figure 3 the terminal velocity is smaller than the velocity before the two bubbles collide and coalesce.
The effects of the Eotvos number on the evolution of bubble velocity are shown in Figure 4. One can observe that the bubble velocity as well as the magnitude of its oscillation increases with the Eotvos number. Moreover, it has been found that in all cases () the two bubbles would collide and coalesce. The lower bubble always rises at a higher speed than the upper one because it is located at the wake behind the upper bubble. This leads to a smaller drag force for the lower bubble. Since the relative velocity of the bubbles is nonzero, the distance between the two bubbles keeps decreasing with time. Eventually the bubbles come to contact and merge and soon after form a larger bubble with twice the volume as the initial bubble.
The motion of four bubbles under gravity is also studied in this work. The initial distance between bubbles is 2.5. Figure 5 shows the contours of order parameter for at different time to display the entire evolution with time. As shown in Figure 5(b), the two uppermost bubbles begin to coalesce firstly after some time. Then the newly merged bubble and the bubble next to it begin to collide and coalesce again, as shown in Figure 5(c). Finally, the coalescence process occurs between the merged bubble and the bottommost bubble, which leads to a spherical-cap shaped bubble with an interface inside it, as shown in Figures 5(e) and 5(f). This is due to high inertia at the time of coalescence, so the liquid is unable to squeeze out completely in time through the narrow gap between the two bubbles, thus getting trapped inside the gas phase. As time goes on, the liquid “pops” out from the top surface of the bubble, and the interface eventually becomes a part of bubble boundary, as shown in Figures 5(g) and 5(h). Similar phenomenon has also been reported by Takada et al. [8] and Gupta and Kumar [16]. In the work of Gupta and Kumar [16], the liquid “pops” out from the bottom surface of the bubble instead. Figure 6 shows the time history of velocity of the four bubbles. For better understanding of the effects of the coalescence on the velocity, the marks of “(a)–(h)” in Figure 6 are labelled which are corresponding to the times of Figures 5(a)–5(h). It is observed that there are several oscillations in the velocity before it becomes constant. Obviously, the oscillations labelled by (b), (c), and (e) in Figure 6 are caused by the bubble coalescence which takes place three times during the rising process of the four bubbles. The last oscillation labelled by (g) in Figure 6 is due to the interface inside the bubble moving to the top surface.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
Figures 7 and 8 show the time history of the rising velocity of four bubbles for and , respectively. Similarly, the bubble coalescence is displayed by showing the contours of order parameter at different times in the figures. Obviously, the magnitude of oscillations in the velocity is much larger than that of . It is also expected that the deformation of the merged bubble is much enhanced for larger Eotvos numbers. Furthermore, it can be observed that in both cases the uppermost bubble has the maximum shape deformation, due to the highest drag experienced by this bubble. Furthermore, it has been found that in all cases () the four bubbles merge into one bubble in a similar manner. In other words, the sequence of bubble coalescence among the four bubbles is similar to that of . The first bubble coalescence always takes place between the two uppermost bubbles. And the last coalescence takes place between the merged bubble and the bottommost bubble at all time. Similar observation is obtained even for larger initial bubble distances, as shown in Figure 9, which shows the merging process of four bubbles separated by center to center distance of 5 initially.
(a)
(b)
(c)
(d)
In order to present more insight into the rising process of multiple bubbles under gravity, the terminal velocities of two and four bubbles for different Eotvos numbers are shown in Figure 10. For the purpose of comparison, the results of a single bubble under the same conditions are also shown in the figure. It is unexpected that for very low Eotvos numbers such as , the terminal velocities of two and four bubbles are a little larger than those of a single bubble, as can be seen in Figure 10. The reason is not clear. However, for high Eotvos numbers such as the terminal velocities are similar for all cases, indicating that the influence of the wall boundaries on the terminal velocity is insignificant. Moreover, the bubble coalescence strongly affects the time evolution of bubble velocity but has little influence on the terminal velocity.
In this work the motion of the bubbles which are initially placed in a horizontal orientation is also studied. For the case of two bubbles, the bubble coalescence does not happen at all time for the Eotvos number that was concerned. Thereby the present work primarily focuses on the cases of four bubbles and eight bubbles.
Figure 11 shows the contours of order parameter of four bubbles for . Because of wall boundaries the two outermost bubbles are rising at a lower speed than that of the middle bubbles, as shown in Figure 11(b). Then the outermost bubbles are traveling into the wake left behind the middle bubbles, as shown in Figures 11(c) and 11(d). Eventually, two larger bubbles are formed due to the coalescence process which takes place at the same time for the left and right pairs of bubbles, respectively, as shown in Figures 11(e) and 11(f). The effects of Eotvos number on the motion of four bubbles are also investigated in this work. In all cases one can observe the similar process of bubble evolution. There are always two bubbles left eventually. For reason of space other results are not shown here.
(a)
(b)
(c)
(d)
(e)
(f)
However, things could be a little more complex for the cases of eight bubbles. Figure 12 shows the evolution of the rising bubbles under gravity for . As shown in Figures 12(b) and 12(c), firstly three pairs of bubbles are colliding nearly at the same time. Then the merged bubble in the middle is rising faster than other bubbles, as shown in Figures 12(d) and 12(e). After a while this bubble is leading. Meanwhile, the other two merged bubbles are becoming close enough to contact and coalesce again, as can be seen in Figure 12(f). Obviously, this bubble is located at the wake behind the leading one. Thereby after a period of trailing it catches up with the leading bubble at last and of course they merge into a larger bubble again, as shown in Figure 12(h). There are three bubbles left eventually in this case. However, as Eotvos number increases the situation is different, such as , as can be seen in Figure 13. At first only two pairs of outermost bubbles coalesce at the same time when rising under gravity, as shown in Figures 13(b) and 13(c). Due to the wall boundaries the two merged bubbles are behind other bubbles for a long time, as shown in Figures 13(d)–13(f). It can also be observed that before the coalescence takes place, the bubbles which are in two columns are rising in an oscillatory manner. Eventually, there are only two bubbles left after the coalescence process, which is much different from the case of .
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
It should be stated that numerical results have shown that for larger Eotvos numbers () the evolution is qualitatively similar to that of . However, a little difference has been found, as can be seen in Figure 14. This figure shows the rising pattern of eight bubbles for , which is the largest Eotvos number taken in this work. After the outermost bubbles coalescence they are not behind all the other bubbles instead. They are traveling into the wake left behind the uppermost bubbles, as can be observed from Figures 14(d)–14(f), which is different from the case of as shown in Figure 13. Eventually, two bubbles are formed after the bubble coalescence process, which is similar to the case of . However, as shown in Figure 14(h) there exists an interface inside the merged bubble, indicating that some liquid is trapped inside the bubble during the bubble coalescence process. This phenomenon was not observed for small Eotvos numbers.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
Figure 15 shows the terminal velocities of four and eight bubbles for different Eotvos numbers. For the case of eight bubbles, only the results of totally two merged bubbles are taken into account. For the purpose of comparing, the results of one bubble and two bubbles under the same configuration are shown in Figure 15 to indicate the effects of multiple bubbles motion on the rising velocity. Surprisingly, for the Eotvos number ranging from 1 to 12 () that was simulated in this work the terminal velocity is similar for the cases of two, four, and eight bubbles, indicating that the influence of the bubble coalescence on the terminal velocity is insignificant. However, these results are smaller than the results of the case of one bubble, as shown in Figure 15. This is mostly due to the influence of wall boundaries and the interactions between bubbles. Furthermore, as the Eotvos number increases the difference of terminal velocity between one bubble and multiple bubbles keeps increasing.
5. Conclusion
In this work the lattice Boltzmann method is used to simulate the motion of multiple bubbles under gravity. Benchmark studies have been conducted to validate the code. Roughly speaking, two kinds of initial bubble arrangements are taken into account in the simulations: vertical arrangement and horizontal arrangement. In both cases the focus is on the bubble coalescence pattern and the (averaged) rising velocity. The effects of the Eotvos number are studied in detail.
For vertical arrangement, the bubbles always coalesce in a similar manner for Eotvos number ranging from 1 to 12. The uppermost bubble deforms the most because of the maximum drag. In addition, after the coalescence a larger bubble is formed eventually. For more than two bubbles, the first coalescence takes place between the uppermost bubble and the bubble next to it at all time, and the last coalescence always takes place between the merged bubble and the bottommost bubble. The terminal velocity of two or four bubbles is a little larger than that of one bubble under the same conditions. The reason is not clear. However, the terminal velocity is similar for all cases when the Eotvos number is large.
Things are much different for the case of horizontal arrangement. It has been shown that in the case of four bubbles the outermost bubbles always travel into the wake of the middle bubbles, and the coalescence always takes place which leads to two larger bubbles. The situation is more complex for the case of eight bubbles. For very low Eotvos number the bubble coalescence takes place five times and three bubbles are left at last, which are one merged bubble and two initial bubbles. However, as Eotvos number increases there are two merged bubbles left resulting from the bubble coalescence which takes place six times. Furthermore, it has also been observed that the terminal velocity is similar for the cases of two, four, and eight bubbles. Comparison is also been made with the results of one bubble, indicating the presence of the well-known phenomenon of hindered rise for multiple bubbles.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work was supported by the National Program on Key Basic Research Project of China (no. 2011CB706501), the Major Program of the National Natural Science Foundation of China (no. 11132008), and the National Natural Science Foundation of China (no. 11272302).