界面

A lattice Boltzmann method for incompressibletwo-phase flows with large density differencesT.Inamuro *,T.Ogata,S.Tajima,N.KonishiDepartment of Chemical Engineering,Graduate School of Engineering,Kyoto University,Katsura Campus,Kyoto 615-8510,JapanReceived 13June 2003;received in revised form 20December 2003;accepted 12January 2004Available online 21February 2004AbstractA lattice Boltzmann method for two-phase immiscible fluids with large density differences is proposed.The difficulty in the treatment of large density difference is resolved by using the projection method.The method can be applied to simulate two-phase fluid flows with the density ratio up to 1000.To show the validity of the method,we apply the method to the simulations of capillary waves,binary droplet collisions,and bubble flows.In capillary waves,the an-gular frequencies of the oscillation of an ellipsoidal droplet are obtained in good agreement with theoretical ones.In the simulations of binary droplet collisions,coalescence collision and two different types of separating collisions,namely reflexive and stretching separations,can be simulated,and the boundaries of the three types of collisions are in good agreement with an available theoretical prediction.In the bubble flows,the effect of mobility on the coalescence of two rising bubbles is investigated.The behavior of many bubbles in a square duct is also simulated.Ó2004Elsevier Inc.All rights reserved.AMS:65M99;76T10Keywords:Lattice Boltzmann method;Two-phase flow;Capillary wave;Binary droplet collision;Bubble flow1.IntroductionRecently,the lattice Boltzmann method (LBM)has been developed into an alternative and promising numerical scheme for simulating multicomponent fluid flows.Gunstensen et al.[1]developed a multi-component LBM based on the two-component lattice gas model.Shan and Chen [2]proposed an LBM model with mean-field interactions for multiphase and multicomponent fluid flows.Swift et al.[3]devel-oped an LBM model for multiphase and multicomponent fluid flows using the free-energy approach.He et al.[4]proposed a new lattice Boltzmann multiphase model using the idea of level-set for multiphase flow.*Corresponding author.Tel.:+81-75-753-5791;fax:+81-75-761-4947.E-mail address:inamuro@cheme.kyoto-u.ac.jp (T.Inamuro)/locate/jcpJournal of Computational Physics 198(2004)628–6440021-9991/$-see front matter Ó2004Elsevier Inc.All rights reserved.doi:10.1016/j.jcp.2004.01.019T.Inamuro et al./Journal of Computational Physics198(2004)628–644629 Inamuro et al.[5,6]also proposed an LBM for multicomponent immisciblefluids with the same density. Ginzburg and Steiner[7]have proposed a generalized LBM model to simulate free-surfaceflows by modifying the above model by Gunstensen et al.T€o lke et al.[8]also have proposed an LBM for immiscible binaryfluids with variable viscosities and density ratios based on the model by Gunstensen et al.and applied the method to the simulation of binaryflows through porous media.Sankaranarayanan and Sundaresan[9]and Sankaranarayanan et al.[10,11]have developed an implicit LBM for multicomponent flows based on the above Shan and Chen model and used the method for the analysis of drag,lift,and virtual mass forces in bubbly suspensions.The above LBM models have great advantages over conven-tional methods for multiphaseflows.They do not track interfaces,but can maintain sharp interfaces without any artificial treatments.Also,the LBM is accurate for the mass conservation of each component fluid.While most of the above LBM models for multiphase and multicomponentfluidflows are based on heuristic ideas with no direct connection to kinetic theory,Luo[12,13]and Luo and Girimaji[14,15]have rigorously derived the LBM models for multiphasefluids from the Enskog equation and for multicom-ponentfluids from the corresponding kinetic equations and provided a unified framework to treat the LBM models for multiphase and multicomponentfluids.On the other hand,diffuse-interface methods are used for many applications to interfacial phenomena [16–19].Each of the above LBM models may be regarded as one of the diffuse-interface methods.In the LBM models,the pressure tensor is introduced through the collision term.Since the LBM is classified as a mesoscopic approach to the simulation offluid dynamics,the combination of LBM and diffuse-interface methods is suitable for the simulation of two-phasefluidflows.Although the LBM is a promising method for multicomponentfluidflows,one of disadvantages is that all above schemes are limited to small density ratios less than10due to the numerical instability in the interface with large density ually the density ratio of liquid–gas systems is larger than100, e.g.,the density ratio of water to air is about1000:1.Thus,the development of an LBM for two-phase fluids with large density ratios is required.Teng et al.[20]used the total variation diminishing with artificial compression(TVD/AC)scheme to the lattice Boltzmann multiphase model in order to stabilize the computation for large density ratios up to100,but the method has been applied only to the problems with infinitesimalflows such as phase separation.The numerical instability in the interface with large density ratios is mainly caused by spurious velocities near the interface which become larger as the density ratio increases and do not satisfy the continuity equation.In the simulation offlows with large density differences,a key to the problem is to assure the continuity equation near the interface at each time step.The aim of the present paper is to propose a heuristic LBM for two-phasefluids with large density differences.The difficulty in the treatment of large density differences is resolved by using the projection method[21].In the projection method the continuity equation in the interfacial region is satisfied at every time step.Two particle velocity distribution functions are used.One is used for the calculation of an order parameter which distinguishes two phases,and the other is used for the calculation of a predicted velocity of the two-phasefluid without a pressure gradient.The current velocity satisfying the continuity equation can be obtained by using the relation between the velocity and the pressure correction which is determined by solving the Poisson equation.In order to show the validity of the method,we apply the method to the simulations of capillary waves,binary droplet collisions,and bubbleflows.The paper is organized as follows.In Section2we propose an LBM for two-phasefluids with large density differences and derive governing equations for macroscopic variables by using the asymptotic theory.In Section3we present numerical examples,capillary waves,binary droplet collisions,and bubble flows.In capillary waves,the oscillation of an ellipsoidal droplet is calculated.In the simulations of binary droplet collisions,coalescence collision and two different types of separating collisions,namely reflexive and stretching separations,are simulated.In the bubbleflows,the effect of mobility on the coalescence of two bubbles is investigated.Concluding remarks are given in Section4.2.Numerical method2.1.Two-phase lattice Boltzmann methodNon-dimensional variables defined in Appendix A are used as in[23].In the LBM,a modeledfluid, composed of identical particles whose velocities are restricted to afinite set of N vectors c i(i¼1;2;...;N), is considered.Thefifteen-velocity model(N¼15)is used in the present paper.The velocity vectors of this model are given byc1;c2;c3;c4;c5;c6;c7;c8;c9;c10;c11;c12;c13;c14;c15½¼0100À1001À111À11À1À100100À1011À11À1À11À1000100À1111À1À1À1À11264375:ð1ÞIn the present paper we describe an athermal model for two-phasefluid.Two particle velocity distribution functions,f i and g i,are used.The function f i is used for the calculation of an order parameter which distinguishes two phases,and the function g i is used for the calculation of a predicted velocity of the two-phasefluid without a pressure gradient.Although we can use any of LBM multiphase models,we use the model by Swift et al.[3]in the present work.The evolution of the particle distribution functions f iðx;tÞand g iðx;tÞwith velocity c i at the point x and at time t is computed by the following equations:f iðxþc i D x;tþD tÞ¼f iðx;tÞÀ1f ½f iðx;tÞÀf ciðx;tÞ ;ð2Þg iðxþc i D x;tþD tÞ¼g iðx;tÞÀ1s gg iðx;tÞÂÀg ciðx;tÞÃþ3E i c i a1qoo x blo u bo x a&þo u ao x b'!D x;ð3Þwhere f ci and g ciare functions of Chapman–Enskog type in which the variables x and t enter only throughmacroscopic variables and/or their derivatives[22],s f and s g are dimensionless single relaxation times,D x is a spacing of the cubic lattice,D t is a time step during which the particles travel the lattice spacing,and the other variables,q,l,and u,and constants E i are defined below.The last term in the right-hand side of Eq.(3)represents the viscous stress tensor divided by the density of the two-phasefluid.In Eqs.(2)and(3),thefunction f ci and g ciare determined so as macroscopic variables satisfy desired equations as shown below.Inthis sense,the present method is a heuristic LBM for two-phasefluids.The order parameter/distinguishing two phases and the predicted velocity uÃof the multicomponent fluids are defined in terms of the two particle velocity distribution functions as follows: /¼X15i¼1f i;ð4ÞuüX15i¼1c i g i:ð5ÞThe functions f ci and g ciin Eqs.(2)and(3)are given byf c i ¼H i/þF i p0"Àj f/o2/o x2aÀj f6o/o x a2#þ3E i/c i a u aþE i j f G abð/Þc i a c i b;ð6Þ630T.Inamuro et al./Journal of Computational Physics198(2004)628–644g c i ¼E i1þ3c i a u aÀ32u a u aþ92c i a c i b u a u bþ32s gÀ12D xo u bo x aþo u ao x bc i a c i b!þE ij gqG abðqÞc i a c i bÀ23F ij gqo qo x a2;ð7ÞwhereE1¼2=9;E2¼E3¼E4¼ÁÁÁ¼E7¼1=9;E8¼E9¼E10¼ÁÁÁ¼E15¼1=72;H1¼1;H2¼H3¼H4¼ÁÁÁ¼H15¼0;F1¼À7=3;F i¼3E iði¼2;3;4; (15)ð8ÞandG abð/Þ¼92o/o x ao/o x bÀ32o/o x co/o x cd ab;ð9Þwith a;b;c¼x;y;z(subscripts a;b,and c represent Cartesian coordinates and the summation convention is used).In the above equations,j f is a constant parameter determining the width of the interface,j g is a constant parameter determining the strength of the surface tension,and d ab is the Kronecker delta.In Eq.(6),p0is given byp0¼/o wo/Àw¼/T11Àb/Àa/2ð10Þwithwð/;TÞ¼/T ln/1Àb/Àa/2;ð11Þwhere w is the bulk free-energy density,and a,b,and T are free parameters determining the maximumand minimum values of the order parameter/.It is noted that f ci is the same as that of the model bySwift et al.[3]except that Eq.(6)includes no quadratic term of theflow velocity u a.The last two terms of Eq.(7)represent the effect of interface tension as described below.The term including s g represents the negative counterpart of the viscous stress tensor appearing in the original LBM,and this term is needed in order to vanish the original viscous stress tensor.The followingfinite-difference approxi-mations are used to calculate thefirst and second derivatives(o/=o x a,o u b=o x a,o q=o x a,and o2/=o x2a )inEqs.(6),(7)and(9):o k o x a %110D xX15i¼2c i a kðxþc i D xÞ;ð12Þo2k o x2a %15D xX15i¼2kðx"þc i D xÞÀ14kðxÞ#:ð13ÞThe density in the interface is obtained by using the cut-offvalues of the order parameter,/ÃL and/ÃG,forthe liquid and gas phases with the following relation:T.Inamuro et al./Journal of Computational Physics198(2004)628–644631q ¼q G ;/</ÃG ;D q2sin /À/ÃD /Ãp þ1h i þq G ;/ÃG 6/6/ÃL ;q L ;/>/ÃL ;8>><>>:ð14Þwhere q G and q L are the density of gas and liquid phases,respectively,D q ¼q L Àq G ,D /ü/ÃL À/ÃG ,and üð/ÃL þ/ÃG Þ=2.It is essential to introduce the cut-offvalues of the order parameter in the present method.The maximum and minimum values of /are changed a little during calculations,and then the change of the density of each phase is magnified in proportion to the density ratio of the two-phase fluid.In order to keep the density of each phase constant,we introduce the cut-offvalues in the definition of the density of the two-phase fluid.The viscosity l in the interface is obtained byl ¼q Àq G q L Àq G ðl L Àl G Þþl G ;ð15Þwhere l G and l L are the viscosity of gas and liquid phases,respectively.The interfacial tension r is obtained byr ¼j g Z 1À1o q o n 2d n ð16Þwith n being the coordinate normal to the interface [16,26].It is noted that the flux J a of the order parameter /is given byJ a ¼Àh M /o l c o x a ¼Àh M o P ab o x b ;ð17Þwhere h M is the mobility obtained below,and l c and P ab are the chemical potential and the pressure tensor of the fluid used for the calculation of the order parameter /which are given byl c ¼o w o /Àj f o 2/o x 2a;ð18ÞandP ab ¼p 0"Àj f /o 2/o x 2c Àj f 2o /o x c 2#d ab þj f o /o x a o /o x b :ð19Þ2.2.Pressure correctionSince u Ãis not divergence free (r Áu ü0),the correction of u Ãis required.The current velocity u which satisfies the continuity equation (r Áu ¼0)can be obtained by using the following equations:Sh u Àu ÃD t¼Àr p q ;ð20Þr Ár pq ¼Sh r Áu ÃD t ;ð21Þ632T.Inamuro et al./Journal of Computational Physics 198(2004)628–644where Sh ¼U =c is the Strouhal number and p is the pressure of the two-phase fluid.It is noted that D t ¼Sh D x .The Poisson equation (21)can be solved by various methods.In the present paper,we solve Eq.(21)in the framework of ly,the following evolution equation of the velocity distribution function h i is used for the calculation of the pressure p :h n þ1i ðx þc i D x Þ¼h n i ðx ÞÀ1s h ½h n i ðx ÞÀE i p n ðx Þ À13E i o u Ãa o x a D x ;ð22Þwhere n is the number of iterations and the relaxation time s h is given bys h ¼1q þ12:ð23ÞThe pressure is obtained byp ¼X15i ¼1h i :ð24Þ2.3.Algorithm of computationWe now summarize the algorithm of computation.Step ing Eqs.(2)and (3),compute f i ðx ;t þD t Þand g i ðx ;t þD t Þ,and then compute /ðx ;t þD t Þandu Ãðx ;t þD t Þwith Eqs.(4)and (5).Also,q ðx ;t þD t Þis calculated with Eq.(14).Step ing Eq.(22)with Eqs.(23)and (24),compute p ðx ;t þD t Þ.The iteration is repeated untilj p n þ1Àp n j =q <e is satisfied in the whole domain.Step pute u ðx ;t þD t Þusing Eq.(20).Step 4.Advance one time step and return to Step 1.It is found in preliminary calculations that using the present method we can simulate multiphase flows with the density ratio up to 1000.erning equations for macroscopic variablesApplying the asymptotic theory [23–25]to Eqs.(2),(3)and (22),we find that the asymptotic expansions of the macroscopic variables with respect to a small parameter k ¼O ðD x Þ,/0¼/ð0Þþk /ð1Þ,q 0¼q ð0Þþk q ð1Þ,u 0a ¼ku ð1Þa þk 2u ð2Þa ,p 0¼k 2p ð2Þþk 3p ð3Þ,and P 0ab ¼P ð0Þab þkP ð1Þab satisfy Sh o /0o t þu 0a o /0o x a ¼h M o 2P 0ab o x a o x b ;ð25Þo u 0a o x a¼0;ð26ÞSh o u 0ao t þu 0b o u 0a o x b ¼À1s g q 0o p 0o x a þ1q 0o o x b l o u 0b o x a "þo u 0a o x b !#þo o x b j g q 0o q 0o x a o q 0o x b Ào q 0o x c o q 0o x c d ab !;ð27Þwhere h M is the mobility given byh M ¼s f À12 D x :ð28ÞT.Inamuro et al./Journal of Computational Physics 198(2004)628–644633Eq.(25)is the phase-field advection–diffusion equation (the Cahn–Hilliard equation plus advection)for /.Eqs.(26)and (27)are the continuity equation and the Navier–Stokes equations for incompressible two-phase fluid with an interfacial tension,if we choose s g ¼1.The last term of the right-hand side in Eq.(27)represents the interfacial tension given by Eq.(15).It is pointed out that large mobilities overly damp flows [19].The mobility given by Eq.(28)is of O ðD x Þ,unless s f is close to 1=2,but as s f approaches 1=2,Eq.(2)becomes unstable due to numerical instability.An idea of decreasing the mobility is to add the term of E i C ðo P ab =o x b Þc i a D x to the function f c i given by Eq.(6),where C is a constant.The mobility h M in this case is given byh M ¼s f À12À13C D x :ð29ÞBy choosing a proper value of C ,we can decrease the mobility even with s f ¼1.3.Numerical examplesTo show the validity of the present method,we apply the method to the simulations of capillary waves,binary droplet collisions,and bubble flows.3.1.Capillary waveThe oscillation of a spherical liquid droplet in a gas phase under the action of capillary force is con-sidered.Initially,an ellipsoidal droplet is placed at the center of a cubic domain.The surface of the droplet is given by r ¼R þ0:1R cos 2h ,where r is the distance from the center,R is the radius of the droplet,and h is the polar angle.The deviation corresponds to the second mode of oscillations.The density ratio of the liquid to the gas is q L =q G ¼50(q L ¼50,q G ¼1).The viscosities of the droplet and the gas are l L ¼8Â10À3D x and l G ¼1:6Â10À4D x ,respectively.The periodic boundary condition is used on all the sides of the domain.The quarter of the domain is calculated using the symmetry with x and y axes.The quarter domain is divided into a 40Â40Â80cubic lattice.The parameters in Eq.(10)are freely chosen,and for this example their values are a ¼1,b ¼6:7,and T ¼3:5Â10À2;it follows that the maximum and minimum values of the order parameter are /max ¼9:714Â10À2and /min ¼1:134Â10À2.The cut-offvalues of the order parameter are /ÃL ¼9:2Â10À2and /ÃG ¼1:5Â10À2.The other parameters are fixed atj f ¼0:5ðD x Þ2,s f ¼1,s g ¼1,C ¼0,and e ¼10À5.We change the radius of the droplet R and the interfacial tension r given by Eq.(16).Table 1Comparison of calculated angular frequencies x D x of the oscillation of the droplet with theoretical ones Rr r Theory Present Error (%)15D x14:43D x 0:2038D x 3.294Â10À3 3.162Â10À3)4.020D x19:25D x 0:05107D x 1.070Â10À3 1.034Â10À3)3.420D x19:25D x 0:2012D x 2.124Â10À3 2.068Â10À3)2.620D x19:25D x 0:7797D x 4.182Â10À3 3.957Â10À3)5.425D x 24:06D x 0:1980D x 1.508Â10À3 1.513Â10À3+0.3The theoretical angular frequencies are obtained by x D x ¼ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi8r =q L r 3p ,where r is the averaged radius ofthe initial ellipsoidal droplet [27].634T.Inamuro et al./Journal of Computational Physics 198(2004)628–644The angular frequency x of the oscillation of the droplet is obtained from the time variation of the radius of the droplet.The calculated angular frequencies x for various R and r are shown in Table 1in comparison with the theoretical ones [27].It is seen from Table 1that the calculated angular frequencies are in good agreement with the theoretical ones within the errors of 6%.As for spurious velocities induced in the interface,the maximum velocity in the gas phase is 2.1Â10À2for a stationary spherical droplet with R ¼20D x and r ¼0:2012D x [j g ¼2:5Â10À4ðD x Þ2].Since the width of the interface is thin ($3D x )and the interfacial tension is large in this case,the spurious velocities are relatively large.But they do not affect the behavior of the oscillation of the droplet,because the inertia of the gas phase is very small and the oscillation of the droplet is determined only through the ratio between the inertia of the droplet and the interfacial tension.Also,we find that the spurious velocity becomes smaller as the interfacial tension r decreases and as the lattice points in the interface increase.In numerical examples of bubble flows below,spurious velocities are much smaller than the above case (themaximumputational domain and binary dropletcollision.Fig.2.Time evolution of droplet shape for We ¼20:2and B ¼0(t ütV =D ).T.Inamuro et al./Journal of Computational Physics 198(2004)628–644635spurious velocity <5Â10À4),since much smaller interfacial tensions are used and slightly more lattice points are distributed in the interface.3.2.Droplet collisionTwo liquid droplets with the same diameter D are placed 2D apart in a gas phase,and they collide with the relative velocity V (see Fig.1).The density ratio of the liquid to the gas is q L =q G ¼50(q L ¼50,q G ¼1).The viscosities of the droplet and the gas are l L ¼8Â10À2D x and l G ¼1:6Â10À3D x ,respectively.The dimensionless parameters for binary droplet collisions are the Weber number We ¼q L DV 2=r ,the Reynolds number Re ¼q L DV =l L ,and the impact parameter B ¼X =D ,where X is the distance from the center of one droplet to the relative velocity vector placed on the center of the other droplet (see Fig.1).The periodic boundary condition is used on all the sides of the domain.The half of the domain is calculated using the symmetry with y ¼L y =2.The half domain is divided into a 192(or 128)Â48Â96cubic lattice.The pa-rameters in Eq.(10)are freely chosen,and for this example their values are a ¼1,b ¼6:7,and T ¼3:5Â10À2;it follows that the maximum and minimum values of the order parameter are /max ¼9:714Â10À2and /min ¼1:134Â10À2.The cut-offvalues of the order parameter are /ÃL ¼9:2Â10À2and /ÃG ¼1:5Â10À2.The other parameters are fixed at V ¼0:1,j f ¼0:5ðD x Þ2,s f ¼1,s g ¼1,C ¼0,e ¼10À5and D ¼32D x ,and j g is changed in the range of 20<We <80.The Reynolds number is fixed at Re ¼2000.Fig.3.Time evolution of droplet shape for We ¼20:2and B ¼0(t ütV =D ).636T.Inamuro et al./Journal of Computational Physics 198(2004)628–644Fig.4.Time evolution of droplet shape for We ¼39:7and B ¼0(t ütV =D).Fig.5.Time evolution of droplet shape for We ¼79:7and B ¼0:813(t ütV =D ).T.Inamuro et al./Journal of Computational Physics 198(2004)628–644637Fig.2shows the calculated results of time evolution of droplet shape for We ¼20:2and B ¼0.After two droplets collide head-on,they form a disk-like droplet.Owing to the large curvature at the cir-cumference of the disk-like droplet,there is a pressure difference between its inner and outer regions caused by surface tension.Thus,the disk contracts radially inward and pushes the liquid outward from its center forming a long cylinder with rounded ends.Then the cylinder oscillates until a spherical droplet is formed.This type of collision is called ‘‘coalescence collision’’.The fluid velocity fields at y ¼L y =2are shown in Fig.3.The complicated gas flows as well as liquid flows inside the droplets are clearly found.Fig.4shows the calculated results for We ¼39:7and B ¼0.The time evolution of droplet shape is similar to the previous case up to the formation of a long cylinder with rounded ends.In this case,however,the cylinder breaks into two droplets and a smaller satellite droplet in the middle.This type of collision is called ‘‘reflexive separation collision’’.Fig.5shows the calculated resultsforFig.6.Calculated results classified into three types of collisions.The solid and broken curves represent the theoretical prediction of the boundaries between three types of collisions[28].putational domain of two rising bubbles.We ¼79:7and B ¼0:813.Since the two droplets collide at the high impact parameter,only a portion of them contacts directly,and the remaining portions of the droplets tend to move in the direction of their initial velocities and consequently stretch the region of the interaction.Finally,the droplet breaks into two primary droplets and small satellite droplets.This type of collision is called ‘‘stretching separation collision’’.The computation time for the case of Fig.5required about 47h on an AMD AthlonXP 1800+PC machine.The number of iterations in Step 2of the algorithm of computation is large in a very early stage of calculation due to initial large acceleration (e.g.,over 50iterations were required for t Ã<1in the calcu-lation of Fig.5),but becomes smaller than 10after that time.As a result,51%of the total computation time was spent for Step 2in this case.We have calculated the binary droplet collisions for various Weber numbers and impact parameters,and classified the results into the above-mentioned three types of collision in the We –B plane as shown in Fig.6.It is seen that the reflexive separation collisions appear in the region of low impact parameters and high Weber numbers over a critical value,and the stretching separation collisions occur at high impact pa-rameters.The coalescence collisions occur between the two regions.In the figure the theoretical predictions of the boundaries of the three types of collisions [28]are also plotted,and the present calculated results are in good agreement with the theoreticalpredictions.Fig.8.Time evolution of bubble shapes (left)and velocity vectors and density contours on y ¼L y =2(right)for M ¼1Â10À5,E ¼10,q L =q G ¼50;and h M ¼0:5D x (t ütV =D ,where V is the averaged velocity of gas phase in the lower right result).3.3.Bubble flowFirst,we calculate the coalescence of two rising bubbles.In the calculation of bubble flows,the gravi-tational force is considered by adding the term À3E i c iz ð1Àðq L =q ÞÞg D x ,where g is the gravitational ac-celeration to the right-hand side of Eq.(3).Two bubbles with the same diameter D are placed ð4=3ÞD apart in a liquid inside a rectangular domain (see Fig.7)and is released at time t ¼0.The density ratio of the liquid to the gas is q L =q G ¼50(q L ¼50,q G ¼1).The dimensionless parameters for this phenomenon arethe Morton number M ¼g l 4L ðq L Àq G Þ=q 2L r 3and the E €otv €o s number E ¼g ðq L Àq G ÞD 2=r .The periodic boundary condition is used on all the sides of the domain.The domain is divided into an 80Â80Â160cubic lattice.The parameters in Eq.(10)for this example are a ¼1,b ¼1,and T ¼2:93Â10À1;it follows that the maximum and minimum values of the order parameter are /max ¼4:031Â10À1and/min ¼2:638Â10À1.The cut-offvalues of the order parameter are /ÃL ¼3:80Â10À1and /ÃG ¼2:75Â10À1.The other parameters are fixed at s f ¼1,s g ¼1,e ¼10À6,D ¼30D x ,l L =l G ¼50,j f ¼0:05ðD x Þ2and j g ¼1Â10À5ðD x Þ2.The Morton number and the E €otv €o s number are M ¼1Â10À5and E ¼10,respec-tively.The mobilities of h M ¼0:5D x and 0:01D x are used.Figs.8and 9show the calculated results for h M ¼0:5D x and 0:01D x ,respectively.It is seen that as the time passes,the lower bubble catches up with the upper bubble.However,if we look at the results carefully,it is seen that the shapes of the lower bubbles for the two cases are different.It is noted that twobubblesFig.9.Time evolution of bubble shapes (left)and velocity vectors and density contours on y ¼L y =2(right)for M ¼1Â10À5,E ¼10,q L =q G ¼50,and h M ¼0:01D x (t ütV =D ,where V is the averaged velocity of gas phase in the lower right result of Fig.8).。