Direct Numerical Simulation of Bubble Dynamics Using Phase-Field

Article pubs.acs.org/IECR

Direct Numerical Simulation of Bubble Dynamics Using Phase-Field Model and Lattice Boltzmann Method Shuli Shu†,‡ and Ning Yang†,* †

State Key Laboratory of Multiphase Complex Systems, Institute of Process Engineering, Chinese Academy of Sciences, P.O. Box 353, Beijing 100190, People’s Republic of China ‡ University of Chinese Academy of Sciences, Beijing 100049, China ABSTRACT: Knowledge of bubble dynamics provides a microscale or mesoscale basis to further understand the macroscale complex behavior of bubbly flow or extract the constitutive correlations for continuum-based multifluid computational fluid dynamics (CFD) models. This study attempts to investigate the bubble dynamics from a systematic and multiscale perspective, that is, progressively probing the behavior of a single bubble, a bubble pair, and a bubble swarm with a Multiple-Relaxation-Time (MRT) lattice Boltzmann method (LBM). The MRT-LBM method is used to solve the phase-field model and can accurately capture the interface evolution under different flow conditions. This study also incorporates the force term in a different way, and the model is validated by classical numerical tests for a static bubble and capillary wave. Numerical tests show that the model is able to simulate the gas−liquid flows with a large kinematic viscosity ratio (up to 1000) and the bubbles of lower Morton number and higher Reynolds number. The simulation results of the terminal rise velocity of air bubbles are comparable with the classical experimental correlations. Interactions of two bubbles vertically or horizontally aligned are simulated, reproducing the trailing bubble catching up with the leading one for a vertically aligned bubble pair and the repeatedly bouncing and attracting phenomenon for a horizontally aligned bubble pair. The behavior of bubble pairs in a bubble swarm is analyzed with the model to understand the different bubble−bubble interactions and their influence on micromixing.

■

INTRODUCTION Bubble columns are frequently utilized in chemical and process industries, because of the excellent characteristics of heat and mass transfer. While the scaleup and optimization requires an extensive understanding of mixing and fluid dynamics of gas− liquid flows, the bubble dynamics provides a microscale or mesoscale basis to further understand the macroscale complex behavior and derive the constitutive laws of continuum computational fluid dynamics (CFD) models for simulating the macroscale behavior, such as gas holdup distribution, bubble size distribution, and liquid circulation. Theoretical and experimental analysis offers two classical approaches to study bubble dynamics. For example, Clift et al.1 generalized a regime map of bubble shapes with respect to three dimensionless parameters, i.e., the Eötvös number (Eo), the Morton number (Mo), and the Reynolds number (Re) for bubbles in unhindered gravitational motion through liquids, and proposed several correlations for terminal velocity and drag coefficient of isolated bubbles. Fan and Tsuchiya2 then presented wider applicable correlations of rise velocity and drag coefficients for isolated bubbles and investigated the bubble wake dynamics in liquid and liquid−solid suspensions. However, theoretical analysis is limited to low Re number or unbounded flow conditions, and correlations from experimental measurement are difficult, if not impossible, to reach the microscale and mesoscale flow details around individual bubbles. Hence, so far, no complete knowledge has been provided by experiments or theoretical analysis, because of the complexity of gas−liquid flow in chemical reactors. It becomes essential to study the details of the bubble dynamics at microscale and mesoscale levels under more-complicated cases, using direct numerical simulation (DNS). © 2013 American Chemical Society

DNS of bubble dynamics has attracted more and more attention with the advance of computer science and numerical algorithms of interface tracking. Different methods have been developed for gas−liquid two-phase flow in the past. These numerical methods can fall into two classes:3 interface tracking and interface capturing. While the former mainly refers to the front tracking (FT) method,4 the latter involves the volume of fluid (VOF) method5 and the level set (LS) method.6 For the FT method, the interface within one grid can be accurately captured; however, mass is not always conserved, and it is also difficult to deal with bubble breakage or coalescence or extend the twodimensional (2D) model to a three-dimesnional (3D) simulation. The VOF method can maintain mass conservation but cannot guarantee numerical accuracy, because of the difficulty in computing local curvature.7 The LS method is easy to implement, but it cannot retain the mass conservation. One common disadvantage of the above-mentioned methods is the need to solve the pressure Poisson equation, which is time-consuming and not suitable for parallel computation. In recent years, the lattice Boltzmann method (LBM) has been receiving more and more attention, because of its capability to model interfaces, and it is more suitable for massive computation.8 LBM bridges the macroscale flow and microscale molecular motion by solving the mesoscopic kinetic Special Issue: Multiscale Structures and Systems in Process Engineering Received: Revised: Accepted: Published: 11391

December 16, 2012 February 24, 2013 February 28, 2013 February 28, 2013 dx.doi.org/10.1021/ie303486y | Ind. Eng. Chem. Res. 2013, 52, 11391−11403

Industrial & Engineering Chemistry Research

Article

viscosity, because of the numerical instability. It is still challenging to simulate bubbles at high Re numbers, low Mo numbers, and two phases with large kinematic viscosity ratio by LBM. In this paper, a MRT-LBM model is used to solve the PF model and simulate bubbles at high Re number (up to 1000), low Mo number (air−water system), and high viscosity ratio for two phases. The large kinematic viscosity ratio is up to 1:1000 for the case of static bubbles. The motivation lies in the following aspects. First, the LBM model itself, as a mesoscopic approach, is suitable and promising to probe the microscale and mesoscale behavior of complex gas−liquid flow. This work can serve as a basis to further incorporate the mass and heat transfer in the next step. Second, force is incorporated in the multiphase LBM model in a different way with previous PF models as far as we know. Third, this work investigates the bubble dynamics from a systematic and multiscale perspective, that is, simulating the behavior of a single bubble, a bubble pair, and a bubble swarm, respectively. Bubble dynamics such as the bubble shape and terminal rise velocity of single bubble, and the interactions of two bubbles and the bulk motion of multiple bubbles, are investigated. The model is validated by some classical numerical tests of the order parameter distribution, the pressure difference between static bubbles of different diameter and its surrounding liquid of different surface tension with kinematic viscosity ratio of two phases up to 1000, and the capillary wave. Then, 2D simulations are carried out to capture the variation of bubble shapes under different flow conditions and portray the behavior of a single bubble rising under gravity. The terminal velocity of a single rise bubble is compared with experimental correlations. Two-dimensional (2D) simulations for vertically or horizontally aligned bubble pairs and a 3D simulation for a bubble swarm are carried out to investigate the bubble−bubble interactions under different flow conditions.

equation. Numerically, LBM is an explicit method without timeconsuming iterations and is easy to implement. LBM can simulate multiphase flow with complex interface phenomenon without special treatments for interfaces. There are three main multiphase LBM models, i.e., the color model,9 the Shan−Chen model,10 and the free-energy model or phase-field (PF) model.11 The color model is the first LBM multiphase model, but the recolor step is time-consuming and reduces the computational efficiency; another problem is the anisotropic surface tension that the model generates.12 The Shan−Chen model has been widely used for its simplicity and high computational efficiency.13 Sankaranarayanan et al.14 established some benchmark tests by applying the Shan−Chen model to bubbly flow and the bubble rise velocity agreed well with empirical correlations. Sankaranarayanan et al.15 used an implicit Shan−Chen model to obtain the closure for drag and virtual mass coefficients. Sankaranarayanan and Sundaresan16 proposed the closure relations for the lift coefficients for 2D and 3D bubbles in periodic boxes by implicit LBM. Yu and Fan3 incorporated the technique of adaptive mesh refinement (AMR) into the Shan−Chen model; this way, the interface can be accurately captured using smaller spatial steps around the interface. However, the Shan−Chen model was reported to have several drawbacks, e.g., the local momentum and energy conservation cannot be guaranteed.17 The PF model developed by Swift et al.11 can ensure the local momentum conservation and is thermodynamically consistent. Inamuro et al.18 developed a new PF model to simulate gas−liquid system at a density ratio of ∼1:1000, but the pressure Poisson equation was required to solve in the model. Lee and Lin19 proposed a more stable PF model at high density ratio by solving the pressure evolution equation. However, all of these PF models cannot completely recover the Cahn−Hilliard equation.20 Zheng et al.20 developed a PF model that can simulate a gas−liquid system with a density ratio of ∼1:1000. The PF model proposed by Zheng et al.20 can recover the Cahn−Hilliard equation with second-order accuracy. Ghosh et al.21 used the PF model20 to investigate the interfacial behavior of bubbles for different magnitudes of physical properties and the bubble interactions by 2D simulation. All of the explicit LBM methods for multiphase flow were based on lattice Bhatnagar−Gross−Krook (BGK) model and limited to low Re number, relatively high Mo number conditions, which hardly appear in practical industrial flows. For example, Mo is usually in the order of 10−11 and Re is up to 1000 or more in real cases. Although the implicit method developed by Sankaranarayanan et al.15 can simulate bubbles at larger Re (up to 400) and Mo (as low as 10−9) numbers, these implicit schemes may decrease the computational efficiency.22 To overcome the numerical stability problem at low kinematic viscosity in the lattice BGK model, Lallemand and Luo23 developed a multiple-relaxation-time (MRT) lattice Boltzmann model that was proven to be more stable at high Re number or low kinematic viscosity conditions than the lattice BGK method. Recently, Yu and Fan22 developed a MRT-based Shan−Chen LBM model to enhance the numerical stability at high Re number (up to 1000) and low Mo number (at the order of 10−11). Yu et al.24 further incorporated AMR and MRT into the Shan−Chen model to investigate the multiple bubble interactions. Fakhari and Rahimian25 developed a MRT-based PF model, demonstrating its ability to simulate gas−liquid twophase flows at high Re number (up to 1000) and low Mo number (low to 10−10). Hazi et al.17 pointed out that the original Shan−Chen and PF models cannot simulate phases with different

2. MODEL DESCRIPTION LBM, as a mesoscopic method, can be recovered to the macroscopic fluid dynamic equations through the multiscale Chapman−Enskog expansion. Thus, LBM can be simply viewed as a special method for solving the macroscopic flow equation as well. The PF model is an interface capturing method in DNS. In this section, we first introduce the PF model for interface capturing within the framework of macroscopic fluid dynamic equations, and then we replace the macroscopic equations with LBM models. 2.1. Phase-Field (PF) Model. In the PF model, the Cahn− Hilliard equation is used to capture the gas−liquid interface with finite width interface, whereas the one-field Navier−Stokes equations are used to obtain the flow details, as shown below: Cahn−Hilliard equation:20,25,26 ∂ϕ + ∇·(ϕ u) = θM∇2 μϕ ∂t

(1)

Navier−Stokes equations:

11392

∂ρ + ∇·(ρ u) = 0 ∂t

(2)

∂ρ u + ∇·(ρ uu) = −∇·P + μ∇2 u + Fb ∂t

(3)

dx.doi.org/10.1021/ie303486y | Ind. Eng. Chem. Res. 2013, 52, 11391−11403


Article

where θM is the mobility, μϕ the chemical potential, P the pressure tensor, μ the viscosity, Fb the body force, and ϕ the order parameter. The density (ρ) is defined as ρ = (ρ1 + ρg)/2; here, ρl and ρg denote the density of the liquid and gas, respectively. To reduce the spurious velocity in the calculation, a potential form of −∇·P is adopted:19,20 −∇·P = −ϕ∇μϕ − ∇p0

where p is the pressure and formulated as p = A(3ϕ4 − 2ϕ2ϕ*2 − ϕ*4 ) − kϕ∇2 ϕ −

(14)

2.2. Phase-Field-Based LBM for Multiphase Flow. 2.2.1. Single Relaxation Time. Since MRT can be derived from the LBM method, using a single relaxation time (e.g., the lattice BGK model), the lattice BGK model is introduced first in this section. LBM is originated from the Boltzmann equation,

(4)

Hence, eqs 2 and 3 can be rewritten as ∂ρ + ∇·(ρ u) = 0 ∂t

∂f + ξ·∇x f + a ·∇ξ f = Ω( f ) ∂t

(5)

∂ρ u + ∇·(ρ uu) = −∇(p0 + ϕμϕ ) + μϕ ∇ϕ + μ∇2 u + Fb ∂t

∫ ε dV = ∫ (ψ (ϕ) + 2k (∇ϕ)2 + 13 ρ ln ρ) dV

(6)

∂f f eq − f + ξ·∇x f + a ·∇ξ f = ∂t τc

(7)

(8)

Here, A is also related to the surface tension σ and interface width W. ϕ* = (ρl − ρg)/2 and ±ϕ* are constants corresponding to two equilibrium states. The chemical potential then could be deduced from eq 7:20 μϕ =

∂ε ∂ε − ∇· = A(4ϕ3 − 4ϕ*2ϕ) − k∇2 ϕ ∂ϕ ∂∇ϕ

F⎛ξ − u⎞ a∇ξ f ≈ a∇ξ f eq = − ⎜ 2 ⎟ f eq ρ ⎝ cs ⎠

(9)

The parameters k and A can be easily obtained from macroscopic parameters (i.e., surface tension and interface width). The k and A are computed as follows: k=

A=

3 ⎛ Wσ ⎞ ⎟ ⎜ 8 ⎝ ϕ*2 ⎠

(17)

It has not yet been utilized for the phase-field multiphase LBM simulation and therefore employed in this study for its simplicity. Discretizing the Boltzmann equation with the BGK model with respect to time, location, and velocity spaces leads to the lattice BGK equation30,31 for single-phase flow:

(10)

3⎛ σ ⎞ ⎜ ⎟ 4 ⎝ Wϕ*4 ⎠

(16)

where τc denotes the relaxation time and f eq is the equilibrium distribution function. The early study of LBM simulation dealt with single-phase flows without external force fields. To take the force term into consideration, Shan and Chen10 modified the equilibrium distribution function. A more straightforward method is to put a force term in the evolution equation, but it is difficult to discretize the term a∇ξ f directly in LBM. Buick and Greated27 analyzed the force models that originated from the lattice gas cellular automaton model28 and the force model of Shan and Chen,10 and they proposed a mixed force model by adding an additional term into the evolution equation. He et al.29 proposed an approximation for the force term a∇ξ f with the following formula for single-phase simulation:

Here, k is related to the surface tension (σ) and the interface width (W), ε represents the free-energy density function, and ψ(ϕ) represents the bulk free energy density and has great influence on the interface behavior. ψ(ϕ) is chosen to be a double-well form in this study:20,26 ψ (ϕ) = A(ϕ2 − ϕ*2 )2

(15)

where f denotes the distribution function, ξ the velocity, a the acceleration, and Ω(f) the collision operator. However, the application of the Boltzmann equation is restricted due to the complex nonlinear integration term Ω(f), i.e., the collision operator. To simplify the collision operator, a linear BGK model or single relaxation time model has been proposed and has become the most widely used simplified collision model. The Boltzmann equation with the BGK model is

p0 is the static pressure (p0 = ρcs2), and cs is the speed of sound. Numerically, −∇(p0 + ϕμϕ) can be treated as a pressure term, and μϕ∇ϕ + Fb can be treated as a force term. The chemical potential (μϕ) is required to numerically solve eqs 1, 5, and 6 and is related to the free energy E, which is defined by the expression20 E=

ρ k |∇ϕ|2 + 2 3

fi (x + eiΔt , t + Δt ) − fi (x, t ) (11)

=

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

The profile of order parameter along the direction x perpendicular to the interface can be given analytically as26 ⎛ 2x ⎞ ϕ = ϕ* tanh⎜ ⎟ ⎝W ⎠

fi = fi −

(12)

Pij = ϕμϕ − ε = pδij − k[(∇ϕ) δij − ∇i ϕ∇j ϕ]

1 ⎡ (ei − u) ·F ⎤ eq 1 ⎡ (e − u) ·F ⎤ eq ⎥ fi ⎥ f i = fi − ⎢ i 2 ⎢ 2 2⎣ 2⎣ ρcs ρcs ⎦ ⎦ (19)

This analytical solution will be used to validate the simulation later. The pressure tensor in eq 3 can be calculated from20 2

τ

⎛ 1 ⎞⎟ (ei − u)·F eq + ⎜1 − f i (x, t )Δt ⎝ 2τ ⎠ ρcs 2 (18)

⎡ e ·u (e ·u)2 u2 ⎤ f ieq = wiρ⎢1 + i 2 + i 4 − 2 ⎥ 2cs cs 2cs ⎦ ⎣

(13) 11393

(20)



Article

Here, ei denotes the discrete velocity and wi is the weight coefficient. τ is the dimensionless relaxation time. ei is defined as follows in the D2Q9 and D3Q15 lattice models.

The coefficient Bi is defined as ⎧ w−1[ρ − (1 − w )(ρc 2 + ϕμ )/c 2] i = 0 0 s s ⎪ 0 ϕ Bi = ⎨ ⎪(ρcs 2 + ϕμϕ )/cs 2 i>0 ⎩

D2Q9 model: ⎧(0, 0) i=0 ⎪ ei = ⎨ e( ±1, 0), e(0, ±1) i = 1−4 ⎪ ⎪ e( ±1, ±1) i = 5−8 ⎩

2.2.2. Multiple Relaxation Time. Despite its simplicity, application of the lattice BGK model is limited, because of its numerical instability at low kinematic viscosity. Compared to the single relaxation time for all physical quantities in the lattice BGK model, different relaxation times are used for different physical quantities in the MRT-LBM model to approximate the complex collision operator in the Boltzmann equation. The MRT-LBM model has a clear physical background and is more numerically stable.23 The MRTLBM equation is expressed as31

(21)

D3Q15 model: ⎧(0, 0) i=0 ⎪ ei = ⎨ e(± 1, 0, 0), e(0, ± 1, 0), e(0, 0, ± 1) i = 1−6 ⎪ ⎪ e(± 1, ± 1, ± 1) i = 7−14 ⎩

(22)

fi ̅ (x + eiΔt , t + Δt ) − fi ̅ (x , t )

Here, e = Δx/Δt and represents the lattice velocity, and cs = e/√3. Δx and Δt are the unit lattice space and time steps, respectively. w0 = 4/9, w1−4 = 1/9, and w5−8 = 1/36 for the D2Q9 model. w0 = 2/9, w1−6 = 1/9, and w7−14 = 1/72 for the D3Q15 model. D2Q9 and D3Q15 lattice models are implemented for 2D and 3D simulations in this work, respectively. To recover the Cahn−Hilliard equation in the PF model (eq 1), the governing equation in the lattice BGK model can be rewritten as gi(x + eiΔt , t + Δt ) = gi(x , t ) +

⎛ 1 ⎞ = Λij[f jeq − fj̅ ] + ⎜Iij − Λij⎟Sj ⎝ 2 ⎠

f ̂ = M· f ̅

(23)

where gi represents the distribution function for the order parameter, and the equilibrium distribution function geq i (x,t) is given as32

fi ̂ (x + eiΔt , t + Δt ) eq ⎛ 1 ⎞ = fi ̂ (x , t ) + Λ̂ij[f ĵ − fĵ ] + ⎜Iij − Λ̂ij⎟Sĵ ⎝ 2 ⎠

(25)

where the parameter θ is related to the mobility θM and θM = (τg − 1/2)θΔt. To recover the macroscopic fluid flow equations in the PF multiphase model (eqs 2 and 3), the governing equation in the lattice BGK model (eq 18) is rewritten as

gî (x + eiΔt , t + Δt ) = gî (x , t ) + Λ̂ij[gĵeq – gĵ ]

ni̅ (x + eiΔt , t + Δt ) = ni̅ (x , t ) +

(32)

eq

where ĝi = Mijgj. The equilibrium moment ĝ is obtained by multiplying the transformation matrix and eq 24 together. The details of ĝeq are given in the Appendix. The macroscopic variable ϕ then can be obtained from ϕ = ∑igi. The collision matrix Λ̂ij is related to the dimensionless relaxation time τg in eq 23:

t ) − ni̅ (x , t )) τf

⎛ 1 ⎞⎟ (ei − u) ·F eq + ⎜⎜1 − ni (x , t )Δt 2τf ⎟⎠ ρcs 2 ⎝

where F = Fb + μϕ∇ϕ and

(31)

eq where f ̂ eq i = Mij f j and represents the equilibrium moment. ̂ Λij = diag(s1,s2,...,sb) and denotes the collision matrix (dialogue matrix) and b is the number of lattice velocity directions (b = 9 for D2Q9 model, b = 15 for D3Q15 model). If si in the collision matrix is set as 1/τ, the MRT-LBM model is reduced to the lattice BGK model. Therefore, the lattice BGK model can be viewed as a special case of MRT. Iij is the identity matrix; Sî then denotes the force term in moment space (Ŝi = MijSj). The MRT-LBM equation for the Cahn−Hilliard equation then can be reconstructed as

(24)

Here, Ai is defined as

(nieq (x ,

(30)

where f ̅ = (f0̅ , f1̅ , ..., fb−1 ̅ ) and b is the number of lattice velocity directions. The transformation matrix M for the D2Q9 and D3Q15 lattice models is given in the Appendix. Then, the MRT-LBM equation is rewritten as31

τg

⎧ w−1[ϕ − (1 − w )θμ /c 2] i = 0 0 ⎪ 0 ϕ s Ai = ⎨ ⎪ θμϕ /cs 2 i>0 ⎩

(29)

Here, Λij is the collision matrix and Iij is the identity matrix. Sj denotes the force term Sj = [(ej − u)·F/(ρcs2)] f eq j Δt. To solve eq 29, the distribution function is transferred into the moment space,

gieq (x , t ) − gi(x , t )

⎛ ϕei·u ⎞ ⎟ gieq (x , t ) = wi⎜Ai + cs 2 ⎠ ⎝

(28)

neq i (x,t)

(26) 32

is given as

⎧ ⎡ e ·u (e ·u)2 u 2 ⎤⎫ nieq = wi ⎨Bi + ρ⎢ i 2 + i 4 − 2 ⎥⎬ 2cs 2cs ⎦⎭ ⎣ cs ⎩ ⎪

⎪

⎪

⎪

(27) 11394

τg =

1 1 = s8 s9

τg =

1 1 1 1 1 = = = = s10 s11 s12 s13 s14

(D2Q9 model)

(D3Q15 model)



Article

The MRT-LBM model for momentum equations is reconstructed as nî (x + eiΔt , t + Δt ) ⎛ 1 ⎞ = nî (x, t ) + Λ̂ij(nĵeq − nĵ ) + ⎜Iij − Λ̂ij⎟Fĵ Δt ⎝ 2 ⎠

(33)

where n̂i = Mijn̅j. The equilibrium moment n̂eq is obtained by multiplying the transformation matrix and eq 27 together. The details of n̂eq and the force term in the moment space are given in the Appendix. The macroscopic variables then are calculated from ρ = ∑i ni̅ ρ u = ∑i ni̅ ei +

1 FΔt 2

The kinematic viscosity is related to the collision matrix Λ̂ij: ⎛1 ⎛ 1⎞ 1⎞ v = cs 2⎜ − ⎟Δt = cs 2⎜τf − ⎟Δt ⎝ 2⎠ 2⎠ ⎝ s8

Figure 1. Contour of order parameter for a static bubble.

and

s8 = s9

(D2Q9 model)

⎛1 ⎛ 1⎞ 1⎞ v = cs 2⎜ − ⎟Δt = cs 2⎜τf − ⎟Δt ⎝ 2 2⎠ s ⎝ 10 ⎠ s10 = s11 = s12 = s13 = s14

and

(D3Q15 model)

The 2D simulation was run on the single core of Intel (R) Xeon E5430 CPU (RAM 16G). It takes ∼8.5 min for 10 000 time steps to update the 100 × 100 computational domain. For the 3D simulation, the computational domain is decomposed in one dimension, and each subdomain is run on a single CPU core. The communications between the adjacent subdomains are implemented by a Message Passing Interface (MPI). In this work, the computation domain is 144 × 144 × 288 for 3D simulations. The code is run on the super computer platform “Mole-8.5” located at Institute of Process Engineering of Chinese Academy of Sciences. Seventy two (72) CPU cores are used in this 3D simulation, and it takes ∼55 min to update the total of ∼6 million nodes for 10 000 time steps.

Figure 2. Order parameter distribution across the bubble interface: simulation versus analytical solution.

W = 5. The relaxation time τf is set as 0.5001 for liquid and 0.6 for gas; therefore, the kinematic viscosity ratio of two phases is 1:1000, while τg is set as 0.51. The Laplace law for a 2D bubble is given as σ Δp = (34) R Figure 3 illustrates the change of pressure difference between liquid and gas with the inverse of bubble radius for simulation results and Laplace law prediction for bubbles of different sizes and surface tension. The simulation is almost consistent with the Laplace law prediction for a bubble of large diameter. When the bubble diameter is 12 ⎩ Kb0Mo

(41)

The adjustable parameters n0, c1, and Kb0 are pertinent to the physical properties of the given system. For an air−distilledwater system, n0 = 1.6, c1 = 1.2, and Kb0 = 14.7. For an air− tap-water system, n0 = 1.6, c1 = 0.8, and Kb0 = 14.7. In this study, simulation is carried out for different bubble sizes to compare with the correlation. The bubble diameter varies over a range of 40−50 lattice units. The Mo number is kept as 2.62 × 10−11 for all cases. ρl is set as 1000, and ρg is set as 1. To eliminate the boundary effect, the width of the computational domain is 5 times greater than the bubble diameter and the length of computational domain is 4 times greater than the width. Periodic boundary conditions are imposed for all directions. The implementation of body force is similar to the previous cases. Since the pathway of the bubble is spiral or oscillating, the rise velocity is calculated before the pathway oscillation when the bubble diameter is >1.4 mm. Figure 11 compares the simulation with the correlation of 11398



Article

is symmetrical. The orientation of bubbles changes in the rising process. Figure 14 further compares the flow field around the rising bubbles at two different time steps. When the bubbles start to

Figure 14. Flow field around the two rising bubbles at different time steps: (a) T* = 10.4 in an inspection zone 88 × 80 and (b) T* = 14.3 in an inspection zone 116 × 103.

repel from each other (Figure 14a), one can clearly observe a pair of symmetric vortexes located in the rear and inner sides of bubbles. The liquid behind the bubbles tends to flow upward to the inner side of the bubble pair. Figure 14b shows the approach of two bubbles. A pair of very small symmetric vortexes near the rear and outer side of bubbles can be observed. The liquid behind the bubbles tends to flow upward to the outer side of the bubble pair. The variation of this microscale flow field and the vortex in the rising procedure influences the micromixing and leads to a change in the vertical velocity gradient and, thus, the direction of the lift force acting on the bubbles. Such a simulation not only reproduces the experimental phenomenon, but can be used to understand the mechanism of micromixing and deduce the constitutive law of lift force. Further study is required in this aspect.

Figure 12. Rise of two bubbles of the same size in an array: the trailing one catches up with the leading one.

6. RISE OF BUBBLE SWARM Bubbles are always present in the form of a swarm in industrial applications. The rise of 50 monodispersed bubbles is simulated with the LBM model. The gas volume fraction is ∼3.5%. The bubble diameter is 2 mm initially, and Eo = 0.54 and Mo = 2.62 × 10−11. The computational domain is 144 × 144 × 288. Initially, bubbles are randomly distributed and the distance between the centers of two nearest bubbles is larger than 40 lattice units. Periodic boundary condition is imposed for all directions. Figure 15a illustrates the snapshots of the global view of the rising bubble swarm at different time steps and Figure 15b shows the bubbles in the zone between the two shade planes shown in Figure 15a. Coalescence of two bubbles can be observed. Figure 15b also shows the motion of three pairs of adjacent bubbles, i.e., pair A, pair B, and pair C. The bubble marked as “3” is shared by pair B and pair C, and the vertical distance between two bubbles in pair B is larger than the horizontal distance of two bubbles of pair C. The two bubbles in pair A and pair B initially get closer respectively. Then, in pair B, the trailing bubble catches up with the leading one and merges into a single one instantaneously. However, the trailing bubble of pair A catches up with the leading one and outstrips the leading one. The trailing bubble of pair A is on the left side of the leading one in the first two insets, moving to the right side of the leading one in the third and then surpassing the leading one in the end. The horizontal distance between the two bubbles of pair C changes little at first and becomes larger with time. Despite the fact that the initial distance between the two

Figure 13. Rising trajectory of two horizontally aligned bubbles of same size: (a) LBM simulation and (b) experimental study.34

the Mo number is equal to 1.06 × 10−9. Periodic boundary conditions are imposed at the top and bottom, and bounceback boundary conditions are applied to the left and right. As shown in Figure 13a, the bubbles do approach and bounce back repeatedly during the rising process, which reproduces the experimental phenomenon. The moving path of the bubble pair 11399



Article

Figure 15. Snapshots of rising bubble swarm at different time steps.

Figure 16. Local flow field around the three pairs of bubbles at different time steps: (A) pair A, (B) pair B, and (C) pair C.

bubbles in pair C is smaller, they almost have no tendency to contact and merge. Apparently the fate of the three pairs of bubbles is not only influenced by the two bubbles in the pair, but also by the complex surrounding the flow field induced by other bubbles outside the pairs.

Figure 16 illustrates the local streamline around the three bubble pairs corresponding to those in Figure 15b. There are some common streamline “string” connecting the two bubbles in pairs A and B at T* = 3.7 and T* = 7.9, implying the attraction tendency between these bubbles in pair. But the situation for pair 11400



Article

A is more complicated. The position variation and final separation of two bubbles in pair A may be caused by the flow instability and other adjacent bubbles outside the inspection zone. No common streamlines exist for the two bubbles in pair C. The lateral motion of bubbles is evident and the number of bubbles in a fixed inspection window changes from time to time. It can be concluded that the interaction between the vertically aligned adjacent bubbles is stronger than that of horizontally aligned adjacent bubbles.

or horizontally aligned bubbles are consistent with literature reports and the repeatedly attraction and repulsion process of two ellipsoidal bubbles is caused by the nearest bubble-induced vortexes. The change of position between vortexes and bubbles may alter the direction of the lift force acting on the two bubbles. Bubble coalescence is observed in a bubble swarm with 50 bubbles of 2 mm. The interaction between vertically aligned bubbles is stronger than those of horizontally distributed bubbles and, hence, induces a greater tendency for bubble coalescence. This study also provides a microscale or mesoscale basis for the development and validation of the stability-constrained multifluid CFD models for gas-liquid two-phase flows.35−37

7. CONCLUSIONS A Multiple-Relaxation-Time lattice Boltzmann method (MRT-LBM) for phase-field models is developed for bubble dynamics simulation. The model can simulate the gas−liquid flow with a large kinematic viscosity ratio and bubbles with higher Reynolds number and lower Morton number. First, the model is validated with analytical solutions of some classical phenomena, such as the order parameter distribution and the pressure difference across the bubble surface and liquid, as well as the capillary wave. Second, the model proves to be able to qualitatively simulate the shape variation of bubbles at different characteristic parameters. The pathway oscillation of a 2-mm air bubble rising in water is observed in the simulation. In contrast to the 2-mm bubble, the rise velocity of a 5-mm air bubble always fluctuates when reaching the terminal velocity. It is essentially due to the constant shape oscillation of the 5-mm air bubble in the rising procedure before the onset of spiral or zigzag motion. The simulation results of terminal velocity of different sizes of air bubbles in air−water systems are in reasonable agreement with experimental correlations. Third, the simulation results of the interaction between two vertically

■

APPENDIX

For 2D simulation (D2Q9 model)

The transformation matrix for the D2Q9 model is ⎡ 1 1 1 1 1 1 1 1 1⎤ ⎥ ⎢ ⎢− 4 − 1 − 1 − 1 − 1 2 2 2 2 ⎥ ⎢ 4 − 2 − 2 − 2 − 2 1 1 1 1⎥ ⎥ ⎢ ⎢ 0 1 0 − 1 0 1 − 1 − 1 1⎥ M = ⎢ 0 − 2 0 2 0 1 − 1 − 1 1⎥ ⎥ ⎢ ⎢ 0 0 1 0 − 1 1 1 − 1 − 1⎥ ⎢ 0 0 − 2 0 2 1 1 − 1 − 1⎥ ⎥ ⎢ ⎢ 0 1 −1 1 −1 0 0 0 0 ⎥ ⎣ 0 0 0 0 0 1 − 1 1 − 1⎦

Equilibrium moment ĝeq: ĝ eq = ϕ(1, − 4 + 6μϕ θ /ϕ , 4 − 9μϕ θ /ϕ , ux , − ux , uy , − uy , 0, 0)T

Equilibrium moment n̂eq: n̂ eq = ρ(1, −2 + 3u 2 + 6ϕμϕ /ρ , 1 − 3u 2 − 9ϕμϕ /ρ , ux , −ux , uy , −uy , ux 2 − uy 2 , uxuy)T

The force term in moment space:

For 3D simulations (D3Q15)

F̂ = (0, 6F ·u, − 6F ·u, Fx , − Fx , Fy , − Fy , 2(Fxux − Fyuy), Fyux + Fxuy)T

⎡ 1 1 1 1 1 ⎢− 2 − 1 − 1 − 1 − 1 ⎢ ⎢ 16 −4 −4 −4 −4 ⎢ 0 1 −1 0 0 ⎢ 0 −4 4 0 0 ⎢ ⎢ 0 0 0 1 −1 ⎢ 0 0 0 −4 4 M=⎢ 0 0 0 0 0 ⎢ 0 0 0 0 0 ⎢ ⎢ 0 2 2 −1 −1 ⎢ 0 0 0 1 1 ⎢ 0 0 0 0 0 ⎢ ⎢ 0 0 0 0 0 ⎢ 0 0 0 0 0 ⎣ 0 0 0 0 0

1 −1 −4 0 0 0 0 1 −4 −1 −1 0 0 0 0

The transformation matrix for D3Q15 model is 1 −1 −4 0 0 0 0 −1 4 −1 −1 0 0 0 0

11401

1 1 1 1 1 1 1 1 1 0 0 1 1 1 1

1 1 1 −1 −1 1 1 1 1 0 0 −1 1 −1 −1

1 1 1 1 1 −1 −1 1 1 0 0 −1 −1 1 −1

1 1 1 −1 −1 −1 −1 1 1 0 0 1 −1 −1 1

1 1 1 1 1 1 1 −1 −1 0 0 1 −1 −1 −1

1 1 1 −1 −1 1 1 −1 −1 0 0 −1 −1 1 1

1 1 1 1 1 −1 −1 −1 −1 0 0 −1 1 −1 1

1⎤ 1⎥ ⎥ 1⎥ −1⎥ −1⎥ ⎥ −1⎥ −1⎥ −1⎥ −1⎥ ⎥ 0⎥ 0⎥ 1⎥ ⎥ 1⎥ 1⎥ −1⎦



Article

Equilibrium moment ĝeq: ⎛ ⎞T 7 7 7 ĝ eq = ϕ⎜1, 3μϕ θ /ϕ − 2, −45μϕ θ /ϕ + 16, ux , − ux , uy , − uy , uz , − uz , 0, 0, 0, 0, 0, 0⎟ ⎝ ⎠ 3 3 3

Equilibrium moment n̂eq: ⎛ ⎞T 7 7 7 n̂ eq = ρ⎜1, − 1 + u 2 + 3ϕμϕ /ρ , 1 − 5u 2 − 45ϕμϕ /ρ , ux , − ux , uy , − uy , uz , − uz , 3ux2 − u 2 , uy2 − uz2 , uxuy , uyuz , uzux , 0⎟ ⎝ ⎠ 3 3 3

The force term in moment space: F̂ = ⎛ ⎞T 7 7 7 ⎜0, 2F · u , − 10F · u , F , − F , F , − F , F , − F , 6F u − 2F · u , 2F u − 2F u , F u + F u , F u + F u , F u + F u , 0⎟ x x y y z z x x y y z z y x x y y z z y z x x z ⎝ ⎠ 3 3 3

■

k = model parameter related to the surface tension and interface width kw = wavenumber Kb0 = parameter in experimental correlation Mij = transformation matrix Mo = Morton number; Mo = g(ρl − ρg)ρl2v4/σ3 n̅i = distribution function for the N−S equations nieq = equilibrium distribution function for the N−S equations n̂i = distribution function in moment space for the N−S equations n̂eq i = equilibrium distribution function in moment space for the N−S equations n0 = parameters in experimental correlation p0 = static pressure p = pressure Δp = pressure difference between gas and liquid P = pressure tensor R = radius of bubble Re = Reynolds number: Re = ud,tde/v si = adjustable parameters in collision matrix Sj = force terms t* = number of computational time steps T* = dimensionless time ud,t = terminal velocity of a bubble W = interface width We = Weber number; We = (ρdeu2)/(σ) Δx = spatial step Δt = temporal step X = position of lattice node along the x-axis perpendicular to the interface LX = length of computational domain along the x-axis

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors wish to thank the long-term support from the National Natural Science Foundation of China (Nos. 21222603, U1162107), the Ministry of Science and Technology of China (2013BAC12B01, 2009CB219906), and the Strategic Priority Research Program of Chinese Academy of Sciences (No. XDA07080301).

■

NOMENCLATURE A = model parameter related to the surface tension and interface width a = acceleration vector b = number of lattice velocity directions c1 = parameter in experimental correlation cs = speed of sound de = equivalent bubble diameter E = free energy of fluid Eo = Eötvös number: Eo = g(ρl − ρg)d2e /σ e = magnitude of unit lattice of velocity vector ei = lattice velocity vector f i = distribution function f eq i = equilibrium distribution function Fb = body force F̂ = force in moment space g = gravitational acceleration constant gi = distribution function for the C−H equation geq i = equilibrium distribution function for the C−H equation ĝi = distribution function in moment space for the C−H equation ĝeq i = equilibrium distribution function in moment space for the C−H equation H = height of computational domain in lattice units Iij = identity matrix

Greek Letters

σ = surface tension ε = free energy density function θ = LBM model parameter related to the mobility θM θM = mobility Λij = collision matrix Λ̂ij = collision matrix in moment space μ = viscosity 11402



Article

μϕ = chemical potential ν = kinematic viscosity ξ = velocity in continuous space ρ = (ρl + ρg)/2 ρg = density of gas ρl = density of liquid τ = dimensionless relaxation time τc = relaxation time τf = dimensionless relaxation time for fluid flow τg = dimensionless relaxation time for C−H equation ϕ = order parameter ±ϕ* = equilibrium order parameter Ω(f) = collision operator wi = weight coefficient

■

(20) Zheng, H. W.; Shu, C.; Chew, Y. T. A lattice Boltzmann model for multiphase flows with large density ratio. J. Comput. Phys. 2006, 218 (1), 353−371. (21) Ghosh, S.; Das, A. K.; Vaidya, A. A.; Mishra, S. C.; Das, P. K. Numerical Study of Dynamics of Bubbles Using Lattice Boltzmann Method. Ind. Eng. Chem. Res. 2012, 51 (18), 6364−6376. (22) Yu, Z.; Fan, L.-S., Multirelaxation-time interaction-potentialbased lattice Boltzmann model for two-phase flow. Phys. Rev. E 2010, 82, (4). (23) Lallemand, P.; Luo, L. S. Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability. Phys. Rev. E 2000, 61 (6), 6546−6562. (24) Yu, Z.; Yang, H.; Fan, L.-S. Numerical simulation of bubble interactions using an adaptive lattice Boltzmann method. Chem. Eng. Sci. 2011, 66 (14), 3441−3451. (25) Fakhari, A.; Rahimian, M. H., Phase-field modeling by the method of lattice Boltzmann equations. Phys. Rev. E 2010, 81, (3). (26) Jacqmin, D. Calculation of two-phase Navier−Stokes flows using phase-field modeling. J. Comput. Phys. 1999, 155 (1), 96−127. (27) Buick, J. M.; Greated, C. A. Gravity in a lattice Boltzmann model. Phys. Rev. E 2000, 61 (5), 5307−5320. (28) Frisch, U.; d’Humières, D.; Hasslacher, B.; Lallemand, P.; Pomeau, Y.; Rivet, J. P. Lattice gas hydrodynamics in two and three dimensions. Complex Syst. 1987, 1, 649−707. (29) He, X. Y.; Shan, X. W.; Doolen, G. D. Discrete Boltzmann equation model for nonideal gases. Phys. Rev. E 1998, 57 (1), R13− R16. (30) He, X.; Chen, S.; Zhang, R. A Lattice Boltzmann Scheme for Incompressible Multiphase Flow and Its Application in Simulation of Rayleigh−Taylor Instability. J. Comput. Phys. 1999, 152 (2), 642−663. (31) Premnath, K. N.; Abraham, J. Three-dimensional multirelaxation time (MRT) lattice-Boltzmann models for multiphase flow. J. Comput. Phys. 2007, 224 (2), 539−559. (32) Huang, J. J.; Shu, C.; Chew, Y. T. Mobility-dependent bifurcations in capillarity-driven two-phase fluid systems by using a lattice Boltzmann phase-field model. Int. J. Numer. Methods Fluids 2009, 60 (2), 203−225. (33) Cheng, M.; Hua, J.; Lou, J. Simulation of bubble−bubble interaction using a lattice Boltzmann method. Comput. Fluids 2010, 39 (2), 260−270. (34) Sanada, T.; Sato, A.; Shirota, M.; Watanabe, M. Motion and coalescence of a pair of bubbles rising side by side. Chem. Eng. Sci. 2009, 64 (11), 2659−2671. (35) Yang, N.; Chen, J.; Ge, W.; Li, J. A conceptual model for analyzing the stability condition and regime transition in bubble columns. Chem. Eng. Sci. 2010, 65 (1), 517−526. (36) Xiao, Q.; Yang, N.; Li, J. Stability-constrained multi-fluid CFD models for gas-liquid flow in bubble columns. Chem. Eng. Sci. 2013. DOI: 10.1016/j.ces.2013.02.027 (37) Yang, N.; Chen, J.; Wang, Y.; Li, J. Multi-scale analysis of gas− liquid interaction and CFD simulation of gas−liquid flow in bubble columns. Chem. Eng. Sci. 2011, 66 (14), 3212−3222.

REFERENCES

(1) Clift, R.; Grace, J. R.; Weber, M. E. Bubbles, Drops, and Particles; Academic Press: New York, 1978. (2) Fan, L. S.; Tsuchiya, K. Bubble Wake Dynamics in Liquids and Liquid−Solid Suspensions; Butterworth−Heinemann: New York, 1990. (3) Yu, Z.; Fan, L.-S. An interaction potential based lattice Boltzmann method with adaptive mesh refinement (AMR) for two-phase flow simulation. J. Comput. Phys. 2009, 228 (17), 6456−6478. (4) Unverdi, S. O.; Tryggvason, G. A front-tracking method for viscous, incompressible, multi-fluid flows. J. Comput. Phys. 1992, 100 (1), 25−37. (5) Hirt, C. W.; Nichols, B. D. Volume of Fluid (Vof) Method for the Dynamics of Free Boundaries. J. Comput. Phys. 1981, 39 (1), 201−225. (6) Sussman, M.; Smereka, P.; Osher, S. A level set approach for computing solutions to incompressible two-phase flow. J. Comput. Phys. 1994, 114 (1), 146−159. (7) Sussman, M.; Puckett, E. G. A Coupled Level Set and Volume-ofFluid Method for Computing 3D and Axisymmetric Incompressible Two-Phase Flows. J. Comput. Phys. 2000, 162 (2), 301−337. (8) Aidun, C. K.; Clausen, J. R. Lattice-Boltzmann Method for Complex Flows. Annu. Rev. Fluid Mech. 2010, 42 (1), 439−472. (9) Gunstensen, A.; Rothman, D.; Zaleski, S.; Zanetti, G. Lattice Boltzmann model of immiscible fluids. Phys. Rev. A 1991, 43 (8), 4320−4327. (10) Shan, X.; Chen, H. Lattice Boltzmann model for simulating flows with multiple phases and components. Phys. Rev. E 1993, 47 (3), 1815−1819. (11) Swift, M.; Osborn, W.; Yeomans, J. Lattice Boltzmann Simulation of Nonideal Fluids. Phys. Rev. Lett. 1995, 75 (5), 830−833. (12) Chen, S.; Doolen, G. D. Lattice Boltzmann method for fluid flows. Annu. Rev. Fluid Mech. 1998, 30, 329−364. (13) Sukop, M. C.; Thorne, D. T. Lattice Boltzmann Modeling: An Introduction for Geoscientists and Engineers; Springer: Berlin, Heidelberg, Germany, 2007. (14) Sankaranarayanan, K.; Shan, X.; Kevrekidis, I. G.; Sundaresan, S. Bubble flow simulations with the lattice Boltzmann method. Chem. Eng. Sci. 1999, 54 (21), 4817−4823. (15) Sankaranarayanan, K.; Shan, X.; Kevrekidis, I. G.; Sundaresan, S. Analysis of drag and virtual mass forces in bubbly suspensions using an implicit formulation of the lattice Boltzmann method. J. Fluid Mech. 2002, 452. (16) Sankaranarayanan, K.; Sundaresan, S. Lift force in bubbly suspensions. Chem. Eng. Sci. 2002, 57 (17), 3521−3542. (17) Hazi, G.; Imre, A. R.; Mayer, G.; Farkas, I. Lattice Boltzmann methods for two-phase flow modeling. Ann. Nuclear Energy 2002, 29 (12), 1421−1453. (18) Inamuro, T.; Ogata, T.; Tajima, S.; Konishi, N. A lattice Boltzmann method for incompressible two-phase flows with large density differences. J. Comput. Phys. 2004, 198 (2), 628−644. (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 (1), 16−47. 11403


Direct Numerical Simulation of Bubble Dynamics Using Phase-Field

Recommend Documents