Numerical Study of Dynamics of Bubbles Using Lattice Boltzmann

Article pubs.acs.org/IECR

Numerical Study of Dynamics of Bubbles Using Lattice Boltzmann Method Sumana Ghosh and Arup K. Das* Institut Jean Le Rond d’Alembert, UMR 7190, Université Pierre et Marie Curie, Paris 6, France

Ajinkya A. Vaidya and Subhash C. Mishra Department of Mechanical Engineering, IIT Guwahati, Guwahati−781039, India

Prasanta K. Das Department of Mechanical Engineering, IIT Kharagpur, Kharagpur−721302, India ABSTRACT: The dynamics of gaseous bubbles inside a tube filled with liquid has been modeled using the lattice Boltzmann method. The diffused interface concept has been used to capture the shape of the complex interface separating two phases having high density ratio. Hydrodynamics of rising bubble inside the tube is studied in detail. Properties like densities of the phases, viscosity of the liquid, and surface tension are varied to evaluate their effects on the final shape as well as on the terminal velocity of the bubble. The volume of the bubble and the diameter of tube are also varied over a wide range to establish the effect of initial conditions on shape and the terminal velocity. Further attempts have been made to study the interaction of multiple bubbles consisting of arrays in horizontal and vertical forms. Shape distortion of one bubble due to the influence of other and merging of two bubbles due to their different uprising velocities are numerically modeled. Finally, numerical simulation is made to model Rayleigh Taylor instability which matches well with the literature. In conventional CFD models of two-phase flow, a set of Navier− Stokes equations are solved and the interface is captured by using volume of fluid (VOF) or level set methods. But interface reconstruction using VOF becomes complicated for threedimensional cases whereas level set method violates mass conservation for its application in large topological changes. Recently, lattice Boltzmann method (LBM) has emerged as an alternate and promising tool for simulation of complex twophase flow problems as compared to the conventional CFD solver for the Navier−Stokes equation. It takes care of the features of the microscale or mesoscale as well as conserves the macroscopic variables. LBM can handle multiphase systems and complex interfaces with ease and efficiency. Frisch et al.10 first proposed microscopic lattice gas cellular automaton and could reproduce the fluid flow situations. McNamara and Zanetti11 used Bhatnagar−Gross−Krook12 collision operator to make the methodology more robust and devoid of statistical noise. Chen and Doolen13 and Succi14 have reported a thorough summary of the LBM approaches and its subsequent development over the years. Shan and Chen,15 Swift et al.,16 and Gunstensen and Rothman17 have proposed several approaches to model interfacial dynamics in two-phase flow using LBM. One of the earliest multicomponent models was developed by Rothman and

1. INTRODUCTION Multiphase flows have been given due attention as they are often encountered in nature as well as in different industrial processes. Applications of multiphase flows are plenty like combustion, chemical reactions, boiling, petroleum refining, and in other heat and mass transfer processes. With the application of engineering knowledge to daily life problems, multiphase flow is becoming more relevant for different attractive fields like filling the fuel tank of a Formula 1 racing car within a few seconds, oil jet splashing within the cylinders of racing engines, etc. The simplest configuration of multiphase flow is the movement of gaseous bubbles due to buoyancy inside a liquid column. The dynamics of a bubble moving in a fluid medium has crucial ramifications in research, industrial processing, and natural processes. Due to its wide application, bubble dynamics in a liquid column has been studied by a number of researchers both experimentally1−3 and theoretically.4,5 Starting from initial experimental investigation of Zukoski6 to the recent effort of Lu and Prosperetti,7 researchers are continuously making efforts to predict the Taylor bubble shape and velocity in a stationary or moving liquid column. Bhaga and Weber8 conducted experiments to find the rise velocity and shapes of a gaseous Taylor bubble in a liquid column. The experimental investigations continued later on, and recently, Mandal et al.9 examined the shape and stability of Taylor bubbles of five different liquid pairs. They observed and reported the effect of tube diameter and tube inclination on the shape and velocity of the Taylor bubble. A lot of numerical works have also been reported to study the interfacial dynamics of the bubble in liquid ambient. © 2012 American Chemical Society

Received: Revised: Accepted: Published: 6364

July 6, 2011 March 26, 2012 April 19, 2012 April 19, 2012 dx.doi.org/10.1021/ie201445d | Ind. Eng. Chem. Res. 2012, 51, 6364−6376

Industrial & Engineering Chemistry Research

Article

Keller18 using the color method. They used two different particle distribution functions to represent two phases. Later Shan and Chen15 used the potential method to simulate two-phase flow. However, these are phenomenological models. Swift et al.16 made an improvement in the model by suggesting a free energy approach. In this model, the equilibrium distribution is defined based on thermodynamics and conservation of energy is satisfied. Lee and Lin19 developed a stabilized scheme of discrete Boltzmann equation for multiphase flows with large density ratio. However, it does not completely recover the lattice Boltzmann equation for interface to the Cahn−Hilliard equation. On the other hand, a large density ratio between the fluids can cause numerical instability. To overcome this difficulty, Inamuro et al.20 proposed a model, based on the free energy method for multiphase flows with large density ratio. In a recent study by Zheng et al.,21 another model for large density ratio of two fluids has been proposed. In their model, the lattice Boltzmann equation for interface recovers the Cahn−Hilliard equation. They have modeled bubble dynamics for a different range of associated nondimensional parameters. But still a proper description of interfacial behavior of the bubbles for different magnitudes of the physical properties is not clearly established. Moreover interactions between the uprising bubbles and their dynamic shapes have not been investigated numerically. This paper thus aims at numerically understanding the physical behavior of the bubbles and their mutual interactions inside a liquid column. The diffused interface based lattice Boltzmann model as proposed by Zheng et al.21 has been taken as the tool for the simulation. The inputs from the numerical model can be utilized to understand the dynamics of different processes associated with two phase flow. In the next section, the model has been described in detail. Extrapolation of the described model can track complex interfacial behavior in bubble column reactor and cooling of nuclear rod bundles. Simulation results for bubble dynamics under different thermo-physical properties are presented next. Ideas from the numerical observations can evolve into a better design of process equipment involving two-phase flow. Interactions of the bubbles in vertical and horizontal array are also investigated. This is a close simulation of coalescence and breakage of bubbles in the bubble column and can be described as the possible option for transition of bubbly flow to slug flow in the column. Finally, salient conclusions and future scopes of the work are highlighted in section 4.

(3)

2 ∂ϕ + ∇⃗ ·(ϕu ⃗) = θM∇⃗ μϕ ∂t

(4)

where P is the pressure tensor, θM is the mobility, μϕ is chemical potential of the fluid, and F⃗b is the body force. The interface is considered to be diffused within a narrow zone and it has been modeled using convective Cahn−Hillard equation (eq 4). In order to explain eq 4, in lattice structure, a modified lattice Boltzmann equation is adopted21 using probability distribution function gi as gi(x ⃗ + ci⃗δt , t + Δt ) − gi(x ⃗ , t ) = (1 − q)[gi(x ⃗ + ci⃗δt , t ) − gi(x ⃗ , t )] + Ωgi]

(5)

where, the last term of eq 5 is the collision operator and can be written as gieq (x ⃗ , t ) − gi(x ⃗ , t )

Ωgi =

τϕ

(6)

In eq 6 τ ϕ is dimensionless relaxation time, ci⃗ is lattice velocity, and q is a constant coefficient to determine the implicit parameter. The value of τ ϕ depends on the physical properties of the fluid pairs and accordingly the value of implicit parameter is chosen. The value of q is related with τ ϕ as 21

q=

1 τϕ + 0.5

(7)

After calculating the directional gi for the lattice, the macroscopic order parameter ϕ is evaluated as ϕ=

∑ gi

(8)

i

Equation 5 recovers the Cahn−Hillard equation with secondorder accuracy as shown in the following equation: 2

∂tϕ + ∇⃗ ·(ϕu ⃗) − θM∇⃗ μϕ + O(δ 2) = 0

(9)

In this equation, mobility θM is defined using lattice properties as

2. MODEL DESCRIPTION A basic model to track the two-phase flow is developed based on lattice Boltzmann methodology. The diffused interface concept is considered to track the interfacial behavior. As the phenomena is associated with flow of two fluids or phases with different densities, their respective densities are considered as ρL and ρH, respectively. To decouple the interfacial dynamics from the bulk motion, modified momentum equations are considered. In these equations, two derived properties n and ϕ are used which can be defined as22 ρ + ρH ρ − ρL ,ϕ= H n= L (1) 2 2 Using the above derived quantities, the flow can be described by mass and momentum conservation equations and an interface evolution equation as22 ∂n + ∇⃗ ·(nu ⃗) = 0 ∂t

2 ∂(nu ⃗) + ∇⃗ ·(nuu⃗ ⃗) = −∇⃗ ·P + μ∇⃗ u ⃗ + Fb⃗ ∂t

⎛ 1⎞ θM = q⎜τϕq − ⎟Γ ⎝ 2⎠

(10)

where, Γ is the mobility controlling parameter for the fluid pair. In the present model, D2Q5 lattice structure has been considered during the propagation of the interfacial information. In D2Q5 structure, equilibrium distribution function of gi can be calculated as21 gieq = Ai + Bi ϕ + Ciϕ ci⃗ ·u ⃗

(11)

The coefficients in eq 11 are taken as B1 = 1, Bi = 0 (i ≠ 1) Ci =

1 2q

Ai = −2Γμϕ

(2) 6365

(12)

dx.doi.org/10.1021/ie201445d | Ind. Eng. Chem. Res. 2012, 51, 6364−6376


Article

As the interface is related to surface tension across it, the surface tension force can be calculated as → ⎯ Fs = −∇⃗ ·P = −ϕ∇⃗μ − ∇⃗p ϕ

(13)

o

where po is the static pressure given by po = ncs2, cs being the speed of sound. Due to incorporation of diffused interface concept and generated surface forces, modified mass and momentum conservation equations can be written as

∂n + ∇⃗ ·(nu ⃗) = 0 ∂t

(14)

∂(nu ⃗) + ∇⃗ ·(nuu⃗ ⃗) ∂t 2 = −ϕ∇⃗μϕ − ∇⃗po + μϕ ∇⃗ϕ + Fb⃗ + μ∇⃗ u ⃗

(15)

The chemical potential μϕ described in eq 15 is computed by using following relationship:21 μϕ = A(4ϕ3 − 4ϕ*2 ϕ) − κ ∇2 ϕ

(16)

where A is the amplitude parameter used to control the interaction of energy between the two phases and κ is the curvature. Due to the incorporation of diffused interface, variation of ϕ near the interface is calculated as ⎛ 2ζ ⎞ ϕ = ϕ* tanh⎜ ⎟ ⎝W ⎠

Figure 1. Schematic diagram of the computational domain.

(17)

where ζ is the coordinate which is perpendicular to the interface and W is the thickness of the interface layer. This can be expressed as21

W=

Lattice Boltzmann implementation of eqs 14 and 15 can be given as fi (x ⃗ + ci̅ δt , t + δt ) = fi (x ⃗ , t ) + Ω fi

2κ /A ϕ*

(18)

Here, the expected order parameter is ρ − ρL ϕ* = H 2

with a collision parameter as Ω fi =

(19)

(∇⃗ϕ2) n p = A(3ϕ − 2ϕ* ϕ − ϕ* ) − κϕ∇ ϕ + + 2 3 2 2

4

2

⎛ ⎞ 3 2 9 ⎜ f ieq = wA u + uαuβciαciβ ⎟ i i + wn i ⎝3ciαuα − ⎠ 2 2

And surface tension of a fluid pair is related with the derived properties as

⎛

⎞

∫ κ ⎜⎝ ∂∂ϕζ ⎟⎠

dζ

⎛ 2ζ ⎞ ⎛ 2ζ ⎞ 3σ tanh3⎜ ⎟ sech2⎜ ⎟ 2 ⎝ ⎠ ⎝W ⎠ W W

(24)

(25)

where the coefficients are as follows:

2

A1 =

(21)

15(ϕμϕ + n/3) 9 n− 4 4

Ai = 3(ϕμϕ + n/3)

Simplifying eqs 16 and 21, one can write μϕ ∇ζ ϕ =

τn

⎛ 1 ⎞ wi + ⎜1 − ⎟ 2τn ⎠ cs 2 ⎝

The equilibrium distribution function in eq 24 can be written as21

(20)

σ=

f ieq (x ⃗ , t ) − fi (x ⃗ , t )

⎡ c u⃗ ⃗ ⎤ × ⎢( c ⃗ − u ⃗) + 2 c ⃗ ⎥(μϕ ∇ϕ + Fb⃗ )δt cs ⎦ ⎣

The pressure tensor is calculated as 4

(23)

w1 =

(22)

4 1 1 , wi = 2,3,4,5 = , wi = 6,7,8,9 = 9 9 36

(26)

τn is the relaxation parameter and can be related with viscosity as21

Thus, the potential form of surface tension related term is independent of the average density and density difference. It is obvious from eq 22 that μϕ∇ζϕ is related with the surface tension coefficient and the width of interface layer due to the incorporation of diffused interface concept. For discretization of mass and momentum conservation, eq D2Q9 lattice structure is used as proposed by Zheng et al.21

μ = cs 2(τn − 0.5)n

(27)

3. RESULTS AND DISCUSSION Motion of a gaseous bubble under the action of buoyancy and the corresponding interfacial configurations are studied using 6366



Article

Figure 2. Numerical validation of the developed model along with the work of Zheng et al:21 (a) initial condition, (b) present simulation at 0.4 s, (c) Zheng et al.21 at 0.4 s, (d) present simulation at 0.6 s, (e) Zheng et al.21 at 0.6 s.

the described model. A spherical bubble is initially kept under the cylindrical liquid pool and allowed to move up due to buoyancy force. During the upward movement of the bubble, it will change its shape inside the domain. The movement of the bubble and its shape has been validated with published numerical data of Zheng et al.21 After the validation of the developed numerical code, the effect of different thermophysical parameters of the liquid−gas pair on bubble shape and its terminal velocity is investigated. Study has also been made to simulate mutual interactions of the bubbles during its presence inside the liquid column. Multiple bubbles or array of bubbles are taken for the simulation of such cases. 3.1. Numerical Validation. Developed axisymmetric numerical model is first validated with the published numerical results of Zheng et al.21 A basic schematic diagram of the air bubble in a water column is shown in Figure 1. Air and water have been taken as the working fluids for the simulation. The domain for validation consists of 120 × 240 lattices. A spherical bubble with radius of 20 lattices is located at the bottom (60, 60) of the liquid column. The densities, surface tension, and mobility are set as ρH = 1000.0, ρL = 1.0, σ = 2.0, Γ = 100.0

Table 1. Comparison of Uprising Bubble Velocity Obtained from Present Numerical Simulation and the Work of Bozzano and Dente23 bubble diameter (mm)

VT (present) (m/s)

VT (Bozzano and Dente (2001)) (m/s)

11 12 13

0.195 0.288 0.321

0.2 0.295 0.33

complexity of the interface configurations as has been observed in other established methodologies. The dimensionless numbers associated with this problem can be defined as the Eotvos, Morton, Reynolds numbers. The functional forms of these numbers are as follows: Eo =

g (ρH − ρL )d 2 σ

,M=

g (ρH − ρL )μH 4 2 3

ρH σ

, Re =

ρH VTd μH (29)

For the present input conditions of the simulation, the dimensionless numbers take the following values: Eo = 39.96, M = 1.52425, Re = 5.984. The shape of the bubble and its location are tracked at different instants for comparison with Zheng et al.21 Figures 2a−e show the comparison of the present numerical simulation (Figures 2b and d) and published results (Figures 2c and e) of

(28)

The proposed model is computationally simple, and for a single Pentium IV (3 GB RAM) processor, the average CPU time is around 15 min. Numerical efforts increase with the 6367



Article

Figure 3. Bubble shape and position for different densities of the discrete phase keeping water as the continuous phase: (a) air−water, (b) liquefied LPG−water, (c) kerosene−water.

Figure 6. Variation of bubble velocity for different density of the continuous phase.

Figure 4. Variation of bubble velocity for different density of the discrete phase.

Zheng et al.21 A good match can be seen from the contours of the bubbles at different time instants. Velocity of an uprising bubble is another important parameter which needs to be validated with the literature. As the bubble velocity is highly dependent on its volume, the

numerical velocity of air bubbles in a water column is calculated for three different volumes and compared with observations of Bozzano and Dente.23 11, 12, and 13 mm diameter spherical air bubbles are placed in a 20 mm diameter tube which is filled with water. A lattice resolution of 100 × 300 is assigned for the

Figure 5. Bubble shape and position for different densities of the continuous phase keeping liquefied LPG as discontinuous phase: (a) liquefied LPG−water, (b) liquefied LPG−chloroform, (c) liquefied LPG−citric acid. 6368



Article

Figure 7. Bubble shape and position for different surface tension of the fluid pair: (a) 0.96, (b) 0.32, (c) 0.045 N/m.

Figure 8. Variation of bubble velocity for different values of surface tension.

Figure 10. Variation of bubble velocity for different values of viscosity or τϕ.

simulation. Simulations are performed, and the terminal velocities are calculated for each input parameter. The velocities are compared with Bozzano and Dente23 after the bubble has traversed 0.4 s up in the liquid column. Periodic boundary conditions are employed at top and bottom walls in order to calculate stable terminal velocities of the bubble. The

numerical values of the velocities for different conditions are reported in Table 1. Velocities of the bubble for identical conditions as reported by Bozzano and Dente23 which are also tabulated along with the present numerical simulation. A good match in the velocity prediction can be observed from the table.

Figure 9. Bubble shape and position for different values of τϕ: (a) τϕ = 0.65, (b) τϕ = 0.75, (c) τϕ = 0.85. 6369



Article

3.2. Effect of Various Thermophysical Properties. Next, the effect of densities of liquid and gas phases, viscosity of the liquid phase, and surface tension on the shape and velocities of the bubble are studied in detail. Variation of all these properties keeping the others properties the same will establish the influence of different thermophysical properties on the dynamics of the bubble. 3.2.1. Effect of Phase Densities. In this section, simulations have been made to identify the effect of discrete phase density on the shape and terminal velocity of the bubble. Water (1000.0 kg/m3) is kept fixed as the continuous phase for air− water, liquefied LPG−water, and kerosene−water combinations in the simulations. The corresponding densities of the discrete phases are 1.2, 560, and 787 kg/m3. All other properties are kept same to only see the effect of discrete phase density on bubble dynamics. Surface tension is taken as 0.05 N/m for all the studies. The diameter and length of the tube are taken as 12 and 30 mm, respectively. A mesh size of 120 × 300 is employed to study the effect of discrete phase density. After 0.4 s from the initial placement of the spherical bubble, we study the shapes and the positions of the bubble of different fluids inside the liquid water column. Figures 3a−c show the shapes and positions of the bubble for (a) air−water, (b) liquefied LPG−water, and (c) kerosene−water combination, respectively. It can be seen from the figures that air bubble in water atmosphere (Figure 3a) moves faster as compared to the other two combinations (Figure 3b and c) due to its higher density difference between the phases. After the same time interval, the movement of the kerosene bubble is smaller compared to the liquefied LPG and air bubble. The shapes of the bubbles are also different due to the variable effect of buoyancy force. The kerosene bubble is more spherical whereas the liquefied LPG bubble is more distorted as compared to the air bubble in water column. Figure 4 shows a variation of bubble velocity with respect to different density of the discrete phase, ρ. The density of discrete phase has been nondimensionalized by the density of liquefied LPG as ρ/ρLPG. It can be clearly observed from the figure that with decrease in density difference between the continuous and discrete phases (ρ − ρwater), the terminal velocity of bubble decreases drastically. Further efforts have also been made to investigate the effect of continuous phase density on the dynamics of the growing bubble. For such investigation, liquefied LPG is kept as the discrete phase throughout and the continuous phase is varied. Water (1000 kg/m3), chloroform (1460 kg/m3), and citric acid (1660 kg/m3) are taken as continuous phases keeping all other properties the same as the previous situation. Figures 5a−c show the bubble contours for different fluid combinations after 0.48 s from the initial placement of the spherical bubble. From the figure, it is again clear that the bubble having high density difference (liquefied LPG−citric acid (Figure 5c) moves faster compared to the other pairs (liquefied LPG−water in Figure 5b and liquefied LPG−chloroform in Figure 5a). The shape of the bubble is also dependent on the density of the continuous phase. It is clear from the figure that the bubble remains spherical when the density difference between the two phases is small. For high density difference, the bubble turns into a shape similar to a Taylor bubble. Effect of continuous phase density on terminal velocity of the bubble is shown in Figure 6. The density of the continuous phase is nondimensionalized using the density of chloroform. It is observed from the figure that with the increase of

Figure 11. Bubble shape and position for different initial bubble size: db = (a) 2, (b) 4, (c) 6, (d) 8 mm.

Figure 12. Variation of bubble velocity for different values of bubble diameter.

continuous phase density higher terminal velocity of the bubble is achieved. 3.2.2. Effect of Surface Tension. To understand the effect of surface tension on the shape of the bubble, keeping all other properties same, simulations have been made for different values of surface tension. As the phenomenon is buoyancy dominated, to understand the effect of surface tension, a bigger 6370



Article

Figure 13. Bubble shape and position for different tube diameter: D = (a) 12, (b) 15, (c) 18 mm.

Figure 14. Variation of bubble velocity for different values of tube diameter.

the tube as the surface tension decreases. Moreover, from the location of bubble front, it is also clear that bubble moves faster in fluid pairs having low values of the surface tensions. The velocities of the bubbles are also found out for three different values of surface tensions (0.96, 0.32, and 0.045) as mentioned earlier. Results are shown in Figure 8. It can be observed from the figure that bubble velocity decreases with the increase of the values in surface tension. With low surface tension coefficient, the bubble not only moves up fast but it also experiences a dominant change in shape.

bubble has been simulated. Simulations have been started from a cylindrical bubble inside the tube. Water and kerosene are selected for continuous and discrete phases, respectively. Simulations have been done for three different values (0.96, 0.32, and 0.045 N/m) of surface tension for the same domain and lattice structures as considered in the previous cases. Figure 7 shows the shape of the bubbles for three different surface tensions after 0.48 s from its initial placement in the column. It can be observed from the results that for higher value of the surface tension (0.96 N/m), the bubble remains cylindrical. The bubble becomes slender and thin inside 6371



Article

Figure 15. Mutual influence of bubbles placed in a horizontal array: t = (a) 0, (b) 0.2, (c) 0.4, (d) 0.6 s.

3.2.3. Effect of Viscosity. In the present model a common viscosity is shared by both the fluids. In this model, it is very difficult to assess the effect of viscosities separately. As the viscosity μ is directly related with the value of τn (eq 27), simulations have been made with different values of τn to show the influence of viscosity on bubble dynamics. Three different values of τn (0.65, 0.75, 0.85) are considered for the simulation of air bubble in water column. A mesh size of 120 × 240 is used for the tube diameter being 12 mm and bubble diameter of 5 mm. The bubble is kept at the center of the tube below the liquid column with length of the tube being 30 mm. Phase contours for three different values of τn are shown in Figure 9. It can be seen from the figure that the interfacial contours change for variation in the liquid viscosity. For lower value of τn, the rear end of the bubble becomes more distorted but the crater at rear end does not penetrate much inside the bubble for higher value of τn. From the figure, it is also clear that the liquid viscosity has minute effect on the terminal velocity of the uprising bubble. As the phenomenon is dominated by buoyancy force the effect of viscosity is not so significant. Effect of viscosity over the terminal velocity of the bubble is also studied for three different value of τn. Figure 10 reports the variation of terminal velocity for variation of τn over a range. It has been observed that the change in velocity is not drastic as

the phenomenon is primarily dominated by buoyancy. But still the effect of different viscosity can be tracked from the results. 3.3. Effect of Bubble Diameter. Simulations have also been made to identify the effect of initial bubble diameter on its shape and uprising velocity. The tube diameter is kept constant as 12 mm while the diameter of the air bubble is varied as 2, 4, 5, and 6 mm respectively. Liquid water is taken as the continuous phase. The domain is discretized into 120 × 300 lattices. The bubbles are allowed to move up under the action of buoyancy force from their initial conditions. After same time interval, their shapes are observed and velocities are calculated. The bubble positions are noted after a time of 0.4 s from its initial placement in the column. The results are shown in Figure 11. It can be observed from the figures that the bubble with the large diameter moves faster as the buoyancy force is proportional to the secondary phase volume. The shape of the bubble is also influenced due to its initial size. The bubble with large volume distorts more from its initial condition as compared to the bubble having smaller diameter. A gas bubble released in a liquid medium will rise and it will move to a point where its buoyancy force, is balanced by the drag force. In turn, the bubble velocity is influenced by its initial diameter. With increase in diameter, the buoyancy force and 6372



Article

drag force increase. But the increase in buoyancy force is more as compared to drag force, and hence, the velocity of bubble increases with diameter of the bubble. This can be clearly seen from Figure 12 where bubble velocity is reported for different bubble diameters. 3.4. Effect of Tube Diameter. Next, we analyze the effect of the tube diameter on the dynamics of air bubbles in water column. For this, the tube diameter is varied as 12, 15, and 18 mm, while the length is kept constant as 24 mm. The bubble diameter is considered to be 10 mm. Lattice resolutions of 120 × 240, 150 × 240, and 180 × 240 are generated to replicate three different tube diameters. The snapshots of the numerical simulations after a time interval of 0.21 s are shown in Figure 13. It can be observed from the figures that tube diameter has a dominant effect on the bubble shape and its mobility. Bubble becomes spherical as the tube diameter increases due to fading effect of no slip boundary condition at the wall. On the other hand for a low tube diameter, the bubble takes a long slender shape due to dominant effect of wall boundary condition. The velocity of the bubble is also influenced by the change of tube diameter as the velocity pattern is different for the surrounding liquid domain in different tubes. It can be observed that the bubble at higher tube diameter climbs more distance as compared to the bubble at smaller tube diameter. This can be also made clear from the velocity plot of the bubbles for different tube diameters. For four different tube diameters, the velocities of the bubbles are depicted in Figure 14. For the 12 mm tube diameter, the terminal velocity is 0.23 m/s, while for the 15 mm tube diameter, the terminal velocity is 0.269 m/s. The corresponding velocity of the bubble inside a tube having 18 mm diameter is found to be 0.295 m/s. It can be concluded that with the increase in tube diameter, the bubble gets flattened and also its velocity increases. With increase in tube diameter the drag force due to wall effect decreases considerably, and hence, the net terminal velocity increases. The effect of wall drag nullifies quickly as the tube diameter increases further. In the case of 25 mm tube diameter, the bubble velocity remains almost similar to the observed velocity of a bubble inside an 18 mm diameter tube. This clearly represents the fading nature of the wall interaction with the bubble inside a large tube. 3.5. Horizontal Array of Bubbles. Having studied the effect of various properties on the shape and velocity, next we study the interactions of the bubbles in a liquid column and understand the underlying physics. Interaction of bubbles and their mutual dynamics are very common in nature. It can lead to typical distortion of the interfaces or merging to two bubbles. First, we study the case of a horizontal array of bubbles. Two air bubbles of different diameters and proper spacing are placed side by side in a water column. 160 × 320 lattices are used to simulate the bubble dynamics in a 2D rectangular domain. The channel width is considered to be 16 mm, and the length is taken as 30 mm. The diameter of the big bubble is taken as 4 mm whereas the small one is considered to be 2.4 mm. It has been observed earlier that bigger bubbles move fast as compared to their smaller counterpart. The objective of placing small and big bubbles together in the liquid domain is basically to study the influence of a big bubble on a small bubble and vice versa. Snapshots are taken at a time interval of 0.2 s. In Figure 15a−d, four snapshots are reported at 0, 0.2, 0.4, and 0.6 s, respectively, for phase and velocity contours in the domain. It can be observed from the figure that initially the

Figure 16. Mutual influence of bubbles placed in a vertical array: t = (a) 0, (b) 0.1 (c) 0.25, (d) 0.4 s.

Figure 17. Velocity vectors of the liquid domain (a) before and (b) after the merging of bubbles.

bubbles are not influenced by each other but as time progresses smaller bubble is more influenced by the bigger bubble. It has been observed that bubbles moves up with high velocity which generates flow in the surrounding liquid. The influence on the 6373



Article

Figure 18. Mutual influence of bubbles placed in matrix: t = (a) 0, (b) 0.15, (c) 0.4, (d) 0.5 s.

surrounding increases with the bubble size and interactions with each other. It is seen from Figure 15 that the shape of the smaller bubble deviates from its original spherical shape. It becomes flat and disklike due to the influence of the differential velocity field of the surrounding liquid. A bigger bubble is also influenced by the smaller bubble in the same manner, but its magnitude of influence is smaller as compared to its counterpart. As a result, the shape of the bigger bubble also distorts but not much from its original spherical shape. As the velocity of the bigger bubble is higher (velocity contours) as compared to the smaller bubble, their mutual interaction fades away after some period and they move individually by following their own dynamics. 3.6. Vertical Array of Bubbles. Next we consider the case in which two bubbles are placed one over the other. Earlier (Figure 11) it has been established that the bigger bubble has higher velocity as compared to the smaller bubble for air−water combination due to variation in the buoyancy effect. Keeping this in mind, a bigger bubble (2 mm diameter) is placed under a smaller bubble (1 mm diameter) to show their mutual dynamics as these bubbles will come closer due to a change in velocity variation. A 120 × 300 lattice structure is used for the simulation. The tube diameter is kept at 12 mm. Snapshots of phase and velocity contours taken at 0, 0.1, 0.25, and 0.4 s are depicted in Figure 16a−d, respectively. It can be observed from the figure that with time, the bigger bubble approaches toward the smaller bubble and eventually merges into a single bubble. In this case the influences of the bubbles on the surrounding liquid interact with each other and result in drainage of liquid film between the bubbles. Later on, the single bubble follows its own dynamics based on its own volume and neighboring conditions. Before merging of the bubbles, it can

be also observed that the individual shapes of the bubbles are different due to their variation in size. To probe in the liquid flow pattern around the merging bubbles, the velocity vectors of the surrounding liquid are depicted in Figure 17. From the figure, it is clear that drainage of liquid film advances the process of merging. After merging of the bubbles, a big trailing vortex can also be observed from the velocity pattern (Figure 17b). 3.7. Simultaneous Horizontal and Vertical Arrays of Bubbles. Having studied the dynamics of horizontal and vertical arrays of bubbles, next we take up the case in which we have both horizontal and vertical arrays together. It will be a close simulation of bubble column. In a rectangular water column of width 24 mm, we place 1 mm diameter air bubbles over 2 mm horizontally arranged bubbles (Figure 18 a). The bubbles move in a 3 × 2 matrix. Both the horizontal (4 mm) and vertical (2 mm) spacing between the bubbles are kept uniform. But for the present simulation, the vertical spacing is considered to be much lower compared to the horizontal spacing in order to only see the effect of bubble merging in columns. Positions and shapes and velocity vectors of the 3 × 2 matrix of bubbles are shown at 0.15, 0.4, and 0.5 s in Figure 18b−d. It can be seen from the figure that the bubbles merge with each other vertically and advance forward obeying their own dynamics. The bubble near the wall moves little faster as compared to the bubble at the center of the conduit. 3.8. Rayleigh Taylor Instability. In order to establish the predictability of the developed model for more complex interfaces, simulations have been made for classical Rayleigh Taylor instability. Rudman,24 Popinet and Zaleski,25 and Sussman and Puckett26 used a similar type of case study to establish their methodologies. A 1 m × 4 m rectangular domain 6374



Article

is dicretized using 64 × 256 lattices. The fluid densities are taken as 1.225 and 0.1694 kg/m3. Viscosities of both the fluids are considered as 0.00313 kg/m S. Figure 19 shows the

■

shows that present model is capable of tracking evolution, breaking, and making of interface separating two fluids having high density difference.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Tel.: +33-130854869. Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Tung, K. W.; Parlange, J. Y. Note on the motion of long bubbles in closed tube influence of surface tension. Acta Mechanica 1976, 24, 313. (2) Kataoka, Y.; Suzuki, H.; Murase, M. Drift-flux parameters for upward gas flow in stagnant liquid. J. Nucl. Sci. Technol. 1987, 24, 580. (3) Polonsky, S.; Shemer, L.; Barnea, D. Relation between the Taylor bubble motion and the velocity field ahead of it. Int. J. Multiphase Flow 1999, 25, 957. (4) Bretherton, F. P. The motion of long bubbles in tubes. J. Fluid Mech. 1961, 10, 166. (5) Tudose, E. T.; Kawaji, M. Experimental investigation of Taylor bubble acceleration mechanism in slug flow. Chem. Eng. Sci. 1999, 54, 5671. (6) Zukoski, E. E. Influence of viscosity, surface tension and inclination angle on motion of long bubbles in closed tubes. J. Fluid Mech. 1966, 25, 821. (7) Lu, X.; Prosperetti, A. Axial stability of Taylor bubbles. J. Fluid Mech. 2006, 568, 173. (8) Bhaga, D.; Weber, M. E. Bubbles in viscous liquids: shapes, wakes and velocities. J. Fluid Mech. 1981, 105, 61. (9) Mandal, T. K.; Das, G.; Das, P. K. Motion of Taylor bubbles and Taylor drops in liquid−liquid systems. Ind. Eng. Chem. Res. 2008, 47, 7048. (10) Frisch, U.; Hasslacher, B.; Pomeau, Y. Lattice-gas automata for the Navier-Stokes equation. Phys. Rev. Lett. 1986, 56, 1505. (11) McNamara, G. R.; Zanetti, G. Use of the Boltzmann equation to simulate lattice-gas automata. Phys. Rev. Lett. 1988, 61, 2332. (12) Bhatnagar, P. L.; Gross, E. P.; Krook, M. A model for collision processes in gases. Phys. Rev. 1954, 94, 511. (13) Chen, S.; Doolen, G. Lattice Boltzmann method for fluid flows. Annu. Rev. Fluid Mech. 1998, 30, 329. (14) Succi, S. The lattice Boltzmann equation for fluid dynamics and beyond; Clarendon Press: Oxford, 2001. (15) Shan, X.; Chen, H. Lattice Boltzmann model for simulating flows with multiple phases and components. Phys. Rev. 1993, 47, 1815. (16) Swift, M. R.; Osborn, W. R.; Yeomans, J. M. Lattice Boltzmann simulation of nonideal fluids. Phys. Rev. Lett. 1995, 75, 830. (17) Gunstensen, A. K.; Rothman, D. H. Microscopic modeling of immiscible flows in three dimensions by a lattice Boltzmann method. Europhys. Lett. 1998, 18, 157. (18) Rothman, D. H.; Keller, J. M. Immiscible cellular automaton fluid. J. Stat. Phys. 1988, 52, 1119. (19) Lee, T.; Lin, C. L. A stable discretization of the lattice Boltzmann equation for simulation of incompressible two-phase flows at high density ratio. J. Comput. Phys. 2005, 206, 16. (20) Inamuro, T.; Yoshino, M.; Inoue, H.; Mizuno, R.; Ogino, F. A lattice Boltzmann method for a binary miscible fluid mixture and its application to a heat transfer problem. J. Comput. Phys. 2002, 179, 201. (21) Zheng, H. W.; Shu, C.; Chew, Y. T. A lattice Boltzmann for multiphase flows with large density ratio. J. Comput. Phys. 2006, 218, 353. (22) Takada, N.; Misawa, M.; Tomiyama, A.; Hosokawa, S. Simulation of bubble motion under gravity by Lattice Boltzmann Method. J. Nucl. Sci. Technol. 2001, 38, 330. (23) Bozzano, G.; Dente, M. Shape and terminal velocity of single bubble motion: a novel approach. J. Comput. Chem. Eng. 2001, 25, 571.

Figure 19. Rayleigh Taylor instability on a 64 × 256 lattice structure along with the numerical results of Popinet and Zaleski.25

evolution of the interface predicted by the numerical simulation at times 0.7, 0.8, and 0.9 s. For comparison, predictions of Popinet and Zaleski25 are also depicted in the figures which show a good match with the complex numerical interface configurations. This shows the potential of the developed methodology for tracking the complex interfacial nature.

4. CONCLUSIONS The lattice Boltzmann method has been used to study the dynamics of a gaseous bubble in a liquid column having high density ratio (≈1000). Results of numerical simulation reveal the following points: (i) Simulations for air−water, liquefied LPG−water, kerosene−water, liquefied LPG−chloroform, and liquefied LPG−citric acid show that the bubble velocity increases with increase in density difference between the phases. (ii) It has been observed from numerical simulation that the velocity of bubble decreases (≈20%) significantly with the increase in surface tension of the fluid pair. On the other hand, a marginal increase (≈3.4%) in bubble velocity is recorded for increase in common fluid viscosity. (iii) The effect of bubble diameter and tube diameter on bubble velocity and its shape has also been investigated in detail. (iv) The developed model has also been tested for mutual interactions of bubbles during their motion in the liquid column. Bubbles placed in horizontal and vertical arrays are simulated using the model. It has been observed that the motion of one bubble is biased by the dynamics of the neighboring bubble. (v) In case of a vertical array, situations (1 mm diameter placed above 2 mm diameter) have been simulated where two bubbles are merged to form a single bubble during their motion. (vi) Finally, simulations have been made successfully for Rayleigh Taylor instability (density ratio 7.3) which 6375



Article

(24) Rudman, M. Volume tracking methods for interfacial flow calculations. Int. J. Numer. Methods Fluids 1997, 24, 671. (25) Popinet, S.; Zaleski, S. A front-tracking algorithm for accurate representation of surface tension. Int. J. Numer. Methods Fluids 1999, 30, 775. (26) Sussman, M.; Puckett, E. G. A coupled level set and volume of fluid method for computing 3D and axisymmetric incompressible twophase flows. J. Comput. Phys. 2000, 162, 301.

6376


Numerical Study of Dynamics of Bubbles Using Lattice Boltzmann

Recommend Documents