Simulating GasâLiquid Flows by Means of a Pseudopotential Lattice

Article pubs.acs.org/IECR

Simulating Gas−Liquid Flows by Means of a Pseudopotential Lattice Boltzmann Method M. R. Kamali and H. E. A. Van den Akker* Department of Chemical Engineering, Delft University of Technology, Delft, Netherlands ABSTRACT: Dispersed gas (vapor)−liquid flow through an inclined microchannel with bends has successfully been simulated, that is, without numerical difficulties, by means of a two-phase Lattice Boltzmann method. Combining in this method the ShanChen1 pseudopotential interaction model with the Yuan and Schaefer2 proposal for dealing with nonideal equations of state makes high density ratios achievable. This approach also allows simulation of gas−liquid flows without explicitly having to track the phase interfaces. Rather, a potential function related to the equation of state for vapor−liquid equilibrium, a coupling strength representing attraction or repulsion between species, and a relaxation time scale take care of microscale and mesoscale phenomena such as phase separation and interfacial tension as well as interphase transport and multiphase flow. In addition, fluid−wall interaction (contact angle) is taken into account by selecting proper potential functions and coupling strengths. As far as the phase behavior is concerned, we assessed our method by studying the phase separation process and by validating against Maxwell’s equilibrium rule. Qualitative validation of our approach of gas−liquid flow has been done with a comparison against experimental data on a single bubble rise. Detailed simulations were carried out for an individual Taylor bubble in a channel, the results of which compared favorably to literature data.

■

INTRODUCTION In the field of chemical (reaction) engineering, the usual way of simulating flow and transport phenomena (and chemical reactions) is in terms of continuum variables and properties and by using the Finite Volume (FV) technique. With FV, the 3-D continuum conservation equations for mass, species, momentum, and heatactually partial differential equationsare approximated by discretized algebraic equations which are solved in a large number of discrete mini control volumes or grid cells. Apart from being known to and favorite among chemical engineers, FV is a direct and appealing translation of the transport equations. The basic assumption is that all variables are conceived as continuum variablesa notion that is very common in chemical engineering. The result of FV simulations are (velocity, pressure, temperature, concentration) fields, or spatial distributions, although in a discretized form. Particularly when the number of transport equations becomes large, such as when several chemically reacting species are of interest and in turbulent single-phase and multiphase flows, or when the computational grid is refined, the above features turn the FV technique into a computationally intensive process with dubious or slow convergence characteristics. In this respect, the Lattice Boltzmann (LB) technique is promising: within a LB code, all information transfer is local in time and in space, to a small set of nearest neighbors only. As a result, LB codes are suitable for massively parallel computations. Further, the LB method is relatively fast due to the low number of operations per grid node and per time step. Consequently, LB methods have started playing an important role as a simulation tool for advanced modeling of both materials properties and flow processes.3 Starting with Eggels4 and Derksen and Van den Akker,5 LB techniques are being used for simulating macroscale turbulent flows.6−9 © XXXX American Chemical Society

Further, LB methods have extensively been invoked for studying the motion of single and multiple solid particles in a fluid, in many cases aimed at finding correlations for drag and lift coefficients.10−12 In addition, LB techniques are increasingly exploited for simulating gas−liquid flows.13−17 The origin of many phenomena in multifluid systems such as phase separation, interfacial tension, and solid surface wettability is in the molecular interactions. Incorporation of these microphysical flow aspects is done rather easily by using the LB approach. The success of LB based methods can largely be attributed to their mesoscopic and kinetic nature, which enables the simulation of the macroscopic multiphase phenomena like interfacial dynamics in terms of the underlying microscopic physics. Typically, interfaces in LB are diffuse: they have a thickness of a few lattice distances. The advantages and challenges of using LB methods have been extensively discussed in previous papers.18−27 There are two main categories of LB models for multiple fluid phases. One category is based on a free energy concept28 which solves the binary system in two steps. In the first step, a set of distribution functions is used to solve the velocity fields while conserving mass and momentum. The calculated velocities are then used for solving an advection−diffusion equation of a certain phase field variable through the Cahn− Hilliard equation29 to find the position of the phase interfaces. Special Issue: Multiscale Structures and Systems in Process Engineering Received: December 5, 2012 Revised: March 9, 2013 Accepted: March 11, 2013

A

dx.doi.org/10.1021/ie303356u | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Industrial & Engineering Chemistry Research

Article

The second category of models is based on considering separate distribution functions for different components in an immiscible mixture. These distribution functions are then coupled through interparticle attraction/repulsion forces. This means that LB particles of the same phase attract each other while repelling the ones in the other phases. These attraction/ repulsion forces are the main sources for the interface formation. This approach is known as the pseudopotential method as proposed by Shan and Chen.1,30,31 One of the most attractive features of the pseudopotential LB method is its ability to automatically detect phase interfaces through the natural density variations across the phase interfaces. This makes these methods more appealing and robust compared to the free energy LB approaches as well as the various FV methods such as the volume of fluid, the level set and the front tracking methods: they all require accurate interface tracking through some extra interface detection algorithms. Such extra algorithms increase the computational costs and limit the application of these methods to rather simplified geometries with just a restricted number of interfaces, for example, bubbles or droplets, in the system. In this paper, we will demonstrate the use of the pseudopotential LB method for the purpose of simulating in 2-D dispersed gas−liquid flow through the complex geometry of an inclined microchannel with two bends in the context of a study as to the use of structured reactors for the Fischer− Tropsch process (see Vervloet et al).32 We will validate our method by applying it to two canonical cases: a 2-D gas bubble freely rising in an unconfined liquid, and the flow of a 2-D Taylor bubble in a microchannel. The presentation of our simulation results is preceded by a detailed discussion on several crucial aspects of this pseudopotential method. Particularly, the role and effect of using different equations of state will be discussed along with their effect on the gas−liquid density ratio. In this context, the coexistence curve obtained by using the Redlich−Kwong equation of state (in the LB framework) will be assessed by means of the Maxwell rule of equal area. Further, we will present snapshots illustrating the process of phase segregation and drop formation starting from a random liquid−vapor field while obeying the Carnahan−Starling equation of state. Most of the paper, however, is devoted to a detailed assessment of the motion of a 2-D Taylor bubble in a channel.

First, the rigid-sphere model is replaced by one that portrays the attractive and repulsive forces between the molecules, such as described by the classical Lennard-Jones potential for nonpolar and spherically symmetric molecules resulting in an intermolecular force. In addition, all molecules do not have the same velocity u but exhibit a velocity distribution f(r,v,t) where the quantity f(r,v,t) dr dv denotes the probable number of molecules being at time t in a volume element dr at position r and having velocities in the range dv about v. The molecular velocity distribution function under equilibrium conditions has been named after Maxwell and Boltzmann and runs as f eq =

ρ ω /2

(2πRT )

⎛ |v − u|2 ⎞ exp⎜ − ⎟ 2RT ⎠ ⎝

(1)

in which u the average velocity of all velocity realizations of the velocity distribution, ω is the number of dimensions, and R is the gas constant. Finally, the notion of the mean free path is abandoned. After all, the Boltzmann equation is introduced to represent the chaotic processes at the molecular scale. It describes how f evolves with time under nonequilibrium conditions, that is, as a result of temperature, velocity, and concentration gradients. This Boltzmann equation runs as ∂t f + v·∇x f + Ftot ·∇v f = (∂t f )collision

(2)

in which ∂t , ∇x , ∇v denote the derivative toward time, space, and velocity, respectively. The above three more relaxed model assumptions lead to the Chapman−Enskog model referred to in section 1.5 of Deen33 and in section 1.4 and Appendix D of Bird et al.34

■

THE CONCEPT OF LATTICE BOLTZMANN Similarly, the mesoscopic description of multicomponent multiphase behavior and the macroscopic description of convective fluid flows may be replaced by a discrete particle based approach. Starting in the 1980s, so-called Lattice Boltzmann (LB) techniques have been developed. The LB methods may be conceived as a bottom-up approach to phase behavior and fluid flow modeling, on the analogy of how kinetic gas theory is at the basis of molecular transport phenomena in simple, dilute gases. In the LB methodology,3,35 fluid flow is mimicked by a system of fictitious or pseudoparticles residing on a regular lattice and exhibiting discrete velocities which can be represented by some velocity distribution function. These pseudoparticles may be considered as coarse-grained groups of fluid molecules carrying some properties of real fluid portions. These pseudoparticles, or fluid portions, are assumed to move at varying speeds in specific, well-defined directions to neighboring lattice sites, or nodes. More specifically, these fictitious particles span a velocity distribution, or probability distribution function, varying in space and time, and therefore also denoted as f(x,v,t). The details of these motions and collisions do not matter. All that matters is an effective and efficient description, as simple as possible, retaining the overall mesoscopic or macroscopic effects, to be realized by reducing the degrees of freedom within this chaotic interaction of fictitious particles without affecting the overall behavior. With a single set of particles, a suitable lattice topology, and proper collision rules, and after a straightforward translation of this chaotic behavior toward the fluid’s velocity and pressure fields, it can be proven that the

■

PARTICLE-BASED METHODS The LB technique is a particle-based method that presumes all fluids can be conceived as consisting of interacting fictitious particles. In classical chemical engineering, continuum models are widely preferred over particle based models, though textbooks on transport phenomena often refer to kinetic gas theory and explain how “random” molecular motions and hardsphere collisions are at the basis of the phenomenological transport coefficients in gases. Deen33 shows that in gases the limiting step for molecular transport is the frequency of the collisions, leading to very similar values for the diffusivities for species, momentum, and heat. The basic model expressing the phenomenological coefficients in terms of the product of a velocity times a path length exploits a rigid sphere model for the molecules as well as the concept of the free mean path during which the velocity of a molecule does not change. More advanced models relax these three rather stringent assumptions of the hard-sphere kinetic gas theory: B



Article

single-phase fluid eventually asymptotically obeys the classical Navier−Stokes equation (at least in the incompressible limit). In the LB methodology, the actual equation to be solved at the microscopic, or mesoscopic, scale is the above Boltzmann equation in which the most complicated part is the “collision” term at the right-hand side representing the microscopic, or mesoscopic, interaction of the fictitious particles. A very common simplification is the Bhatnagar−Gross−Krook (BGK) approach,36 in which the collision or encounter step is conceived as a simple relaxation process with a time constant τ35 and therefore is approximated by (∂t f )collision

f − f eq ≅− τ

tension between the different phases. This is realized by adding an attractive or repulsive tail to the elastic collisions of the fictitious particles making up the two phases. The extra force becomes part of Ftot in eq 5. For a single-component system (such as water−steam), the general expression for this interparticle force in discretized format reads as Fδ (x , t ) = −Gδψδ(x , t ) ∑ wiψδ(x + ciΔt , t )ciΔt i

in which the subscript δ denotes the (single) component δ involved, wi is the weight factor related to the distance from position of lattice node i > 0 to the central node i = 0, ci is the (lattice) vector associated with the lattice direction i (typically connecting the node x to a neighboring node when i is not zero and equal to the null vector for i = 0), ψδ stands for the potential function of the component δ, and Gδ is a constant, named the attraction strength, that represents the mutual interaction of component δ particles. This Gδ is zero for an ideal gas, inherently comprised in the LB formulation and rendering the extra interaction force of eq 6 zero, while a negative value of Gδ expresses attraction in nonideal systems leading to phase separation. As it was originally used by Shan and Chen, the potential function can be an exponential function in the form of ψδ = a exp (−b/ρ). However, other expressions are feasible to represent different equations of state (see further on). The interaction force Fδ expresses the interaction with the neighboring lattice nodes. For nonzero values of Gδ, it is the gradient of the nonideal part of the pressure. In multicomponent systems, the mean field interaction force Fδδ̅ imposed on component δ from itself and all other components δ̅ in the discretized format is defined as

(3)

eq

in which f is very often approximated by using the secondorder Taylor expansion at low Mach numbers.37 In addition, the phase space velocity gradient in the Boltzmann equation is replaced by v − u eq ∇v f ≅ ∇v f eq = − f (4) RT and the Boltzmann equation (eq 2) can be simplified into ∂t f + v·∇x f = −Ftot ·

f − f eq (v − u) eq f − RT τ

(6)

(5)

All external forces such as a pressure gradient or gravity are added to the mesoscale interaction forces between these fictitious particles and comprise the total force per unit of mass Ftot in eq 5. The overall macroscopic moments of the fluid, such as density and velocity, are defined from the distribution functions by integration over velocity space v in ρ = ∫ v f dv and ρu = ∫ vvf dv. The kinematic viscosity ν, which in the continuum description is a measure for the fluid’s resistance to internal stresses, is related to the relaxation time τ in the LB collision term by ν = C2s (τ − 0.5 Δt), while Cs Δt/Δx = (RT)1/2 = 1/√3 denotes the speed of sound. In our simulations, unless stated otherwise, we used reference density ρ0 = 1, discrete space step Δx = 1 and time step Δt = 1 for unit conversion purposes. In the LB approach, the Boltzmann equation, eq 5, is solved on a computational lattice or grid.3,35 Unless stated otherwise, we use the D2Q9 lattice topology in which “2” stands for the dimensionality of the description, while “9” denotes the number of velocity directions seen from the central point i = 0. The most popular lattice in 3-D is D3Q19.

Fδδ ̅ (x , t ) = −ψδ(x , t ) ∑ Gδδ ̅ ∑ wiψδ ̅(x + ciΔt , t )ciΔt i

δ̅

(7)

in which Gδδ̅ is the coupling strength that defines the degree of immiscibility between components, sharpness of the interface, and surface tension. A negative value of Gδδ̅ again stands for an attractive force, while a positive Gδδ̅ value represents a repulsive force between the components. Analogous to the component−component interactions in eq 7, component−wall interactions in the Shan−Chen method are defined as Fδ s(x , t ) = −Gδ sψδ(x , t ) ∑ wS i i(x + ciΔt , t )

■

i

THE PSEUDOPOTENTIAL MULTICOMPONENT TWO-PHASE MODEL The pseudopotential approach, originally due to Shan and Chen, has been extensively used for modeling singlecomponent and multicomponent two-phase problems.18,21,27,38−41 The method considers separate distribution functions for the various components involved, each governed by its own Boltzmann equation. In addition, this method considers an interparticle force that represents the interaction of the components with themselves and with each other. This interparticle force depends on potential functions for each of the components involved which in their turn depend on the density of the components involved through the equation of state. As a further outcome, these interparticle forces provide a means of controlling the phase separation and the interfacial

ψδ s(x + ciΔt , t )ciΔt

(8)

in which the subscript s denotes the solid wall and the variable S is a flag quantity which activates the component−wall interactions in the wall lattice nodes only. By assigning to the wall a certain density ρs and the same EOS as for the component δ, the variable ψδs(x,t) can be defined at solid interfaces. This value together with a proper coupling strength coefficient Gδs defines the total interaction force between component and solid interface. Thus, the wettability of the solid toward a certain component is tuned by ρs and Gδs. For instance, by choosing in a system consisting of the pure component δ the same EOS for the solid wall as for component δ and by putting Gδs = −1 and ρs = ρliq, the solid surface is seen as a liquid and thus the liquid fully wets C



Article

the solid surface and the contact angle is zero. Similarly, choosing Gδs = 1 and ρs = ρvap results in hydrophobicity of the solid and in a nonzero contact angle. By varying these values one may find how the contact angle depends on the density of the wall.41 Yet, a neutral wall boundary condition can be attained by imposing a regular bounce-back rule on the wall nodes while Gδs = Gδδ and while the ψδs values for the wall points are equal to those in the neighboring fluid nodes.

■

THE DENSITY RATIO BETWEEN TWO PHASES The original Shan and Chen method is limited to low density ratios (maximum about 10); for higher values, numerical instabilities will occur. Yuan and Schaefer2 used the potential function in a form which can recover a properly chosen equation of state resulting in an increase in the density ratios between two phases in the system. This concept was developed further by Kupershtokh et al.42 who proposed an alternative implementation of the pseudopotential forces. Pursuing these suggestions, we introduced an expression for the potential function in a single-component environment reflecting some nonideal equation of state (EOS), namely, ψδ = √ 2{Pδ(ρδ)/Cs2 − ρδ}1/2 in which Pδ(ρδ) represents the nonideal EOS. For instance, one can incorporate the Pδ of the Carnahan−Starling (CS) and Redlich−Kwong (RK) EOSs (summarized in Table 1, along with the values for the critical

Figure 1. Liquid and vapor density dependence on the reduced temperature (coexistence curve) for RK EOS. Solid line, curve obtained via the Maxwell equal area rule; circles, single-component flat interface LB simulations by pseudopotential formulation.

liquid to vapor density ratio for a single-component system as obtained from LB simulations can be also determined. As a direct outcome of the pseudopotential force implementation, phase separation occurs. Figure 2 shows in a few snapshots how, starting from a random spatial distribution of vapor and liquid “particles” (not shown), spontaneous phase separation along with droplet formation takes place toward an equilibrium situation corresponding to the EOS imposed. While such a figure has already been presented by Shan and Chen1 for their simulation carried out with an exponential potential function, our Figure 2 has been obtained for the CS EOS, resulting in a density ratio of around 1000. In principle, this figure can be obtained for any other nonideal EOS. It should be added, however, that multicomponent simulations are less stable than single-component simulations; as a result, the flexibility in selecting proper values for the various parameters is smaller and the achievable density ratio is lower in multicomponent simulations.

Table 1. Redlich-Kwong and Carnahan−Starling EOS and the Pertinent Values Used in the LB Simulations values used in the LB simulations

equation of state Carnahan−Starling:

Pδ = ρδ RT

1 + bρδ / 4 + (bρδ / 4)2 − (bρδ / 4)3 (1 − bρδ / 4)3

− aρδ2

Tc = 0.3773a/bR, Pc = 0.0706/b2, ρc = 0.47644/b Redlich−Kwong:

Pδ =

ρδ RT 1 − bρδ

−

aρδ2 T (1 + bρδ )

Tc = (0.2027a/bR)2/3, Pc = 0.0299R1/3a2/3/b5/3, ρc = 0.260/b

a = 1, b = 4, R = 1; Tc = 0.094325, Pc = 0.00441, ρc = 0.11911

■

a = 2/49, b = 2/21, R = 1; Tc = 0.1961, Pc = 0.1785, ρc = 2.73

SIMULATION OF MOTION OF A DROPLET IN AN INCLINED MICROCHANNEL The current pseudopotential LB method was applied for simulating the dispersed two-phase flow in an inclined microchannel with two expanding bends (represented from an industrial microstructured reactor). 32 Two different simulations were carried out to reproduce the effect of a different degree of wettability of the channel wall. Thus, component−wall interactions are tuned such that the static contact angle (obtained from a droplet formation test on a flat solid wall) be 160 degrees (i.e., almost nonwetting) or 45 degrees. The system is single-component single-relaxation (τ = Δt) following the RK EOS at a reduced temperature of TR = 0.85 resulting in a density ratio of around 12 (dimensionless liquid and vapor densities are 6.1 and 0.53, respectively). The simulation is performed in 2-D on a D2Q9 lattice topology and the on-grid bounce-back rule is used for the noslip solid nodes in the domain. As a result of these selections, the inclined walls exhibit a staircase pattern. The grid resolution is chosen such that these staircases may hardly affect the hydrodynamics. In our simulations, the total number of grid points in the domain is 216 × 412. The system is periodic in the z-direction and the global dimensionless gravitational force imposed on the system in the z-direction is 2 × 10−4. Figure 3 shows multiple snapshots of these simulations while liquid

parameters used in the LB simulations) for calculating ψδ. By doing this, one can reach for a single-component system a density ratio of 150 with the RK EOS and beyond 1000 with the CS EOS. We will now first check if the pseudopotential LB method with a nonideal EOS is capable of reproducing the correct liquid and vapor densities as a function of temperature as it is expected from the theory. To this end, we carried out an LB simulation for a single-component system to find, for the equilibrium situation, the densities at either side of a flat liquid−vapor interface in a periodic domain. For the comparison, we calculated the liquid and vapor densities through the Maxwell equal area rule for an isotherm in a P−V diagram. Figure 1 shows such a comparison for a liquid obeying the RK EOS. Figure 1 shows that the liquid density is properly described by the current pseudopotential formulation, while discrepancies in the vapor density are found for temperatures lower than, say, 0.75 TR. Similar discrepancies for a van der Waals EOS were also found by Kupershtokh et al.42 In the remainder of the paper, we only use RK EOS for TR ≥ 0.75. From Figure 1, the D



Article

Figure 2. Snapshots of phase separation at different moments in times for Carnahan and Starling EOS with TR = 0.57 and density ratio of around 1000. Colors: black is liquid phase with dimensionless density of 0.422, white is vapor with dimensionless density of 4.21 × 10−4. Note: this plot is reproducible for any other nonideal EOS.

Figure 3. Snapshots of the single-component Shan and Chen simulation of a liquid droplet motion in an inclined microchannel (45 degrees of inclination) with two expanding bends while the wall has two different wetting properties toward the liquid: the first row for almost nonwetting properties (θ0 = 160) and the second row for wetting properties (θ0 = 45). The liquid follows the RK EOS with a density ratio of around 12. t* = 1e − 4(t − t0)Δt, x* = x/L, z* = z/L.

droplets with different wall wettability are moving in this complex geometry. It is obvious from the simulations (from which these snapshots have been taken) that our scheme works properly and can handle this complex channel geometry with its sharp corners, diameter variations and change of flow direction. Also the staircase pattern of the walls does not give rise to numerical difficulties or instabilities. Although the above results on two-phase single-component flow in microchannels may be appealing and qualitatively correct, we need a more quantitative assessment for the method. For that reason, we carried out two-component gas− liquid simulations, still in 2-D, for two canonical types of flow: the rise of a single free bubble in an unconfined liquid, and the motion of a single Taylor bubble through a straight microchannel. The idea is to focus on various details of the

motion of these two types of bubbles on which literature data are readily available, in order to assess the quality of our simulation tool.

■

UNCONFINED FREE BUBBLE RISE IN 2-D The same pseudopotential LB method was used for simulating the two-phase two-component system of rise of a single bubble in an unconfined liquid in 2-D. While there is no forced flow of the liquid itself, all motions in the liquid phase are just the result of the bubble rise. We assumed the gas-like component A to behave as an ideal gas (ψA = ρA) and the liquid-like component B to obey the van der Waals EOS (ψB = {c1ρB − ((c2ρB)/(1 − c4ρB)) + c3ρB2}1/2 with c1 = 1, c2 = 1, c3 = 0.45, and c4 = 0.1).43 The coupling strengths chosen in this case were GAA = 0, GBB = −0.3, and GAB = 1.8. E



Article

Table 2. Comparison between Bubble Shapes Obtained from Our 2-D LB Simulations (Blue Is Gas), from Bhaga’s Experiment,44,45 and from Simulations (Red Is Gas) Due to Gupta et al.46

For the two components, the component density ratio ρB/ρA amounted to 2.66 while their relaxation times were τA = Δt and τB = 1.9 Δt, respectively. As a result of their interactions, the two components A and B distribute themselves over the two phases. The dimensionless densities of liquid and gas were ρliq = 2.39 and ρgas = 1.74, resulting in a phase density ratio of 1.35. The ratio νliq/νgas of the kinematic viscosities of liquid and gas follows from (τB − 0.5)/(τA − 0.5) and amounted to 2.8. By conducting the Laplace law test, the dimensionless surface tension (with respect to ρ0, Δx, Δt) of this system was found to be σ = 0.24. A simulation was initialized by positioning a stationary square bubble of certain size in a stagnant liquid phase in a periodic domain of nx × ny lattice nodes. In our two-component system, we then let the bubble equilibrate, that is, relax to a 2-D sphere, and the two components redistribute over the two phases according to the parameters specified above. When the bubble had equilibrated, the buoyancy force was turned on. To ensure no net momentum was added to the system and the mass average velocity remained constant, an external body force on the system43 Fext = g(1 − ⟨ρ⟩/ρ)

We carried out several of such simulations for various bubble sizes with the view of arriving at the various bubble shapes specified in Grace’s map45 that relates the shape of a freely rising bubble (3-D) to the Re, Eo, and Mo numbers. According to this map, further substantiated by experiments due to Bhaga and Weber,44 the bubbles can be spherical (S), an oblate ellipsoid (OE) and an oblate ellipsoidal cap (OEC) depending on the pertinent Re, Eo, and Mo numbers. Table 2 presents some of our 2-D simulation results along with real-world experimental data due to Bhaga and Weber and with numerical 2-D simulation results due to Gupta and Kumar.46 As a matter of fact, the cases a, b, and c in Table 2 show we were successful in reproducing the various bubble shapes S, OE, and OEC within the pertinent subdomains of Grace’s map. This success adds substantially to the validity and applicability of our two-component pseudopotential LB code. The fact that the phase density ratio in these free bubble rise simulations was just 1.35, does not take away this success: the code is pseudopotential two-component gas−liquid with the liquid obeying the van der Waals EOS and is perfectly capable of dealing with other EOSs and higher phase density ratios as shown in the simulations resulting in Figures 2 and 3.

■

(9)

TAYLOR BUBBLE FLOW IN A STRAIGHT MICROCHANNEL The second canonical flow with the help of which we have validated the findings of the current pseudopotential multicomponent LB code, is the steady motion of a long air bubble in a microchannel, often called a Taylor bubble. Fairbrother and Stubbs,47 Taylor48 and Bretherton49 were the first to study this type of flow that recently has attracted renewed interest in the context of microfluidics and monolithic reactors. In addition to the Reynolds number Re, the flow is usually dominated by the Capillary number, being defined as Ca = μliqUbub/σ and denoting the ratio of viscous to interfacial forces. In our validation study, we will focus on the thickness of the liquid film between bubble and wall, on pressure and velocity fields, on bubble shape deformation, and on flow pattern inside and outside of the bubble. We are interested in whether or not our code is capable of correctly reproducing these aspects as a function of Ca and Re.

was applied, where g denotes the gravitational acceleration, ρ is the local density of the mixture of the two components (so, ρ = ρA + ρB) at each lattice node, and ⟨ρ⟩ is the same mixture density but now averaged over the whole computational domain. After having turned the buoyancy force on, the bubble started rising until it reached its eventual bubble shape and its pertinent terminal velocity. At the end of the simulation, the equivalent bubble diameter of the eventual bubble was calculated from de = ((4/π)ϕnxny)1/2 where the gas fraction ϕ is given by ϕ = (ρliq − ⟨ρ⟩)/(ρliq − ρgas). In our simulations, de was always pretty close to the equivalent diameter of the square bubble with which we started the simulation. The corresponding Reynolds number Re, the Eötvos number Eo, and the Morton number Mo are defined as Eo = gde2Δρ/σ, Mo = gρliq2νliq2Δρ/σ3, Re = deUT/νliq where Δρ = ρliq − ρgas and UT is the eventual steady-state terminal velocity found in the simulation. F



Article

Details of the Taylor Bubble Simulation. Our simulation relates to a Taylor bubble moving in a 2-D channel of width D and length L = 8D. Inlet and outlet of the channel are periodic which resembles the motion of a train of bubbles in the channel. Provided that the bubble length is large enough (>1.5D) and the slug length to bubble length ratio is sufficiently large, this will mimic the motion of an infinitely large isolated bubble,50 at least as far as film thickness and pressure drop are concerned. The liquid-like component (denoted by B) is described with an exponential expression for the potential function ψB = |1 − e−20ρB| and with GBB = −1.6. Component A is an ideal gas with ψA = ρA and GAA = 0. The interaction between the two components is governed by G AB = +0.2, while the dimensionless surface tension (with respect to ρ0, Δx, Δt) is σ = 0.10349. The liquid to gas density ratio amounted to 8.7 with a dimensionless liquid density of ρliq = 2.0 . The wall (denoted by s) is fully wetted by a liquid film due to the neutral boundary condition described before. We used a single relaxation time as a result of which the kinematic viscosities of gas-like and liquid-like components were the same and the ratio of the dynamic viscosities is equal to the density ratio. Increasing the viscosity ratio by using different relaxation times for the different components may decrease the stability of the method and cause reduction of the spatial accuracy of the simulation.51 Simulating the thin liquid adjacent to the wall requires special attention, as a gas−liquid interface in LB is diffuse, that is, spreads across 3 to 5 grid lines. In our simulations, the position of gas−liquid interface was defined there where the maximum in the absolute gradient of the total density is available which might be associated with interpolations in the spatial distribution of total density. As with Taylor bubble flow at small Ca numbers (Ca < 0.005) the film thickness is smaller than 1% of the channel width, we had to use a fine grid. A halfway bounce-back rule, the wall nodes being located halfway between two grid line (see Figure 4), was imposed on each of

Our simulations were initialized by putting an isolated bubble in the middle of the channel as a rectangle which equilibrates into an elongated bubble. The rectangle length is chosen such that after equilibration the bubble length is not smaller than two times the channel width. After having reached a steady bubble shape, an external body force equivalent to the Hagen− Poiseuille condition was imposed on the system leading to the motion of the bubble within the channel. The simulation was continued until a steady motion of the bubble was attained which usually happened in less than 2 × 105 time steps. The liquid flow was always laminar with Re < 80 and Ca < 1 (in our simulations Re scales with Ca) in the absence of gravity. Film Thickness. Figure 5 presents our simulation results for dimensionless film thickness h/D versus Ca in comparison with experimental literature data as well as with computational findings from the literature. First of all, it is clear that grid resolution is important: the finer our grid, the better our simulations reproduce the experimental data due to Fairbrother and Stubbs47 and Taylor,48 which for Ca < 0.1 are described by Ca0.5. With our very fine grid, our results for Ca < 0.005 also excellently agree with Bretherton’s analytical predictionon the basis of lubrication theorythat film thickness varies with Ca2/3. Our LB simulations confirm that beyond Ca = 0.005 Bretherton’s lubrication theory does not give a satisfactory prediction any longer. For larger Ca numbers, the experimental correlation of Aussillous and Queré52 is well reproduced. Overall, our simulations agree best with the experimental data and the correlation obtained by Han et al.53 Further, our method is at least as accurate as other numerical studies in predicting film thickness. This is shown by comparing our LB simulations with those carried out by Yang et al.54 using the Shan and Chen LB technique, with the Finite Difference study of Reinelt and Saffman,55 and the Finite Element investigation due to Giavedoni and Saita.56 Our method compares favorably with the free-energy based LB simulation study for 0.01 < Ca < 1.0 and Re < 20 conducted by Kuzmin et al.57 in which the two phases had the same density while the viscosity ratio was varied up to as high as 20 by altering the relaxation times. The Bubble to Liquid Velocity Ratio. The ratio of bubble to liquid velocity is another important variable in Taylor bubble motion for which we compare our results with the correlation proposed by Bretherton49 for Ca ≪ 1 and Re ≪ 1: Ubub = 1 + 0.643(3Ca)2/3 Uliq

(10)

in which Ubub and Uliq represent the steady-state velocity of the gas−liquid interface and the steady-state liquid velocity averaged over the channel’s cross-section far away from the bubble, respectively. We plotted the above velocity ratio obtained in our simulations on both the coarse and the fine grids versus Ca as open circles in Figure 6 along with Bretherton’s correlation which is supposed to be valid for Ca < 0.003 only. Our simulation results suggest it may apply in the range 0.003 < Ca < 0.1 as well. The Pressure Field in the Liquid Film. Figure 7 presents the nondimensional pressure change along the channel in the vicinity of the wall (one grid line away) in the vicinity of a moving Taylor bubble motion at two different Ca and Re numbers. The nondimensional pressure was defined as P* = P D/(νliqρliqU̅ liq). With x* = x/D, this figure illustrates that far behind (x* < 3) and ahead of (x* > 5) the Taylor bubble, the

Figure 4. Film thickness measurement in halfway bounce back scheme.

the channel walls, resulting in a second-order accurate no-slip boundary condition. Owing to the fine grid used, the film always contained a number of grid lines (see Figure 4). In successive simulations, the number of grid points was 64, 128, and 256 in the direction normal to the wall (z-axis) and 512, 1024, and 2048 in the flow direction (x-axis); these meshes are denoted as coarse, fine, and very fine, respectively. Due to the halfway bounce-back treatment of the channel walls, the effective channel width was 62, 126, or 254. G



Article

Figure 5. Film thickness variation in variety of Ca numbers: comparison between LB simulation results and different theoretical, experimental and numerical studies;47,48,53,55−57 Closed triangles, open circles and closed circles are the results of present work with grid resolution of 256 × 2048, 128 × 1024, and 64 × 512, respectively.

Ignoring gas phase viscosity is common practice in Taylor bubble studies. According to Afkhami et al.,58 the pressure drop in the liquid film (sandwiched between the no-slip wall and the partial-slip gas−liquid interface) when not ignoring the gas phase viscosity can be calculated from: ⎞ ΔPfilm γ ΔP ⎛ with =− ⎜ 3 3⎟ L − Lslug L ⎝ γ + ε − γε ⎠ μgas ρgas 2h and ε = 1 − γ= = D μ liq ρliq

(11)

in which ΔP/L stands for the total frictional pressure drop in the channel in the absence of the bubble and h still denotes the film thickness. We did not ignore gas phase viscosity in our LB simulations. Thus, the pressure in the liquid film is not constant anymore, as opposed to what was obtained by Kreutzer et al.59 for instance. Our pressure drop in the film region (with ε = 0.84 for the solids dots, ε = 0.74 for the open symbols, and γ = 0.11) matches quite well with eq 11 . The wiggle in the simulated pressure curve close to the tail of the bubble for Ca = 0.08, Re = 8.87 is the result of undulations in the rear meniscus of the bubble (also reported by Kreutzer et al.59). At higher Ca numbers, the undulations are damped. Bubble Shape. Total pressure over a Taylor bubble will affect the bubble shape, while bubble shape may depend mainly on Ca and Re.59 Figure 8 presents bubble shapes for several cases of varying Ca and Re. The bubble shapes in this figure are

Figure 6. Bubble to liquid velocity ratio versus Ca: open circles denote our LB results; solid line represents eq 10.

Figure 7. Dimensionless pressure (total pressure minus external pressure difference imposed) along the channel in the vicinity of the wall. Solid line and filled circles: LB simulation for Ca = 0.08, Re = 8.87. Open circles: LB simulation for Ca = 0.32, Re = 36.6. Broken lines, Hagen−Poiseuille; dash-dot lines, eq 11. Figure 8. Bubble shape for different Ca and Re numbers: (1) Ca = 2.30 × 10−3, Re = 0.27; (2) Ca = 1.72 × 10−2, Re = 1.99; (3) Ca = 7.69 × 10−2, Re = 8.87; (4) Ca = 3.17 × 10−2, Re = 36.59; (5) Ca = 6.25 × 10−1, Re = 71.85.

liquid phase obeys ∂P*/∂x* = 12 for fully developed laminar liquid flow in a 2-D channel. H



Article

investigated the effect of Ca and Re on both undulations (wavelength, amplitude) and shape change by means of a Finite Element method and concluded that Re plays a dominant role with respect to the undulations while essentially Ca affects shape change. Figure 10 shows the bubble interfaces near the rear meniscus as obtained by our LB simulation (solid and broken lines) for a variety of Re and Ca and compares them with those obtained by Giavedoni and Saita for almost similar values of Ca and Re (shown by markers). Generally, there is a good agreement between our simulation results and those found by Giavedoni and Saita. Figure 10 also illustrates that there is a region with constant film thickness far away from the rear meniscus. This constant film thickness persists even when Re increases. This can be also seen more clearly from Figure 8 in which the whole bubble is depicted. Edvinsson and Irandoust61 argued that the film thickness is not constant any more as the Re number increases. Our study and that due to Giavedoni and Saita, however, suggest that a liquid film of constant thickness is produced even when the inertial forces are not negligible. We tend to the following conclusions: for very small Ca (Ca ≈ 0.001), even at strongly different Re (Re = 0.2 versus Re = 50), the effect of Re on bubble shape is minor, while film thickness mainly depends on Ca as long as We = Re·Ca ≪ 1. We do not find a strong undulation at the rear side of the bubble under these conditions. With increasing Re and Ca, the almost perfectly spherical cap at the rear side starts disappearing while both film thickness and the undulation increase. For Ca = 0.62, Re = 74, a gas−liquid interface at the vertex point of the rear end has changed from convex to concave, and the film is rather thick with virtually no undulation. The high gradient in the curvature in the latter case may create some difficulties in numerical approaches that uses the calculated curvature as an input for momentum calculations, as previously reported for instance by Giavedoni and Saita60 and by Kreutzer et al.59 Our LB method, however, does not seem to suffer from such instability issues. A more detailed analysis of the undulations in film thickness at the rear side of the Taylor bubble is presented in Figure 11 that shows how h − hmin, nondimensionalized with D, varies with Ca, h being the constant film thickness in the middle of the bubble or far away from the rear meniscus, and hmin being

found considering the iso-lines of the total density at the averaged value of the gas and liquid phase densities. For the front meniscus of the bubble, Bretherton’s theory at Ca ≪ 1 suggests a spherical cap and a transitional region, labeled AB and BC, respectively, in the inset of Figure 9 (just

Figure 9. Comparison between simulated bubble shape at the tip of the bubble (solid line) with eq 12 (dashed line) for Ca = 0.004. Inset: overall bubble shape obtained from LB simulation.

like in Bretherton’s paper); the transitional interface profile is described by z* = (x*)2 + 0.895(3Ca)2/3

(12)

in which z* = z/D, x* = x/D are nondimensional distances from point C. Figure 9 compares a bubble nose shape as obtained from LB simulations with eq 12 for Ca = 0.004 and Re ≪ 1. The rear side of the Taylor bubble exhibits a more complicated shape, involving undulations in the film thickness, as well as a gradual change (with increasing Re) from a convex to a concave shape at the vertex point of the rear interface, as also shown in Figure 8 and Figure 10. As the Bretherton’s solution to the problem is concerned, there is no unique solution for the rear gas−liquid interface profile contrary to the one in the front. In fact, a unique solution will be only obtained for a certain film thickness found from the solution of the problem in the front meniscus. Except for the limiting case of quasi-static motion of the bubble in a tube, full solution to the problem of bubble motion in the tube is only possible by using numerical methods. Giavedoni and Saita60 systematically

Figure 10. Comparison between our LB simulation results for the shape of the rear of the bubble compared to the Finite Element results of Giavedoni and Saita.60 Lines or broken lines obtained from current LB simulations (counting from the top): (1st) Ca = 0.001, Re = 0.2; (2nd) Ca = 0.01, Re = 2; (3rd) Ca = 0.04, Re = 4; (4th) Ca = 0.08, Re = 9; (5th) Ca = 0.3, Re = 38; (6th) Ca = 0.62, Re = 74. Markers obtained by Giavedoni and Saita: (circels) Ca = 0.001, Re = 50; (squares) Ca = 0.01, Re = 0; (stars) Ca = 0.1, Re = 10; (triangles) Ca = 0.25, Re = 50. I



Article

Figure 11. Amplitude of the wiggle in the rear side of the bubble obtained from LB simulation compared with the results of Giavedoni60 and Bretherton’s eq 13.

the film thickness right at the undulation point. We compare our simulation results with Bretherton’s relation: (h − hmin)/D = 0.1899Ca 2/3

(13)

as well as with simulation results due to Giavedoni and Saita.60 While, generally, our results are in line with these data, our LB simulations for various Ca and Re numbers suggest the amplitude of the undulation does not depend on Re. Velocity Field Inside and Outside the Bubble. With the view of analyzing the velocity field predictions, it is most convenient to use a frame of reference moving with the Taylor bubble. Figure 12 presents the velocity field in terms of velocity vectors and streamlines after having subtracted the bubble velocity from the velocities in the laboratory frame of reference. A few aspects of the flow outside the bubble should be discussed in some greater detail: (i) with increasing Ca, the bubble shape gradually changes from bullet to torpedo while film thickness increases giving room for a bypass flow between bubble and wall; (ii) in the liquid slug, there is a flow circulation with respect to the bubble adding to the (lateral) mixing in the liquid slugmuch in line with the findings due to Heil;62 (iii) at low Ca and Re values, some outliers are visible at the nose and at the rear side of the Taylor bubble, especially at positions where the gas−liquid interface exhibits a strong curvature. These are the spurious currents along curved interfaces typical of LB simulations. In the simulations presented, the parasitic currents are usually very small compared to the average velocity in the computational domain; they only show up at low Ca and Re cases where they are of the order 0.005Δx/Δt. We were not able to detect any effect on aspects such as film thickness and bubble curvature. By not ignoring viscosity inside the Taylor bubble, we were capable of calculating the internal flow field and doing so we found vortices inside the bubble which are indicative of a certain degree of internal mixing. In this respect, our simulations differ essentially from the usual models and simulations.48,62 The internal vortices presented here were previously seen in the study carried out by Afkhami et al.58 on a droplet moving through an immiscible liquid in a tube. While in

Figure 12. Velocity field (plotted in the frame of reference moving with the bubble) and streamlines as simulated for a variety of Ca and Re numbers.

all our simulations a large internal circulation extends over most of the bubble, at low Ca values (when interfacial tension may be dominant) additional smaller vortices are observed at the nose and at the rear side. In that respect, our simulations produce results different from those due to Yu et al.14 J



■

Article

resulted in a very versatile tool for simulating gas−liquid flows at relevant conditions. As far as the Taylor bubble simulations are concerned, our method allows gas viscosity to not be ignored. This has resulted in the prediction of more realistic pressure profiles in the film between Taylor bubble and pipe wall as well as of mixing patterns inside the Taylor bubble. Such results are indicative of the potential of the methodology we described and implemented.

CONCLUSIONS AND OUTLOOK In this study, a pseudopotential multicomponent LB method has successfully been applied to simulate various types of gas− liquid flow systems of both industrial and academic interest. The specific cases we investigated, albeit provisionally in 2-D, were dispersed two-phase flow through an inclined microchannel comprising two bends, the free rise of a single bubble in an unconfined liquid, and an individual Taylor bubble in a microchannel. To this end, we extended the original Shan and Chen formulation of the pseudopotential multicomponent LB methodbeing restricted to density ratios not exceeding, say, 10with ideas proposed by Yuan and Schaefer,2 Yu et al.14 and Kupershtokh et al.42 By exploiting potential functions ψ which reflect nonideal EOSs and by using coupling strengths G which relate to the mutual interactions between components, we are capable of mimicking systems with a density ratio beyond, say, 10, typical of gas−liquid flows. With the Redlich−Kwong EOS it is possible to attain a density ratio of about 150, while for the Carnahan−Starling EOS it can be as high as about 1000 in a single-component system. The extension of Shan and Chen’s LB method for reproducing the wetting/nonwetting behavior of solid walls for simulations of wall-bounded gas−liquid flows was also incorporated in the current method. Combining all these aspects has resulted into a tool suited for simulating gas− liquid flows at industrially relevant conditions such as high(er) values for the phase density ratio and the presence of tube walls. The capabilities of the current method were first of all demonstrated by a drop formation study starting from a random distribution of liquid and vapor. This study (see Figure 2) is similar to that presented before by Shan and Chen1 but now for the higher density ratio of 1000 in a single-component fluid which obeys the Carnahan−Starling EOS. From the simulations of the dispersed gas−liquid flow through the inclined channel we conclude that our code is capable of reproducing physically meaningful results without suffering from numerical instabilities. The selection of proper values for the potential functions, the various coupling strengths, and the relaxation time is crucial in getting smoothly running simulations. The simulations of the free bubble rise and the Taylor bubble motion show that the method produces qualitatively and quantitatively correct results. Not all simulations reported in this paper were carried out for high values of the phase density ratio. Yet, our experience with our code gives us reason to believe that our method is widely applicable as it incorporates essential physical characteristics such as a nonideal equation of state and the interaction with a solid wall. This is realized by means of the use of potential functions ψ and coupling strengths G. This way of mimicking the effect of the mesoscopic physical interactions may be more successful than the phenomenological, more superficial description of macroscopic phenomena of phase interaction in terms of surface tension and slip velocity. A very characteristic feature of this type of two-phase LB simulations is in the automatic generation of the numerous phase interfaces in many types of gas−liquid flows. This is due to calculating interparticle forces made up by interacting potential functions and coupling strengths. Real-life fluids obeying nonideal equations of state can be mimicked by defining appropriate potential functions. All these characteristics have been combined and implemented into a single code. This has

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors are very much indebted to Professor Sankaran Sundaresan (Princeton University) and Dr Jurriaan Gillissen (Delft University of Technology) for many challenging and interesting discussions and for proof reading this paper. The financial support of the Dutch Technology Foundation STW, the Applied Science Division of The Netherlands Organization for Scientific Research (NWO) is highly appreciated. M.R.K. would like to also thank the Islamic Development Bank for the awarded scholarship through the Merit Scholarship Programme for High Technology.

■

REFERENCES

(1) Shan, X.; Chen, H. Lattice Boltzmann model for simulating flows with multiple phases and components. Phys. Rev. E 1993, 47 (3), 1815−1819. (2) Yuan, P.; Schaefer, L., Equations of state in a lattice Boltzmann model. Phys. Fluids 2006, 18, (4). (3) Raabe, D. Overview of the lattice Boltzmann method for nanoand microscale fluid dynamics in materials science and engineering. Modell. Simul. Mater. Sci. Eng. 2004, 12 (6), R13−R46. (4) Eggels, J. G. M. Direct and large-eddy simulation of turbulent fluid flow using the lattice-Boltzmann scheme. Int. J. Heat Fluid Flow 1996, 17 (3), 307−323. (5) Derksen, J.; Van Den Akker, H. E. A. Large eddy simulations on the flow driven by a Rushton turbine. AIChE J. 1999, 45 (2), 209−221. (6) Derksen, J. J.; Kontomaris, K.; McLaughlin, J. B.; Van den Akker, H. E. A. Large-eddy simulations of single-phase flow dynamics and mixing in an industrial crystallizer. Chem. Eng. Res. Des. 2007, 85 (2 A), 169−179. (7) Van Vliet, E.; Derksen, J. J.; Van Den Akker, H. E. A. Turbulent mixing in a tubular reactor: Assessment of an FDF/LES approach. AIChE J. 2005, 51 (3), 725−739. (8) Van den Akker, H. E. A., The Details of Turbulent Mixing Process and their Simulation. In Advances in Chemical Engineering; Elsevier, 2006; Vol. 31, pp 151−229. (9) Vreman, B.; Geurts, B. J.; Deen, N. G.; Kuipers, J. A. M.; Kuerten, J. G. M. Two- and four-way coupled euler-lagrangian large-eddy simulation of turbulent particle-laden channel flow. Flow, Turbulence and Combustion 2009, 82 (1), 47−71. (10) Beetstra, R.; van der Hoef, M. A.; Kuipers, J. A. M. A latticeBoltzmann simulation study of the drag coefficient of clusters of spheres. Computers and Fluids 2006, 35 (8−9), 966−970. (11) Holloway, W.; Yin, X.; Sundaresan, S. Fluid-particle drag in inertial polydisperse gas-solid suspensions. AIChE J. 2010, 56 (8), 1995−2004. (12) Yu, Z.; Fan, L. S. Lattice Boltzmann method for simulating particle-fluid interactions. Particuology 2010, 8 (6), 539−543. K



Article

(37) Wolf-Gladrow, D. A. Lattice-Gas Cellular Automata and Lattice Boltzmann Models: An Introduction; Springer: New York, 2000. (38) Martys, N. S.; Chen, H. Simulation of multicomponent fluids in complex three-dimensional geometries by the lattice Boltzmann method. Phys. Rev. E 1996, 53 (1 SUPPL. B), 743−750. (39) Huang, W.; Li, Y.; Liu, Q. Application of the lattice Boltzmann method to electrohydrodynamics: Deformation and instability of liquid drops in electrostatic fields. Chin. Sci. Bull. 2007, 52 (24), 3319− 3324. (40) Zhang, J. Lattice Boltzmann method for microfluidics: Models and applications. Microfluid. Nanofluid. 2011, 10 (1), 1−28. (41) Kamali, M. R.; Gillissen, J. J. J.; Sundaresan, S.; Van den Akker, H. E. A. Contact line motion without slip in lattice Boltzmann simulations. Chem. Eng. Sci. 2011, 66 (14), 3452−3458. (42) Kupershtokh, A. L.; Medvedev, D. A.; Karpov, D. I. On equations of state in a lattice Boltzmann method. Comput. Math. Appl. 2009, 58 (5), 965−974. (43) 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, 61−96. (44) Bhaga, D.; Weber, M. E. Bubbles in Viscous Liquids: Shapes, Wakes and Velocities. J. Fluid Mech. 1981, 105, 61−85. (45) Grace, J. R. Shapes and Velocities of Bubbles Rising in Infinite Liquids. Chem. Eng. Res. Des. 1973, 51a, 116−120. (46) Gupta, A.; Kumar, R. Lattice Boltzmann simulation to study multiple bubble dynamics. Int. J. Heat Mass Transfer 2008, 51 (21− 22), 5192−5203. (47) Fairbrother, F.; Stubbs, A. E. Studies in electro-endosmosis. Part VI. The ″Bubble-tube″ method of measurement. J. Chem. Soc. (Resumed) 1935, 527−529. (48) Taylor, G. I. Deposition of a viscous fluid on the wall of a tube. J. Fluid Mech. 1961, 10 (2), 161−165. (49) Bretherton, F. P. The motion of long bubbles in tubes. J. Fluid Mech. 1961, 10, 166−188. (50) Ratulowski, J.; Chang, H. C. Transport of gas bubbles in capillaries. Phys. Fluids A 1989, 1 (10), 1642−1655. (51) Kuzmin, A.; Januszewski, M.; Eskin, D.; Mostowfi, F.; Derksen, J. J. Three-dimensional binary-liquid lattice Boltzmann simulation of microchannels with rectangular cross sections. Chem. Eng. J. 2011, 178, 306−316. (52) Aussillous, P.; Quere, D. Quick deposition of a fluid on the wall of a tube. Phys. Fluids 2000, 12 (10), 2367−2371. (53) Han, Y.; Shikazono, N. Measurement of liquid film thickness in micro square channel. Int. J. Multiphase Flow 2009, 35 (10), 896−903. (54) Yang, Z. L.; Palm, B.; Sehgal, B. R. Numerical simulation of bubbly two-phase flow in a narrow channel. Int. J. Heat Mass Transfer 2002, 45 (3), 631−639. (55) Reinelt, D. A.; Saffman, P. G. The penetration of a finger into a viscous fluid in a channel and tube. SIAM J. Sci. Stat. Comput. 1985, 6 (3), 542−561. (56) Giavedoni, M. D.; Saita, F. A. The axisymmetric and plane cases of a gas phase steadily displacing a Newtonian liquidA simultaneous solution of the governing equations. Phys. Fluids 1997, 9 (8), 2420− 2428. (57) Kuzmin, A.; Januszewski, M.; Eskin, D.; Mostowfi, F.; Derksen, J. J. Simulations of gravity-driven flow of binary liquids in microchannels. Chem. Eng. J. 2011, 171 (2), 646−654. (58) Afkhami, S.; Leshansky, A. M.; Renardy, Y., Numerical investigation of elongated drops in a microfluidic T-junction. Phys. Fluids 2011, 23, (2). (59) Kreutzer, M. T.; Kapteijn, F.; Moulijn, J. A.; Kleijn, C. R.; Heiszwolf, J. J. Inertial and interfacial effects on pressure drop of Taylor flow in capillaries. AIChE J. 2005, 51 (9), 2428−2440. (60) Giavedoni, M. D.; Saita, F. A. The rear meniscus of a long bubble steadily displacing a Newtonian liquid in a capillary tube. Phys. Fluids 1999, 11 (4), 786−794. (61) Edvinsson, R. K.; Irandoust, S. Finite-element analysis of taylor flow. AIChE J. 1996, 42 (7), 1815−1823.

(13) 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. (14) Yu, Z.; Hemminger, O.; Fan, L. S. Experiment and lattice Boltzmann simulation of two-phase gas-liquid flows in microchannels. Chem. Eng. Sci. 2007, 62 (24), 7172−7183. (15) Yin, X.; Koch, D. L. The lift force on a bubble in a sheared suspension in a slightly inclined channel. J. Fluid Mech. 2008, 615, 27− 51. (16) Yin, X.; Koch, D. L., Lattice-Boltzmann simulation of finite Reynolds number buoyancy-driven bubbly flows in periodic and wallbounded domains. Phys. Fluids 2008, 20, (10). (17) Gillissen, J. J. J.; Sundaresan, S.; Van Den Akker, H. E. A. A lattice Boltzmann study on the drag force in bubble swarms. J. Fluid Mech. 2011, 679, 101−121. (18) Kang, Q.; Zhang, D.; Chen, S. Displacement of a twodimensional immiscible droplet in a channel. Phys. Fluids 2002, 14 (9), 3203−3214. (19) Briant, A. J.; Wagner, A. J.; Yeomans, J. M. Lattice Boltzmann simulations of contact line motion. I. Liquid−gas systems. Phys. Rev. E 2004, 69 (3 1), 031602−1−031602−14. (20) Zhang, J.; Kwok, D. Y. Contact line and contact angle dynamics in superhydrophobic channels. Langmuir 2006, 22 (11), 4998−5004. (21) Hyväluoma, J.; Koponen, A.; Raiskinmäki, P.; Timonen, J. Droplets on inclined rough surfaces. Eur. Phys. J. E 2007, 23 (3), 289− 293. (22) Wolf, F. G.; Dos Santos, L. O. E.; Philippi, P. C., Modeling and simulation of the fluid-solid interaction in wetting. J. Stat. Mech.: Theory Exp. 2009, 2009, (6). (23) Liu, H.; Zhang, Y., Droplet formation in a T-shaped microfluidic junction. J. Appl. Phys. 2009, 106, (3). (24) Verberg, R.; Pooley, C. M.; Yeomans, J. M.; Balazs, A. C. Pattern formation in binary fluids confined between rough, chemically heterogeneous surfaces. Phys. Rev. Lett. 2004, 93 (18), 184501−1− 184501−4. (25) Horbach, J.; Succi, S., Lattice Boltzmann versus molecular dynamics simulation of nanoscale hydrodynamic flows. Phys. Rev. Lett. 2006, 96, (22). (26) Sbragaglia, M.; Benzi, R.; Biferale, L.; Succi, S.; Toschi, F., Surface roughness-hydrophobicity coupling in microchannel and nanochannel flows. Phys. Rev. Lett. 2006, 97, (20). (27) Harting, J.; Kunert, C.; Hyväluoma, J. Lattice Boltzmann simulations in microfluidics: Probing the no-slip boundary condition in hydrophobic, rough, and surface nanobubble laden microchannels. Microfluid. Nanofluid. 2010, 8 (1), 1−10. (28) Swift, M. R.; Orlandini, E.; Osborn, W. R.; Yeomans, J. M. Lattice Boltzmann simulations of liquid-gas and binary fluid systems. Phys. Rev. E 1996, 54 (5), 5041−5052. (29) Cahn, J. W.; Hilliard, J. E. Free energy of a nonuniform system. I. Interfacial free energy. J. Chem. Phys. 1958, 28 (2), 258−267. (30) Shan, X.; Chen, H. Simulation of nonideal gases and liquid-gas phase transitions by the lattice Boltzmann equation. Phys. Rev. E 1994, 49 (4), 2941−2948. (31) Shan, X.; Doolen, G. Multicomponent lattice-Boltzmann model with interparticle interaction. J. Stat. Phys. 1995, 81 (1−2), 379−393. (32) Vervloet, D.; Kamali, M. R.; Gillissen, J. J. J.; Nijenhuis, J.; van den Akker, H. E. A.; Kapteijn, F.; van Ommen, J. R. Intensification of co-current gas-liquid reactors using structured catalytic packings: A multiscale approach. Catal. Today 2009, 147 (SUPPL), S138−S143. (33) Deen, W. M. Analysis of Transport Phenomena; Oxford University Press: New York, 1998. (34) Bird, R. B.; Stewart, W. E.; Lightfoot, E. N. Transport Phenomena; J. Wiley: New York, 2007. (35) Succi, S. The Lattice Boltzmann Equation for Fluid Dynamics and Beyond; Clarendon Press: Oxford, 2001. (36) Bhatnagar, P. L.; Gross, E. P.; Krook, M. A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems. Phys. Rev. 1954, 94 (3), 511−525. L



Article

(62) Heil, M. Finite Reynolds number effects in the Bretherton problem. Phys. Fluids 2001, 13 (9), 2517−2521.

M


Simulating GasâLiquid Flows by Means of a Pseudopotential Lattice

Recommend Documents

Simulating GasâLiquid Flows by Means of a Pseudopotential Lattice