Study on the Bubble Growth and Departure with A Lattice Boltzmann Method

For the growth and departure of bubbles from an orifice, a free energy lattice Boltzmann model is adopted to deal with this complex multiphase flow phenomenon. A virtual layer is set at the boundary of the flow domain to deal with the no-slip boundary condition. Effects of the viscosity, surface tension, gas inertial force and buoyancy on the characteristics of bubbles when they grow and departure from an orifice in quiescent liquid are studied. The releasing period and departure diameter of the bubble are influenced by the residual gas at the orifice, and the interaction between bubbles is taken into consideration. The relations between the releasing period or departure diameter and the gravity acceleration show fair agreements with previous numerical and theoretical results. And the influence of the gas outflow velocity on bubble formation is discussed as well. For the bubbles growing in cross-flow field, effects of the cross-flow speed and the gas outflow velocity on the bubble formation are discussed, which is related to the application in ship resistance reduction. And optimal choice of the ship speed and gas outflow velocity is studied. Cases in this paper also prove that this high density ratio LBM model has its flexibility and effectiveness on multiphase flow simulations.


Introduction
When the gas is released into the liquid through an orifice, a bubble forms. The growth and departure of the bubble from the orifice is a very common two-phase flow phenomenon in daily life and industrial applications, such as in beverage, pharmacy, marine engineering, environment and aerospace industry (Hepworth et al., 2003;Zhu et al., 2011;Li et al., 2016;Wei et al., 2018). In recent years, the bubble growth and departure is a hot research topic and the physical mechanism in this process is gradually revealed.
Works relating this topic can be classified into theoretical study, experimental study and numerical simulations. The theoretical predictions on the bubble departure size were given by Fritz (1935). The relation between the departure frequency and bubble size were studied by Zuber (1963) and Malenkov (1971). In the quiescent liquid, Jamialahmadi et al. (2001) conducted experiments to observe the bubble motion in a variety of solutions. Zhang and Tan (2003) carried out the experimental study of the bubble growth and departure in the cross flow. The actual physical process of the bubble formation is under the influence of a variety of forces. And phase transition usually happens. Therefore, there is not a perfect theoretical model which can describe this phenomenon accurately. Meanwhile, the cost in experiments is high and there are too many uncontrollable factors (Zhang et al., 2015a;Li et al., 2017). Hence, the numerical methods, which may overcome the above stated shortcomings, develop rapidly. The volume-of-fluid (VOF) method was used to simulate the growth and departure process of a bubble successfully by Islam et al. (2015) and Sunder and Tomar (2013). The Lagrangian smoothed particle hydrodynamics (SPH) method is an effective method to calculate multiphase flow (Zhang et al., 2015b;Cheng et al., 2017). The boundary element method (BEM) is also widely used in multiphase flow simulation (Zhang and Liu, 2015;Wu et al., 2017).
The lattice Boltzmann method (LBM) has developed rapidly in recent decades. The LBM has good computational stability, and it is easy for parallelization. These advantages prompt the LBM to be widely adopted in multiphase flow simulations (Nourgaliev et al., 2003;Wang et al., 2015;Chen et al., 2019a). For the bubble evolved under the gravity effect, Takada et al. (2001) compared their results from the LBM with that of the VOF method and Lee and Lin (2003) devised a pressure evolution equation of the LBM and took the phase transition into consideration. The interface capturing technique was applied to the LBM for simulating multiphase flows with large density ratios, which succeeded in simulating the motion of the rising bubble (Zheng et al., 2006). The rising bubbles passing through obstacles were simulated using the LBM by Alizadeh et al. (2017). The interaction of multiple bubbles with different initial sizes was studied by Cheng et al. (2010).
The studies stated above indicate that the LBM is an appropriate tool for investigating the gas-liquid two-phase flow problem. In this work, we will study the bubble formation under different conditions to explore the mechanism of growth and departure with LBM. In Section 2, the free energy 2D LBM model considering the large density-ratio problem is briefly reviewed for completeness. And we introduce a virtual layer outside the real boundary to deal with the difference of the distribution function at the boundary. Then in Sections 3 and 4, the model is validated by some analytical results, and the growth and departure processes of bubbles in the quiescent liquid and in the cross flow are calculated. Some features of the bubble as well as the flow field under the influence of different parameters are shown and discussed. And the model is associated with the potential application of the drag reduction. Finally, the summary and conclusions are presented in Section 5.

Governing equations
In this work, the fluid domain is regarded as incompressible with a constant temperature. The Navier-Stokes equations are used to govern the motion of the fluid, and the interface capturing equation is the Cahn-Hilliard equation (Zheng et al., 2006). They can be described as follows: ∂ρ ∂t + ∇ · (ρu) = 0; (1) where and are the velocity and the viscosity of the fluid, respectively. is the pressure tensor, is the body force, is the mobility parameter (Cahn et al., 1996), which is taken to be constant, and is the chemical potential de-ρ ϕ rived from the free energy density function (Rowlinson and Widom, 1989). and (Zheng et al., 2006;Takada et al., 2001) are defined as: where and denote the density of fluid 1 and fluid 2, respectively.
The term is relevant to the surface tension force (Zheng et al., 2006), and the force can be written as the function of the chemical potential and the sound speed , The free energy density function is used to derive the chemical potential (Rowlinson and Widom, 1989;Kendon et al., 2001;Zheng et al., 2006), which can be expressed in a closed volume as: where is a constant used to express the intensity of the interaction energy between different phases; and are the order parameter and its state in equilibrium, respectively. is a coefficient in connection with the surface tension and the interface thickness and denotes the control volume. Through the minimum free energy principle (Zheng et al., 2006), the chemical potential function reads K κ where and are expressed as: where denotes the surface tension coefficient and is the thickness of the interface. Here, the equilibrium order parameter is .

Lattice Boltzmann method
The Navier−Stokes equation for fluid flow can be expressed by the discrete Boltzmann equation with single-relaxation time and a forcing term (He and Luo, 1997), as: where is the density distribution function for the momentum, denotes the lattice velocity vector, is the lattice time step, the superscript 'eq' stands for the equilibrium state and is the single relaxation time of the momentum dispersion equation. The surface tension and the gravity are included in the body force as (Guo et al., 2002): where is the weight and is the gravity.
70 CHEN Guo-qing et al. China Ocean Eng., 2020, Vol. 34, No. 1, P. 69-79 c i The nine-velocity (D2Q9) model is used in the present work. The discrete lattice velocity and the weight w i are where the lattice speed with being denoted as the lattice spacing and being denoted as the lattice time step.
The equilibrium distribution function (Zheng et al., 2006) is described as: α β x y A i where and denote the and direction, respectively. And the coefficient reads The Chapman-Enskog expansion is performed on Eq. (10), which corresponds to the Navier-Stokes equation. Eq.
(10) has second order accuracy (He and Luo, 1997), which is enough for the calculations in this work. And the corresponding viscosity coefficient is expressed as: ρ The local macroscopic density and the velocity u are (Zheng et al., 2006) The improved interface capture method can be expressed by a modified lattice Boltzmann equation as (Zheng et al., 2006;Chen et al., 2019b): where and are the distribution function and its corresponding equilibrium state, respectively. is the single relaxation time, and is the lattice velocity vector. w j The five-velocity (D2Q5) model is proved accurate enough to discretize Eq. (19) (Zheng et al., 2006). The velocity vector c j and weight factor are The equilibrium distribution function (Zheng et al., 2006) is defined as: where is a parameter governing the mobility as: ϕ The order parameter is the summation of the distribution functions When calculating the terms and , Zheng et al. (2006) used the three-point and five-point central difference schemes respectively, in a 2D flow. And Cheng et al. (2010) adopted the eighteen-point central difference scheme for a 3D flow, and they both achieved satisfactory results. Nevertheless, the periodic boundary condition is used in their models, whilst the no-slip boundary is needed in our work for some cases. Hence, when dealing with the derivatives of parameters at the boundaries, a virtual layer of particles is added outside the real boundary, as shown in Fig. 1. on the virtual layer (−1) is extrapolated as (Chen et al., 1996) which has a second order accuracy and maintains the overall accuracy of lattice Boltzmann method.

General features
The lattice units are used in our work unless otherwise specified. The periodic boundary condition is applied to the left and right walls of the domain whilst the no-slip and open boundary conditions are used to the bottom and top boundaries, respectively. The density of fluid 1 and fluid 2 respectively are and , and the relaxation times are and (Zheng et al., 2006). The govern- CHEN Guo-qing et al. China Ocean Eng., 2020, Vol. 34, No.  ing parameter of the mobility is , the thickness of the surface layer is and the surface tension is . The interface capturing ability and Laplace law are firstly validated. The order parameters near the interface and the 2D Laplace law can be expressed as follows (Jacqmin, 1999;Zheng et al., 2006): and are coordinates of the bubble center, is the bubble radius, and denotes the pressure difference of interface. The results are shown in Fig. 2 and Fig. 3, indicating that the numerical results agree well with the theoretical ones, so the model has good interface capture ability and satisfies the Laplace law.
Whereafter, as shown in Fig. 4, we check the bubble growth and departure in a domain with the radius of the orifice for convergence analysis. Three different resolutions are chosen, namely , and . The result reveals that our model is convergent as the resolution increases.
We then focus on the general feature of bubbles growth Bo Re and departure in a quiescent liquid, as presented in Fig. 5 The size of the computational domain is . The center of the orifice is set at (31, 0) with the radius . The Bond number is and the Reynolds number is . The gas enters into the liquid with a velocity of . Growth and departure of the first and the second bubbles are recorded in the figure. The Bond number and Reynolds number are defined as: g where is the gravity acceleration.
The streamlines are depicted in Fig. 5, the influences of the bubble growth and departure on the flow field can be observed. At the beginning, the surface tension attaches the bubble to the orifice, while the buoyancy and the gas inertia force promote the growth and expansion of the bubble. And the viscous force is not obvious at this time, due to the small bubble velocity as well as the velocity gradient. The bubble growth induces two small vortices formed on both sides at the base of the bubble, as shown in Fig. 5a. When the bubble reaches a relatively large volume, the neck appears, which is labeled as the end of the bubble growth stage and the beginning of the departure stage. At this moment, the velocity near the neck region gets larger and so is the velocity gradient, hence the viscous force is obvious and it affects the pinch off of the bubble with surface tension. The forces make the flow near the orifice much more complicated and induce two vortices to form near the neck of the bubble, as demonstrated in Fig. 5b. Then, the bubble further grows and rises under the influence of buoyancy with the neck shrinking continuously. Soon after, the bubble breaks from the orifice. Meanwhile, small vortices near the orifice disappear and those near the neck develop into two big vortices and move with the rising bubble during this process. In the end, the bubble gradually becomes round due to the surface tension, as shown in Figs. 5c and 5d.  The second bubble forms when the first one rises up under the influence of the buoyancy. The disturbance to the flow field from the first bubble results in relatively large velocity and velocity gradient near the orifice, which makes the viscous effect become obvious at this area. Two small vortices are induced at the top of the second bubble, as shown in Fig. 5d. Every two vortices rotate in the opposite directions. Since the velocity of liquid at this region is upwards, this kind of viscosity presents some kind of pulling effect, which helps the formation of the second bubble. Then the second bubble continues to grow under the buoyancy and the inertia of gas. Evolution of the vortices is complex due to the viscosity. From Figs. 5d and 5e, we find that these two small vortices move down along both sides of the second bubble. Then the two small vortices move outwards with the expansion of the bubble and new vortices form on the top of the second bubble as shown in Fig. 5f. As the first bubble has gone far away, its influence on the second one is weak, see Figs. 5f and 5g. Then the necking of the second bubble begins. The dynamics of the second bubble is very similar to that of the first one, if comparing Fig. 5b and Fig. 5h. The growth and departure process of the following bubbles presents almost the same features as the second bubble, meaning that the process can be regarded as periodic henceforth.
By studying the growth and departure characteristics of the first and second bubbles, we find that at different stages, the bubble dynamics is dominated by different forces. The evolution of the bubble during its growth and departure is affected by the buoyancy, surface tension, gas inertia force and viscous force together. At the growth stage, the buoyancy and the gas inertia promote the growth and expansion of the bubble, the effect of the viscous force is not obvious, and the surface tension attaches the bubble to the orifice. The departure stage of the bubble starts when the neck appears. At this moment, the necking speed is governed by the viscous force and surface tension. The viscous force impedes the pinch off of the bubble while the surface tension accelerates the necking process. The velocity and the velocity gradient of the bubble at the necking region are relatively large during this stage, resulting in more obvious viscous effect. However, the surface tension is much stronger. Therefore, the departure happens finally.
What is different for the first bubble compared with its followers is that during its growth stage, the velocity as well as the velocity gradient of the liquid above the bubble are small, hence the viscous force there is weak. Whilst this viscous effect is considerable for the second bubble at the early stage of its growth, which helps pull the bubble out of the orifice. Meanwhile, there is still some residual gas near the orifice as shown in Fig. 3c, which helps the growth of the following bubbles. Hence, from the second bubble, it takes less time to pinch off. The generating time of the second bubble is chosen as the releasing period, which is defined as the time that it takes from the beginning of the formation until pinched off from the orifice, in this paper. Then, we focus on the bubble diameter and releasing period of the first and second bubble. The releasing period of a bubble is denoted as . And denotes the departure diameter of bubble . The subscripts '1' and '2' represent the first bubble and second bubble, respectively. The results are depicted in Fig. 6. The period of the first bubble is 20% longer than that of the second bubble for the gas outflow velocity . And this difference increases linearly with the increasing . So in the following cases, the releasing period of the bubbles is defined as that of the second bubble. Nevertheless, the difference of the diameters for these two bubbles fluctuates near zero no matter how large the gas velocity is, due to the existence of the residual gas near the orifice, which occupies a part of the volume of the gas entering the liquid for the first bubble. Therefore, we can ignore the change of the bubble size.
g v 3.2 Effect of the gravity acceleration and the gas outflow velocity on the departure diameter and the releasing period Relation of the gravity acceleration and the departure diameter for the bubble detaching from a horizontal plate has been studied in previous literature. Fritz (1935) found the relation of . In numerical simulations, the relation of was achieved by fitting the results of bubble detachment with small liquid-gas density ratio (Yang et al., 2001), and of boiling with high density radio was obtained by Sun et al. (2013). For bubble growth and departure in a quiescent liquid in this paper, the relation is depicted in Fig. 7. A relation of is achieved through fitting the data.
T D b g The releasing period also has a relation with and . In Zuber's paper (Zuber, 1963), it reads Substituting the relation stated above into Eq. (30), one may get . Through changing the g T g g T ∼ g −0.681 gravity acceleration , a fitting curve between and is obtained in Fig. 8. The bubble releasing period decreases with the increase of gravity acceleration as , which agrees well with the prediction of Eq. (30).
From Figs. 7 and 8, we find that and as well as and all follow exponential decaying relationships. The viscous force impedes the pinch off of the bubble whilst the surface tension promotes the necking. The buoyancy increases with the gravity acceleration , which induces larger velocity and velocity gradient in the liquid near the neck of the bubble, and finally resulting in more obvious viscous effect. However, influence of the surface tension is even stronger. So the larger g speeds up the process of bubble growth and departure, making the departure time and diameter smaller . In addition, these two forces make the absolute values of the slope of the fitting curves decrease as g increases.
Then we study the effect of the gas outflow velocity on the bubble departure diameter and the releasing period. The bubble forms faster with larger volume as the gas outflow velocity increases, as shown in Figs. 9 and 10. Fitting curves of the data show the relations and . Hence and also fall in the exponential relations with . And the slope of the fitted curves decrease v Fig. 6. Differences of the diameter and releasing period between the first and second bubbles with respect to the gas outflow velocity , and the relevant fitting lines in a quiescent liquid.  Fig. 7. Relation between the bubble departure diameter and the gravity acceleration with the relevant fitting curve, for the growth and departure of bubbles in a quiescent liquid.
T g Fig. 8. Relation between the releasing period and the gravity acceleration with the relevant fitting curve, for the growth and departure of bubbles in a quiescent liquid. v with increase of , which is mainly due to the effect of the gravity. When the volume of bubble is big, the effect of gravity is obvious and the inertia force of gas is weak. In this section, a stationary cross flow is considered preexisted in the domain. All boundary conditions and other parameters remain the same as the case in Fig. 5. The periodic boundary condition guarantees that the velocity of fluid is constant at the beginning. When the fluid flows across the bubble, the fluid quickly returns to its initial state due to the large density difference. Meanwhile, the size of the computational domain in x direction is large enough to guarantee that the first bubble has no influence on the second one, so we may postulate that the periodic boundary has no effect on the flow field near bubbles. When the gas is released from an orifice on the bottom boundary, the size of the computational domain is , the velocity of the cross flow is . The evolution of the bubble and the flow field are shown in Fig. 11.
In Fig. 11, the bubble loses its symmetry due to the drag force of the cross flow from the beginning. The bubble shape and the distribution and evolution of the vortex are consistent with the results in Sun et al. (2013). The drag force promotes the bubble growing rightwards and the surface tension attaches the bubble to the orifice, which induces an unsymmetrical small bubble adheres to the orifice. The bubble grows and expands continuously under the ef-T v Fig. 10. Relation between the releasing period and the gas outflow velocity with the relevant fitting curve, for the growth and departure of bubbles in a quiescent liquid.  Fig. 9. Relation between the bubble diameter and the gas outflow velocity with the relevant fitting curve, for the growth and departure of bubbles in a quiescent liquid. fect of the drag force, buoyancy and gas inertia force. When the bubble is elongated, the bubble blocks the flow field to the right side of the bubble. Hence, the velocity gradient appears there, resulting in an elliptical vortex which attaches to the bubble, as shown in Fig. 11b. Moreover, the vortex grows with the bubble until the departure happens, as shown in Figs. 11a−11c, as the area blocked by the bubble is growing. In the end, the bubble gradually deforms from elliptical to round shape due to the surface tension and rises to the upper right direction due to the buoyancy and the drag force, as shown in Figs. 11c−11e. After the departure of the bubble, the vortex is also 'blown away' by the cross flow, moving along the bottom boundary, as can be seen in Fig. 11d. The second bubble begins to grow after the departure of the first one. The disturbance to the flow field due to the first bubble still exists as shown in Fig. 11d. Nevertheless, the disturbance has less effect on the second bubble compared with the case of the quiescent liquid, since the cross flow dominates the flow pattern of the domain, see from  Then, we study the effect of cross-flow velocity and gas outflow velocity on the bubble departure diameter and releasing period. The results are presented in Figs. 12 and 13.
In Fig. 12, for brevity, the times shown are reduced by a factor of 100 from the real time steps used in the simulation. We find that the departure diameter is more sensitive to the gas outflow velocity v. The difference of between every two adjacent lines is 0.7 or 0.6. For larger u and smaller v, decreases more obviously as u increases. When , that is the Reynolds number , the reduction amplitude of is small and its change is no longer monotonous with the increase of u. In Fig. 13, the bubble releasing period decreases with the increase of u and v. For , it is going to generate the bubble with a departure diameter larger than 11, which induces large surface tension and viscous force, the process of bubble formation becomes more sensitive to the environment.

Bubble growth from the top boundary
In recent years, the drag reduction of ships has gradually become a hot research topic due to the demand of energy conservation. The air lubricating technology achieves this goal by spraying large amount of bubbles onto the bottom of a ship. The produced air 'cushion' between the water and the ship effectively reduces the frictional resistance of the ship. This process can be simplified to the bubble growing from an orifice on the top boundary in a cross-flow field. Here, we assume that the ship is still and the water flows. The computational domain has the size of . The periodic boundary condition is applied to the left and the right of the domain whilst the non-slip boundary condition is used to the bottom and the top. The orifice is set at of the top wall. The gas outflow velocity is . We define a speed ratio as . In this case, we set . The motion of the bubble is presented in Fig. 14.
At the beginning, the surface tension attaches the bubble to the orifice, the buoyancy pushes the bubble onto the wall. Meanwhile, the bubble grows downwards and rightwards due to the drag force and gas inertia force. There is no obvious vortex appears on the left side of the bubble due to the strong cross flow as mentioned above. When the bubble grows large enough, an obvious neck forms. Then the bubble departure happens, see Fig. 14c. The departed bubble then moves along the wall with a flat oval shape, due to the buoyancy and the cross flow, as shown in Figs. 14d and 14e. Then the second bubble follows the pace of the first one. During the whole process, the disturbance of the first bubble on the second bubble is negligible, due to the relatively large cross-flow velocity. We then check the influence of the speed ratio λ on the bubble. The gas releasing velocity is kept as . Relevant results are given in Fig. 15.
When the velocity of the cross flow is small, the bubble is stretched into a very long shape before departing from the orifice. Because of the buoyancy, the bubble adheres to the D b u v Fig. 12. Contour of the bubble departure diameter with cross-flow velocity and gas outflow velocity . λ λ wall and covers a large area. As the speed ratio increases, the gas separates from the orifice earlier and the space between bubbles is larger, which weakens the effect of drag reduction. Therefore, if this kind of air lubricating technology is applied to high speed ships, in order to keep a smaller speed ratio , a larger gas outflow speed is necessary, which consumes more energy. Hence, in engineering applications, one should find an optimal ship speed and gas outflow velocity, in order to realize the goal of energy saving and emission reduction. The relations between the departure diameter or releasing period and cross-flow velocity and gas outflow velocity for this case with an upper wall are investigated, as shown in Figs. 16 and 17.  Fig. 16, D b gets larger as increases while decreases with the increase of . In the cases of the speed ratio , the bottom right corner of the figure, the bubble has small size, and for cases at the top left corner of the figure, the bubble size is large when the speed ratio . Meanwhile, the same size of a bubble can be obtained by increasing both the gas outflow velocity and the cross-flow velocity linearly in these two regions. There is also an unstable region between these two regions. For the drag reduc-tion purpose, it is better to have a smaller v in order to reduce the energy consumption, and hope that the bubble size is large enough to cover the bottom of the ship for reducing the resistance. From Fig.16, there are some cases in the unstable region that for a certain u, one can find a relatively small v. Hence, the optimal velocities probably come from this relatively unstable region. From Fig. 17, the releasing period of the bubble is more sensitive to the cross-flow velocity u. For brevity, the times shown are reduced by a factor of 100 from the real time steps used in the simulation. When , that is the Reynolds number , the releasing period T of bubble is hardly affected by the gas outflow velocity . The bubble covers the hull surface better as shown in Fig. 15, under large gas outflow velocity condition. And the departure of bubble is dominated by the buoyancy and the drag force, while the gas inertia force hardly affects the departure process of the bubble, since the neck of the bubble is far from the orifice. However, the releasing period T of bubble decreases obviously with the increase of gas outflow velocity for . Because the bubble does not have a long attachment volume before leaving the orifice, the neck of bubble is near the orifice. So the gas iner-  78 CHEN Guo-qing et al. China Ocean Eng., 2020, Vol. 34, No. 1, P. 69-79