Axisymmetric lattice Boltzmann model for multiphase flows with large density ratio
aa r X i v : . [ phy s i c s . c o m p - ph ] O c t Axisymmetric lattice Boltzmann model for multiphaseflows with large density ratio
Hong Liang a , Yang Li b , Jiangxing Chen a , Jiangrong Xu a, ∗ a Department of Physics, Hangzhou Dianzi University, Hangzhou, 310018, China b Department of Mathematics, Hangzhou Dianzi University, Hangzhou, 310018, China
Abstract
In this paper, a novel lattice Boltzmann (LB) model based on the Allen-Cahn phase-field theory is proposed for simulating axisymmetric multiphaseflows. The most striking feature of the model is that it enables to handle mul-tiphase flows with large density ratio, which are unavailable in all previousaxisymmetric LB models. The present model utilizes two LB evolution equa-tions, one of which is used to solve fluid interface, and another is adopted tosolve hydrodynamic properties. To simulate axisymmetric multiphase flowseffectively, the appropriate source term and equilibrium distribution functionare introduced into the LB equation for interface tracking, and simultane-ously, a simple and efficient forcing distribution function is also delicatelydesigned in the LB equation for hydrodynamic properties. Unlike many ex-isting LB models, the source and forcing terms of the model arising fromthe axisymmetric effect include no additional gradients, and consequently,the present model contains only one non-local phase field variable, which in ∗ Corresponding author.
Email address: [email protected] (Hong Liang); [email protected] (Jiangrong Xu)
Preprint submitted to Elsevier October 23, 2018 his regard is much simpler. In addition, to enhance the model’s numericalstability, an advanced multiple-relaxation-time (MRT) model is also appliedfor the collision operator. We further conducted the Chapman-Enskog anal-ysis to demonstrate the consistencies of our present MRT-LB model withthe axisymmetric Allen-Cahn equation and hydrodynamic equations. A se-ries of numerical examples, including static droplet, oscillation of a viscousdroplet, breakup of a liquid thread, and bubble rising in a continuous phase,are used to test the performance of the proposed model. It is found that thepresent model can generate relatively small spurious velocities and can cap-ture interfacial dynamics with higher accuracy than the previously improvedaxisymmetric LB model. Besides, it is also found that our present numer-ical results show excellent agreement with analytical solutions or availableexperimental data for a wide range of density ratios, which highlights thestrengths of the proposed model.
Keywords:
Lattice Boltzmann method, axisymmetric flows, multiphaseflows, phase field
1. Introduction
Multiphase fluid flows are ubiquitous in nature and are of considerableinterest in both scientific and engineering fields. In spite of improving ex-perimental studies of these multiphase phenomena, numerical modeling be-comes an increasingly important approach with the rapid development ofcomputer technology and gradual enrichment of computing methods. Andalso, it can provide a convenient access to physical quantities, such as vari-ational interface shapes, fluid velocity across interface, pressure distribution2nside and outside bulk phase, which are usually difficult to measure exper-imentally. Nonetheless, it still remains an intractable task to the model-ing of multiphase flows and further develop efficient numerical algorithmsthat can accurately describe physical phenomena behind such flows. Thereasons behind the challenges lie in the complexity of interface dynamicsamong multispecies fluids, density and viscosity jumps across the interface,and surface tension force modeling. Several numerical methods to date havebeen proposed for simulating multiphase flows, which can be roughly dividedinto two categories: sharp-interface methods, and diffusion-interface meth-ods. The former methods commonly include the volume-of-fluid method [1],front-tracking method [2] and level set method [3]. In this type of traditionalmethods, one might solve continuum mechanics equations coupled with asuitable technique to track the phase interface. The interface needs to becaptured manually based on some complex ad-hoc criteria, and fluid prop-erties in these methods such as density and viscosity vary sharply at theinterface. Therefore, they in this regard are not suitable for handling multi-phase flows with large interface topological change [4, 5]. As for the lattermethods, one replaces the sharp interface with a transition region acrosswhich fluid physical properties are allowed to change smoothly. This featuremakes them more potential for simulating complex interfacial dynamic prob-lems. Among diffusion-interface approaches, the phase-field method [6] andlattice Boltzmann (LB) method [7, 8] are two popular ones. In the phase-fieldmethod, the thermodynamic behavior of a multiphase system is described bythe free energy as a function of an order parameter, which is used to capturephase interface. The interfacial governing equation for the order parameter is3ormulated as the convective Cahn-Hilliard [6, 9] or Allen-Cahn (AC) [10, 11]equation, which is usually solved by using a finite difference like scheme inphase field modelling. Therefore the phase field method will inherit somecertain weaknesses rooted in the conventional numerical methods.Alternatively, the LB method [7, 8] has received great attention for mod-elling hydrodynamic phenomena and in particular for multiphase flows. Itis a mesoscopic method based on the kinetic equation for the particle distri-bution function, connecting the bridge between the macroscopic continuummodel and molecular dynamics. Due to its mesoscopic nature, the LB methodhas several distinct advantages over the traditional numerical methods, suchas the simplicity of algorithmic, nature parallelization and easy implemen-tation of complex boundary. Particularly, the intermolecular interactions ina multiphase system can be incorporated directly in the framework of LBmethod, while they are difficult to handle in traditional numerical methods.Historically, from different physical pictures of the interactions, several typesof LB models have been established for simulating multiphase flows, whichinclude the color-gradient model [12], pseudo-potential model [13], free en-ergy model [14], and phase-field based model [15–18]. Some advanced LBmodels based on these original models have also been proposed in successionand interesting readers can refer to good reviews [19, 20] of this field and thereferences therein.In practice, there exists many multiphase fluid problems that display ax-ial symmetry. Examples include head-on collision of binary droplets [21, 22],bubble rising in a continuous phase [23, 24], droplet formation in micro-channel [25, 26], droplet splashing on a solid surface or wetting liquid film [27,48], and so on. The most natural manner to simulate such axisymmetric flowsis to apply a three-dimensional (3D) LB multiphase model with suitablecurved boundary conditions. This treatment, however, does not take anyadvantages of the axisymmetric property of flow. If we recognize the factthat 3D axisymmetric flows can reduce to the two-dimensional (2D) onesin meridian plane, an effective approach in this regard is to develop quasi-2D LB models for simulating these flows. Up to now, some scholars havemade an effort to construct effective axisymmetric models based on the LBmethod, whereas most of them are proposed for single-phase flows within theisothermal and thermal systems [29–34]. The first axisymmetric LB modelfor multiphase flows was attributed to Premnath and Abraham [35], who in-troduced some suitable source terms in the Cartesian coordinate multiphasemodel of He et al. [15]. The source terms are used to account for the axisym-metric contribution of inertial, viscous and surface tension forces, while theycontains many complex gradients of density and velocity and thus underminethe simplicity of the model. In addition, the model has a drawback inheritedin He’s model that the available highest density ratio is limited to 15. Andalso based on the multiphase model of He et al. [15], Huang et al. [24] pre-sented an improved axisymmetric LB model, in which a mass correction stepis imposed in every numerical iteration. They applied this model to the sim-ulation of bubble rising, and found rather good agreement with experiments.However, the highest density ratio in the simulation is limited to 15.5, andthe model would undergo numerical instability with the increasing densityratio, as they stated. To remove this limitation, Mukherjee and Abraham [36]proposed an axisymmetric counterpart based on the high-density-ratio model5f Lee and Lin [16], in which some source terms relating to the density andvelocity gradients are also included. To improve numerical stability, they alsoused a three-stage mixing discretization scheme, which unfortunately couldinduce the violation of the mass conservation [37]. Later, Huang et al. [38]put forward a hybrid LB model for axisymmetric binary flows, where a fi-nite difference scheme is used to solve the convective Cahn-Hilliard equationfor interface capturing and a multiple-relaxation-time (MRT) LB scheme isadopted for solving the Navier-Stokes (NS) equations. The densities of binaryfluids are supposed to be uniform in the given NS equations, and thereforetheir model in theory can only be applicable for density-matched two-phaseflows. The extension of the popular pseudo-potential model [13] to the ax-isymmetric version was conducted by Srivastava et al. [39]. Recently, Liang et al. [40] proposed a Cahn-Hilliard phase-field based LB model for axisym-metric multiphase flows. One distinct feature of this model is that the addedsource terms representing the axisymmetric effect contain no gradients, thussimplifying the computation. The model is also demonstrated to be accu-rate for simulating multiphase flows with moderate density ratios, while it isunable to tackle large-density-ratio cases. More recently, Liu et al. [41] de-veloped a color-gradient LB model for axisymmetric multicomponent flows.This model can deal with binary fluids with high viscosity ratio, whereas thedensity ratio considered is very small.As reviewed above, several types of LB models have been proposed foraxisymmetric multiphase flows, and all these models are restricted to theflows with small or moderate density contrasts. Generally, the density ratiocan approach 1000 for a realistic liquid-vapor system, and to develop a high-6ensity-ratio multiphase model is very attractive in LB community. In thispaper, we intend to present a simple, accurate and also robust LB model foraxisymmetric multiphase flows, which can tolerate large density contrasts.The proposed LB model is based on the AC phase field theory, which involvesa lower-order diffusion term compared with the Cahn-Hilliard equation, andthus is expected to achieve a better numerical accuracy. Modified equilibriumdistribution functions and simplified forcing distribution functions are alsoincorporated in the model to recover the correct axisymmetric AC and NSequations. Besides, unlike most of available LB models [24, 35, 36, 38, 39],the introduced source terms arising from the axisymmetric effect contain nogradients in the present model. The rest of this paper is organized as follows.In Sec 2, the governing equations and their axisymmetric formulations arefirst given, and then a novel axisymmetric LB model based on the AC phase-field theory are proposed. Numerical validations for the present model canbe found in Sec 3, and at last, we made a brief summary in Sec. 4.
2. GOVERNING EQUATIONS AND MATHEMATICAL MODEL
The AC equation contains only at most a second-order gradient diffusionterm, which can be more efficient and less dispersive in solving the phaseinterface compared with the commonly used Cahn-Hilliard theory [42–45].Therefore the interface tracking equation in this study is built upon the ACphase field theory. The original AC equation was not globally conservative,and recently reformulated into the conservative form [11] based on the workof Sun and Beckermann [10]. This specific formulation is commonly referred7s the conservative phase-field model and will be adopted here. Then theconservative AC equation for interface tracking can be written as [11, 42, 46], ∂φ∂t + ∇ · ( φ u ) = ∇ · [ M ( ∇ φ − λ n )] , (1)where φ is the order parameter acting as a phase field to distinguish differentfluids, u is the fluid velocity, M is the mobility, n is the unit vector normalto the interface and is calculated as n = ∇ φ/ |∇ φ | , λ can be expressed by λ = − φ − φW , (2)where W is the interfacial thickness, the phase field φ takes 1 and 0 in thebulk regions of the liquid and vapor fluids, respectively, and φ = 0 . ∇ · u = 0 , (3a) ∂ ( ρ u ) ∂t + ∇ · ( ρ uu ) = −∇ p + ∇ · (cid:2) νρ ( ∇ u + ∇ u T ) (cid:3) + F s + G , (3b)where p is the hydrodynamic pressure, ν is the kinematic viscosity, G is thepossible external force, F s is the surface tension force. According to Ref. [47],there exists several treatments in terms of surface tension force, which couldgive rise to different performances. To reduce the spurious velocity, in thiswork we choose the widely used the potential form F s = µ ∇ φ [17, 40, 43],where µ is the chemical potential given by µ = 4 βφ ( φ − φ − . − k ∇ φ, (4)where the parameters β and k can be determined by the surface tension σ and the interfacial thickness, i.e., k = 1 . σW , β = 12 σ/W .8e now performed the coordinate transformation to derive the governingequations of isothermal multiphase flow in the axisymmetric system. Thetransformation is given by ( x, y, z ) → ( r, θ, z )with the relations x = r cos θ , y = r sin θ , where r , z , θ denote the radial,axial and azimuthal directions, respectively. The flows are assumed to haveno swirl here and thus we can set the azimuthal velocity and azimuthalcoordinate derivatives to be zero. In this case, the resulting macroscopicequations in the axisymmetric framework can be expressed by ∂ t φ + ∂ α ( φu α ) + φu r r = ∂ α ( M ∂ α φ ) + M ∂ r φr − M ∂ α ( λn α ) − M λ∂ r φr |∇ φ | (5)and ∂ α u α + u r r = 0 , (6a) ∂ t ( ρu β )+ ∂ α ( ρu α u β ) = − ∂ β p + ∂ α νρ ( ∂ α u β + ∂ β u α )+ ˜ F sβ + G β + νρ ( ∂ r u β + ∂ β u r ) r − ρνu r δ βr r − ρu r u α r , (6b)where α, β = { r, z } , δ is the Kronecker function, ˜ F sβ is the modified surfacetension force given by ˜ F sβ = ˜ µ ∇ φ , and ˜ µ = µ − k∂ r φ/r . From Eqs. (5) and(6), it can be clearly found that some additional source terms are generateddue to the axisymmetric effect. These source terms contain some gradients onthe fluid velocity, and seem to be implemented complicatedly. To our knowl-edge, most of previously proposed axisymmetric LB models [24, 29, 32, 34–36, 38, 39] for hydrodynamic flows were constructed based on the governingequation (6), which thus makes them more complex than the standard modelsin the Cartesian coordinate system. Alternatively, we convert the governing9quations (5) and (6) in another form. With some algebraic manipulations,the derived macroscopic equations can be presented as ∂ t ( rφ ) + ∂ α ( rφu α + M φδ αr ) = ∂ α [ M ∂ α ( rφ ) − M rλn α ] (7)and ∂ α ( ru α ) = 0 , (8a) ∂ t ( rρu β )+ ∂ α ( rρu α u β ) = − ∂ β ( rp )+ ∂ α [ rνρ ( ∂ α u β + ∂ β u α )]+ r ( ˜ F sβ + G β )+( p − ρνr u r ) δ βr . (8b)The present LB model will be built upon the macroscopic equations (7) and(8), which utilizes two LB evolution equations, one of them is used to solvethe axisymmetric AC equation, and another for solving the axisymmetricNS equations. It will be demonstrated below that the introduced sourceterms representing the axisymmetric contributions contain no gradients inthe present model. Based on the collision operators, the LB approach can be roughly di-vided into four categories, including the MRT model [48], the two-relaxation-time model [49], the single-relaxation-time model [50], and the entropic LBmodel [51]. To date, these models all have their own impressive versatilityin simulating hydrodynamic flows, while the MRT model in comparison hasits superiority in terms of stability and accuracy, and thus will be used inthe present LB modelling for multiphase flows. The LB evolution equationof the MRT model for the axisymmetric AC equation can be proposed as f i ( x + c i δ t , t + δ t ) − f i ( x , t ) = Ω i ( x , t ) + δ t F i ( x , t ) , (9)10here the collision operator Ω i is defined byΩ i ( x , t ) = − ( M − S f M ) ij [ f j ( x , t ) − f eqj ( x , t )] , (10)where f i ( x , t ) is the phase-field distribution function associated with the dis-crete velocity c i at position x and time t , f eqi ( x , t ) is the equilibrium distri-bution function. To match the target equation, we design a new equilibriumdistribution function as f eqi = ω i (cid:20) rφ + c iα ( rφu α + M φδ αr ) c s (cid:21) , (11)where c s is lattice sound speed, ω i are the weighting coefficients. ω i and c i depend on the choice of the discrete-velocity lattice model. For plane flows,the popular D2Q9 lattice model [17, 18, 50, 52] is adopted here, and ω i canbe then given by ω = 4 / ω − = 1 / ω − = 1 /
36, and c i is defined as c i = (0 , c, i = 0 , (cos[( i − π/ , sin[( i − π/ c, i = 1 − , √ i − π/ π/ , sin[( i − π/ π/ c, i = 5 − , (12)where c = δ x /δ t denotes the lattice speed with δ x and δ t representing the lat-tice spacing and the time step, respectively, and c s = c/ √
3. By convention, δ x and δ t in the study of multiphase flows are set to be length and time units,i.e., δ x = δ t = 1. In Eq. (10), M is the transformation matrix, which is usedto project the distribution function in the discrete-velocity space onto theones in the moment space. Based on the D2Q9 lattice structure, it can be11iven by [48] M = − − − − − − − − − − − − − − − − − −
10 0 − − −
10 1 − − − − , and S f in Eq. (10) is a diagonal relaxation matrix, S f = diag ( s f , s f , s f , s f , s f , s f , s f , s f , s f ) , (13)where 0 ≤ s fi <
2. To simplify the LB algorithm, one part of the diffusionterm in Eq. (7) is treat as the source term here, and then a novel discretesource term is introduced as F i = [ M − ( I − S f M ] ij R j , (14)where I is the unit matrix, and R i is defined by R i = ω i c iα [ ∂ t ( rφu α + M φδ αr ) + c s rn α λ ] c s . (15)In the present model, the phase field variable φ is derived by the summationof the distribution function f i , and then is computed by φ = 1 r X i f i . (16)12hysically, the density distribution in a multiphase system is consistent withthat of the phase field variable. To satisfy this property, we take the linearinterpolation scheme to determine the fluid density, ρ = φρ l + (1 − φ ) ρ g , (17)where ρ l and ρ g denote the liquid and gas fluid densities, respectively. We alsocarried out the Chapman-Enskog analysis to demonstrate the consistency ofthe present model with the target equation (7). Applying the multiscaleexpansions to Eq. (9), it is found in Appendix A that the axisymmetric ACequation can be recovered exactly from the present model, and the mobilityis given by M = c s δ t ( τ f −
12 ) , (18)where 1 /τ f = s f = s f . In comparison with the standard LB model [42] forAC equation, it is found that the introduced source terms to account for theaxisymmetric effect include no additional gradient in the present model.We now give a discussion on implementation of the MRT model. Gen-erally, the MRT-LB equation (9) can be solved in two steps, including thecollision process, f + i = f i ( x , t ) − ( M − S f M ) ij [ f j ( x , t ) − f eqj ( x , t )] + δ t F i ( x , t ) (19)and the propagation process, f i ( x + c i δ t , t + δ t ) = f + i . (20)To reduce the matrix operations, it is wise that the collision process of MRTmodel is conducted in the moment space. By premultiplying the matrix M ,13he equilibrium distribution function in moment space is derived as Mf eq = (cid:18) rφ, − rφ, rφ, rφu z c , − rφu z c , rφu r + M φc , − rφu r + M φc , , (cid:19) T , (21)where u r and u z are the radial and axial components of velocity. Similarly,the discrete source term R i in the moment space can also be obtained as MR = (0 , , , mR , − mR , mR , − mR , , T , (22)where mR and mR are given by mR = ∂ t ( rφu z ) + c s rn z λc , mR = ∂ t ( rφu r + M φ ) + c s rn r λc . (23)At the end of this subsection, we also would like to stress that the presentMRT model for the axisymmetric AC equation can reduce to the SRT coun-terpart when the relaxation factor s fi in Eq. (13) equals to each other, andthus the SRT model is only one special case of the MRT model. In thiscase, the more freedom in the choice of the relaxation factors can providemore potential for the MRT model to achieve better numerical accuracy andstability. The LB evolution equation with a generalized collision matrix for theaxisymmetric NS equations can be written as g i ( x + c i δ t , t + δ t ) − g i ( x , t ) = − (cid:2) M − S g M (cid:3) ij (cid:2) g j ( x , t ) − g eqj ( x , t ) (cid:3) + δ t G i ( x , t ) , (24)where g i is the hydrodynamic distribution function, g eqi is the correspondingequilibrium distribution function. To incorporate the axisymmetric effect, a14odified form of the equilibrium distribution function is used [40], g eqi = rpc s ( ω i −
1) + rρs i ( u ) , i = 0 rpc s ω i + rρs i ( u ) , i = 0 (25)with s i ( u ) = ω i " c i · u c s + ( c i · u ) c s − u · u c s . (26)For fluids exposed to forces, the discrete lattice effects should be taken intoaccount, when the forces are introduced in LB approach [53], and then thediscrete force term in the MRT framework for hydrodynamics is defined by G i = (cid:20) M − ( I − S g M (cid:21) ij T j , (27)where T i is the forcing distribution function, and S g is a non-negative diagonalmatrix given by S g = diag ( s g , s g , s g , s g , s g , s g , s g , s g , s g ) . (28)Different from other axisymmetric LB multiphase models [24, 35, 36, 38–41],a much simplified forcing distribution function is constructed in this model, T i = ω i (cid:20) c iα F α c s + u α ∂ β ( rρ ) c iα c iβ c s − ρu r (cid:21) , (29)where F α is the total force and is given by F α = r ( e F sα + G α ) + ( p − ρνr u r ) δ αr . (30)Taking the first-order moment of the distribution function, the fluid velocityin this model can be evaluated as rρu α = X i c iα g i + 0 . δ t F α , (31)15hich can be further recast explicitly, u α = P i c iα g i + 0 . δ t [ r ( ˜ F sα + G α ) + pδ αr ] rρ + δ t r − νρδ αr . (32)The hydrodynamic pressure p can be calculated in a particular manner. Asshown in Appendix C , it can be evaluated as p = c s − ω " r X i =0 g i + δ t u · ∇ ρ + ρs ( u ) − δ t s p ω ρu r r , (33)where s p is a parameter relating to the relaxation times, s p = s g + s g − s g s g s g s g .We also conducted the Chapman-Enskog analysis on the present model foraxisymmetric NS equations, and it is demonstrated in Appendix B that thecorrect hydrodynamic equations in the cylindrical coordinate system can alsobe derived from the present model, and the relation between the kinematicviscosity and the relaxation factor can be presented as ν = c s δ t ( τ g −
12 ) , (34)where 1 /τ g = s g = s g = s g . Generally, the relation s g = s g is satisfied inprevious MRT model for hydrodynamics. Here this new constraint for therelaxation factors is derived to account for the axisymmetric effect. In a mul-tiphase system, the fluid viscosity is no longer uniform due to its sharp jumpat the liquid-gas interface. For simplicity, the popular linear interpolationscheme is used to determine the viscosity at the interface [15], ν = φν l + (1 − φ ) ν g , (35)where ν l and ν g represent the kinematic viscosities of the liquid and gasphases, respectively. As for the MRT hydrodynamic model, the collision16rocess can also be implemented in moment space. With some algebraic ma-nipulations, the equilibrium and forcing distribution functions in the momentspace can be respectively given by Mg eq = r (cid:18) , p + ρ u c s , − p + ρ u c s , ρu z c , − ρu z c , ρu r c , − ρu r c , ρ ( u z − u r ) c , ρu r u z c (cid:19) T (36)and MT = (cid:18) ru α ∂ α ρ, ρu r , − u α ∂ α rρ − ρu r , F z c , − F z c , F r c , − F r c ,
23 ( u z ∂ z rρ − u r ∂ r rρ ) ,
13 ( u z ∂ r rρ + u r ∂ z rρ ) (cid:19) T . (37)In practice, the time derivative and the spatial gradients should be dis-cretized with suitable difference schemes. In this paper, we adopt the explicitEuler scheme to calculate the time derivative as follows, ∂ t χ ( x , t ) = χ ( x , t ) − χ ( x , t − δ t ) δ t . (38)Following the work of Liang et al. [40], the gradient and the Laplacian opera-tor are discretized with the following four-order isotropic difference schemes, ∇ χ ( x , t ) = X i =0 ω i c i [8 χ ( x + c i δ t , t ) − χ ( x + 2 c i δ t , t )]6 c s δ t , (39)and ∇ χ ( x , t ) = X i =0 ω i [16 χ ( x + c i δ t , t ) − χ ( x + 2 c i δ t , t ) − χ ( x , t )]6 c s δ t , (40)where χ represents an arbitrary variable. The evaluation of the gradientterms at fluid nodes neighbouring to boundary, using Eqs. (39) and (40),typically requires unknown information at the ghost nodes, and we use themirror symmetric rule to derive the unknown values for fluid nodes neigh-bouring to solid wall, while applying the axis-symmetric means for the axial17odes. In addition, the boundary conditions for the distribution functionsshould also be specially treat, since the singularity arises at the symme-try axis of r = 0. To avoid this problem, we set the first lattice line at r = 0 . δ x , and apply the symmetry boundary condition for axial bound-ary. For other solid boundary, we impose the no-slip bounce back boundarycondition. The detailed treatments on the used boundary conditions anddiscretization schemes can be also referred to Refs. [31, 40].
3. NUMERICAL VALIDATION FOR AXISYMMETRIC LB MUL-TIPHASE MODEL
In this section, we will test the accuracy and stability of the proposedaxisymmetric LB model by using several basic tests. These typical numericalexamples include the static droplet, oscillation of a droplet, breakup of aliquid thread, and bubble rising in a continuous phase, where a very boardrange of density ratios are considered to highlight the advantage of the presentLB model.
A benchmark problem of the static droplet is first used to verify thepresent LB model for axisymmetric multiphase flows. Initially, a stationarydroplet is located at the computational domain with the size of N z × N r =200 × ,
0) and its radius ( R ) occupies50 lattice units. The left and right boundary conditions for both f i and g i are set to be periodic, while the symmetry boundary condition is applied atthe axial line and the no-slip boundary condition is imposed at the upper18oundary. To be smooth across the interface, the phase field variable isinitialized by φ ( z, r ) = 0 . . " R − p ( z − + r W , (41)and the corresponding distribution of the density field based on Eq. (17) canbe given by ρ ( z, r ) = ρ l + ρ g ρ l − ρ g " R − p ( z − + r W , (42)where ρ l and ρ g in the test are set to be 1000 and 1, corresponding to alarge density ratio of 1000. Some other simulation parameters are given asfollows: W = 4 . σ = 0 . M = 0 .
01 and ν l = ν g = 0 .
1. Considering thecomponents of the equilibrium in moment space, we fix s f = s f = s f = 1 asusual, and the relations s f = s f and s f = s f = s f = s f are satisfied. Forsimplicity, the values of s f and s f in the simulation are set to be 1. Therelaxation factors s f and s f are given by the value of 1 /τ f , which is furtherdetermined by the mobility M based on Eq. (18). As for the ones in thematrix S g , the relaxation times s gi are chosen to be unity expect for s g , s g and s g , which are adjusted according to the fluid viscosity.Figure 1(a) depicts the steady droplet shape obtained by the present LBmodel, together with the initial one. It is found that they match to each otherwith a high accuracy. Further, we also plotted the numerical prediction of thedensity profile across the interface in Fig. 1(b), where the analytical solutiongiven by Eq. (42) is as well presented for the comparison. It can be clearlyobserved that numerical result agrees well with the analytical solution, whichindicates that the present model can correctly solve the phase interface. The19
50 100 150 200−100−50050100 z r (a) r ρ Analytical solutionLB simulation ρ g ρ l (b) Figure 1: static droplet test with density ratio ρ l /ρ g = 1000: (a) the velocity distributionof the whole domain at the equilibrium state and the solid and dashed lines respectivelyrepresent the equilibrium shape of the droplet and its initial shape; (b) density profileacross the interface obtained from LB simulation and corresponding analytical solution. existence of the spurious velocities is a common undesirable feature in manynumerical methods for multiphase flows. If they have comparable magnitudesas the real fluid velocity, numerical instability or some unphysical phenomenacould take place, and thus achieving small spurious velocities is very signi-ficative. The spurious velocities also rooted in the LB method arise fromthe imbalance between discrete forces in the interfacial region and cannot betotally eliminated to round-off in the framework of the LB approach [54]. Wealso displayed the spurious velocities generated by the present LB model. Thedistribution of the velocity field in the whole computational domain is shownin Fig. 1(a). It can be obviously found that the spurious velocities indeedexist in the vicinity of the interface, and the maximum magnitude of the spu-rious velocities computed by | u | max = ( √ u + v ) max is about 1 × − . We20lso compared the spurious velocities generated by different LB approaches.The axisymmetric color-gradient LB multiphase model proposed by Liu etal. [41] can only deal with multiphase flows with moderate density ratio andthey reported that it can derive the spurious velocities at the level of 10 − . Inaddition, it is found that the maximum amplitude of spurious velocities in animproved pseudo-potential model [55] has the order of 10 − . In comparisonswith these common LB models, we can conclude that the present phase-fieldbased LB model can derive relatively low spurious velocities. The droplet oscillation is a classic example that is widely used to validateaxisymmetric multiphase LB model for simulating dynamic problem [35, 39–41]. The liquid droplet will exhibit oscillatory behavior, if it is initiallydistorted from equilibrium spherical shape into an elliptical one. We intendto compare the droplet oscillation frequency obtained from the LB simulationwith the analytical solution reported by Miller and Scriven [56]. Accordingto their analysis, the theoretical prediction for the oscillation frequency ofthe n th mode can be given by ω n = ω ∗ n − . α p ω ∗ n + 0 . α , (43)where ω ∗ n is Lamb’s natural resonance frequency, ω ∗ n = s n ( n − n + 1)( n + 2) R e [ nρ g + ( n + 1) ρ l ] σ, (44)where R e is the radius of droplet at the equilibrium state. In Eq. (43),the parameter α is included to account for the viscosity contribution and is21xpressed by [56] α = (2 n + 1) ρ l ρ g √ ν l ν g √ R e [ nρ g + ( n + 1) ρ l ]( ρ l √ ν l + ρ g √ ν g ) . (45)The mode of the oscillation is denoted by n , which is set to be 2 for an initialellipsoidal shape considered here.The initial computational setup consists of an ellipsoidal droplet with thehalf -axial and -radial lengths denoted by R z and R r , placed in a domainwith the size of N z × N r = 300 × φ ( z, r ) = 0 . . " R e − p ( z − z c ) /R z + ( r − r c ) /R r W , (46)where ( z c , r c ) = (150 ,
0) is the center coordinate of ellipsoidal droplet, andthe equilibrium droplet radius R e can be calculated as ( R r R z ) / based onthe mass conservation of the liquid droplet. The density ratios used in thistest range from 10 and 100, and we find that the model is still stable for avery large density ratio of 1000. Considering a very low oscillation frequencyat a large density ratio, it is generally measured difficultly, and thus theresult in this case is not presented. Some remaining physical parameters inthe simulation are fixed as W = 4, σ = 0 . M = 0 . ν l = ν g = 0 .
1, andthe relaxation factors in the matrices S f and S g are set to be those in thelast test.Figure 2 shows the snapshots of an oscillating droplet at different times,where the density ratio is set as 100 and the initial droplet size is given by R z = 90 and R r = 30. As can be seen from Fig. 2, the droplet changes22 =0.0 t=8600 t=15200 t=19300 Figure 2: Evolution of the shape of an oscillating droplet with ρ l /ρ g = 100, R r = 30, R z = 90, σ = 0 . t R r / R e R e =43.3R e =56.7R e =68.7R e =79.7 Figure 3: Time evolution of the half-axis length R r at different equilibrium radii with ρ l /ρ g = 100, σ = 0 .
3. And the half-axis length R r is normalized by the equilibrium radius R e . t R r / R e ρ l / ρ g =100 ρ l / ρ g =50 ρ l / ρ g =20 ρ l / ρ g =10 Figure 4: Time evolution of the half-axis length R r at different density ratios with R r = 30, R z = 90, σ = 0 .
3. And the half-axis length R r is normalized by the equilibrium radius R e . able 1: Comparison between the computed oscillating frequency ω LB and correspondinganalytical values ω ana at different equilibrium radii for ρ l /ρ g = 100, σ = 0 . R z = 90. R e ω LB ( × − ) 4.8934 3.4074 2.5964 2.1313 ω ana ( × − ) 5.3764 3.5864 2.6908 2.1533 E r = | ω LB − ω ana | ω ana × t = 8600. Suchchange continues until the droplet reaches its equilibrium spherical shape.The above droplet behaviors are consistent with the expectation. We alsoconducted a quantitative study of droplet oscillation problem. The variationof the half-axis length R r versus time is plotted in Fig. 3, where the resultswith other initially given values of R r = 45, 60 and 75 are also presentedto examine the effect of droplet size. It can be found from Fig. 3 thatthe amplitude of the oscillation normalized by the corresponding equilibriumradius R e fluctuates around 1 for all cases, and as expected, the dimensionlessmaximum amplitude is reduced for a larger droplet size. Further, from Fig.3 we can compute the droplet oscillating frequency and showed the resultsin Table 1. For a comparison, the theoretical solution for predicting theoscillating frequency was also presented. It is found that they agree well ingeneral, with the maximum relative error of around 8.9%.The influence of the density ratio on the oscillating frequency is alsoinvestigated here. We use the present model to simulate this case with fourdifferent density ratios ρ l /ρ g = 10, 20, 50 and 100, and the droplet size is fixedto be R z = 90 and R r = 30. Figure 4 depicts the evolution of the oscillating25 able 2: Comparison between the computed oscillating frequency ω LB and correspondinganalytical values ω ana at different density ratios for R z = 90, R r = 30, and σ = 0 . ρ l /ρ g
10 20 50 100 ω LB ( × − ) 15.2504 11.0000 7.0439 4.8934 ω ana ( × − ) 16.5854 12.0000 7.5895 5.3764 E r = | ω LB − ω ana | ω ana × To show the capacity of the current model in simulating large interfacetopological changes, in the subsection we consider a fascinating problem ofthe breakup of a liquid thread into satellite droplets. The breakup of liquidfilaments is of long-standing importance not only in its own interests, but alsoin the fact that numerous practical applications, such as gene chip arraying,ink-jet printing and microfluidics, depend critically on the knowledge of thebreakup mechanisms [57, 58]. The earliest experimental study on the breakupof liquid filament was attributed to Plateau [59], and later Rayleigh [60]26erformed the linear stability analysis on this problem, who argued that acylindrical liquid thread of radius R is unstable, if the wavelength ˜ λ of adisturbance on thread surface is longer than its circumference 2 πR . Usually,this instability criterion can also be described by using the wave number k ,i.e., k = 2 πR/ ˜ λ . When k is smaller than 1, a liquid thread is unstable andcan separate into small droplets, while it is stable for the case of k >
1. Somepreliminary numerical experiments with k > k <
1, and compare the present numerical resultswith some available data.The simulations were carried out in a N z × N r = ˜ λ ×
200 lattice domainwith the same boundary conditions used in the static and oscillating droplettests. An initial disturbance imposed on thread surface is given by settingthe phase field profile as φ ( z, r ) = 0 . . R + d − r ) W , (47)where the radius R is taken as 60, and d is the perturbation function given by d = 0 . R cos(2 πz/ ˜ λ ). To examine the wave number effect, several cases withdifferent wavelengths ˜ λ =420, 500, 600, 800, 1000, 1200, 1800 were simulated,which correspond to various wave numbers k =0.90, 0.75, 0.63, 0.47, 0.38,0.31, 0.21. In the simulations, two different density ratios of ρ l /ρ g = 100 and10 are considered and the interfacial tension σ is fixed as 0.3. We would liketo stress that the model can also tackle this case with a high density ratio of1000, while it takes plenty of time for the system to reach the breakup stateunder the identical surface tension force. In this case, the moderate density27 .09.412.312.513.013.415.1 (a) (b) (c) Figure 5: Evolutions of the breakup of a liquid thread at different wave numbers with thedensity ratio of ρ l /ρ g = 100, (a) k = 0 .
63, (b) k = 0 .
31, (c) k = 0 .
21. Time has beennormalized by the capillary time p R ρ l /σ . ratio is merely considered. Here some other used physical parameters aregiven as W = 4 . M = 0 .
01 and ν l = ν g = 0 .
1. The parameters in both S f and S g are chosen as those in the last test. We first present the resultsfor the density ratio of 100. Figure 5 depicts the time evolution of a liquidthread with three typical wave numbers, where time has been normalizedby the capillary time p R ρ l /σ . It can be observed from Fig. 5 that forall cases, the interfacial perturbation grows continuously at early stage, andthen the liquid thread in middle region becomes more and more thin, whileits ends are enlarged with time. As time goes on, the liquid thread breaks upat two thin linking points, forming into a liquid ligament as well as a main28roplet. Afterwards, the liquid ligament shrinks constantly with the action ofthe dominant surface tension, until it reaches an equilibrium spherical state.Lastly, a pair of liquid droplets including the main droplet and the satelliteone can be observed in the system. Actually, the formed liquid ligament canbe treat as the secondary liquid thread, and whether its breakup occurs ornot depends on the new wavelength. If the wavelength is sufficiently largesuch that the instability criterion can be satisfied, the secondary breakupof the ligament can take place, leading to the formation of multiple satellitedroplets. The above behaviors of the thread breakup obtained by the presentmodel are qualitatively consistent with the previous studies [61, 62].To further investigate the density ratio effect, we also simulate the breakupof a liquid thread with the density ratio of ρ l /ρ g = 10, where other physi-cal parameters remain unchanged. The snapshots of the breakup of a liquidthread with three different wave numbers are shown in Fig. 6. It is foundthat compared with the case of ρ l /ρ g = 100, the liquid thread exhibits simi-lar behaviors at early time for all the wave numbers. The initial disturbanceimposed on the thread surface increases in time, which then results in theformation of a thin ligament and a main droplet. Afterwards, some distinctbehaviors of a liquid ligament are observed for different wave numbers. Forthe case of k = 0 .
62 or 0.31, the filament contracts continuously into onesatellite droplet. This phenomenon is a familiar spectacle in simulation re-sults with the density ratio of 100. Whereas for the case of k = 0 .
21, the liquidligament shrinks at first, and then pinches off due to the Rayleigh instability,which leads to the relaxation of two daughter droplets. The two daughterdroplets move toward each other, and eventually merge into a larger satellite29 .08.910.211.011.411.624.2 (a) (b) (c)
Figure 6: Evolutions of the breakup of a liquid thread at different wave numbers withthe density ratio of ρ l /ρ g = 10, (a) k = 0 .
63, (b) k = 0 .
31, (c) k = 0 .
21. Time has beennormalized by the capillary time p R ρ l /σ . k = 0 .
63, 0.31,0.21 are 12.3, 18.9, 27.1, while they for the density ratio of 10 are 10.2, 16.8,23.9. Therefore it is concluded that increasing density ratio can slow downthe breakup of a liquid thread. Furthermore, we also give a quantitativestudy on the thread breakup, and showed in Fig. 7 the main and satellitedroplet sizes with two density ratios and various wave numbers. For compar-isons, some available literature results, including the finite element results ofAshgriz and Mashayek [61], the analytical solutions and experimental data ofLafrance [63], as well as the experimental data by Rutland and Jameson [64],are also presented. It can be seen that both the main and satellite dropletsizes decrease with the increasing wave number, and also, it is found thatthe satellite droplet size is reduced for a larger density ratio, whereas thetrend is just the reverse for the main droplet due to the mass conservation.The numerical predictions for the droplet sizes obtained from the present LBsimulations are further found to be in good agreement with these availabledate in general. 31 k R * Ashgriz [61], Num., MainAshgriz [61], Num., SatelliteLafrance [63], Anal., MainLafrance [63], Anal., SatelliteLafrance [63], Exp., MainLafrance [63], Exp., SatelliteRutland [64], Exp., MainRutland [64], Exp., SatelliteLB simulation, MainLB simulation, SatelliteLB simulation, MainLB simulation, Satellite
Figure 7: The terminal droplet sizes at different wave numbers. Here the red circle andblue rhombus respectively represent the cases of ρ l /ρ g = 10 and ρ l /ρ g = 100, and thedroplet radius R ∗ has been normalized by R .
4. Summary
Numerical modeling of axisymmetric multiphase flows with large densitycontrasts is still an intractable task in the framework of the LB method. Inthis paper, we propose a simple and efficient LB model for axisymmetricmultiphase system, which is capable of handling large density differences.The proposed LB model is built upon the conservative phase-field equa-tion, which involves a lower-order diffusion term as opposed to the widelyused Cahn-Hilliard equation. Therefore the present model in theory canachieve better numerical stability and accuracy in solving phase interfacethan the Cahn-Hilliard type of axisymmetric LB models. Two LB evolutionequations are utilized in the current model, one of which is used for cap-turing phase interface and the other for solving fluid velocity and pressure.32 new equilibrium distribution function and some discrete source terms aredesigned in the LB equation for interface capturing. Meanwhile, a simplerforcing distribution function is also proposed in the LB equation for hy-drodynamics. Different from most of previous axisymmetric LB multiphasemodels [24, 35, 36, 38, 39], the added source terms accounting for the ax-isymmetric effect contain no gradients in the present model. The presentmodel is also equipped with an advanced MRT collision operator to enhanceits numerical stability. We conducted the Chapman-Enskog analysis on thepresent MRT-LB equations, and it is demonstrated that both the conser-vative AC equation and the incompressible NS equations in the cylindricalcoordinate system can be recovered correctly from the present model. Tovalidate the present model, we first simulated a steady problem of the staticdroplet with a high density ratio of 1000, and it is found that the presentmodel can accurately solve the density field, and also achieve low spuriousvelocities with the order of 10 − . To further assess the current model, twodynamic benchmark problems of a droplet oscillation and a liquid threadbreakup are considered. It is shown that the present model can provide sat-isfactory predictions of the droplet oscillation frequency and the generateddaughter droplet sizes for a board range of density ratios. At last, we simu-lated the realistic buoyancy-driven bubbly flow with a large density ratio of1000. Some fascinating bubble dynamics are successfully reproduced by thepresent model, and the numerically predicted terminal bubble patterns showgood agreements with the experimental date and some available numericalresults. It is also reported that the present model can describe bubble inter-facial dynamics with a higher accuracy than the existing axisymmetric LB33odel [24]. The present method is developed based on the standard orthog-onal MRT model, and its extension to the non-orthogonal MRT one can beconducted directly. It is expected that the non-orthogonal model can retainthe numerical accuracy while simplifying the implementation [69, 70]. Inconclusion, we anticipate that our present numerical method will be usefulfor many practical and sophisticated problems. Acknowledgments
This work is also financially supported by the National Natural ScienceFoundation of China under Grants Nos. (11602075, 51776068 and 11674080),and Natural Science Foundation of Zhejiang Province (No. LR17A050001).
Appendix A. Chapman-Enskog analysis of axisymmetric LB modelfor the Allen-Cahn equation
The Chapman-Enskog analysis is now carried out to demonstrate the con-sistency of the LB evolution equation (9) with the axisymmetric AC equation.We first introduce the multiscaling expansions for the distribution function,time derivative, space gradient, and discrete source term as f i = f (0) i + ǫf (1) i + ǫ f (2) i + · · · , (A.1a) ∂ t = ǫ∂ t + ǫ ∂ t , ∇ = ǫ ∇ , F i = ǫF i (1) + ǫ F i (2) , (A.1b)where ǫ is a small expansion parameter. Applying Taylor expansion to Eq.(9) about x and t , one can derive the resulting continuous equation D i f i + δ t D i f i + · · · = − δ t ( M − S f M ) ij ( f j − f eqj ) + F i , (A.2)34here D i = ∂ t + c i · ∇ . Substituting the expansions into Eq. (A.2), we canobtain the following equations in consecutive order of the parameter ǫ , ǫ : f (0) i = f ( eq ) i , (A.3a) ǫ : D i f (0) i = − δ t ( M − S f M ) ij f (1) j + (cid:20) M − ( I − S f M (cid:21) ij R j (1) , (A.3b) ǫ : ∂ t f (0) i + D i f (1) i + δ t D i f (0) i = − δ t ( M − S f M ) ij f (2) j + (cid:20) M − ( I − S f M (cid:21) ij R j (2) , (A.3c)in which D i = ∂ t + c i · ∇ . Rewriting Eqs. (A.3) in vector form and premul-tiplying the matrix M on both sides of them, we can derive the correspondingequations in the moment space ǫ : m (0) f = m ( eq ) f , (A.4a) ǫ : ˆ D m (0) f = − S f ′ m (1) f + ( I − S f MR (1) , (A.4b) ǫ : ∂ t m (0) f + ˆ D ( I − S f m (1) f + δ t D ( I − S f MR (1) = − S f ′ m (2) f +( I − S f MR (2) , (A.4c)where S f ′ = S f /δ t , ˆ D = MD M − , D = ∂ t I + ∇ · diag ( c , c , · · · , c ), m (1) f = ( m (1) f , m (1) f , · · · , m (1) f ) and Eq. (A3.b) has been used to derive theresulting equation (A.4c). Substituting Eq. (21) to Eq. (A.4b), we canobtain several equations related to the target governing equations, ∂ t ( rφ ) + ∂ z ( rφu z ) + ∂ r ( rφu r + M φ ) = 0 , (A.5a) ∂ t ( rφu z ) + ∂ z ( c s rφ ) = − cs f ′ m (1) f + c (1 − s f mR (1)1 , (A.5b) ∂ t ( rφu r + M φ ) + ∂ r ( c s rφ ) = − cs f ′ m (1) f + c (1 − s f mR (1)2 . (A.5c)35imilarly, the substitution of Eq. (21) into Eq. (A.4c) yields ∂ t ( rφ )+ ∂ z c (1 − s f m (1) f + ∂ r c (1 − s f m (1) f + δ t ∂ z c (1 − s f mR (1)1 + δ t ∂ r c (1 − s f mR (1)2 = 0 , (A.6)where the unknown m (1) f and m (1) f can be determined from Eqs. (A.5b-A.5c),and then their expressions can be given by cs f ′ m (1) f = − ∂ t ( rφu z ) − ∂ z ( c s rφ ) + c (1 − s f mR (1)1 , (A.7a) cs f ′ m (1) f = − ∂ t ( rφu r + M φ ) − ∂ r ( c s rφ ) + c (1 − s f mR (1)2 . (A.7b)Substituting Eqs. (A7) to Eq. (A6), we have ∂ t ( rφ ) = ∂ z c s δ t ( 1 s f −
12 ) (cid:2) ∂ z ( rφ ) − rλn (1) z (cid:3) + ∂ r c s δ t ( 1 s f −
12 ) (cid:2) ∂ r ( rφ ) − rλn (1) r (cid:3) , (A.8)where Eq. (23) has been applied. Combining Eq. (A.5a) at t time scale andEq. (A.8) at t time scale, we can derive the recovered equation, ∂ t ( rφ ) + ∂ α ( rφu α + M φδ αr ) = ∂ α ( M ∂ α rφ ) − M ∂ α ( rλn α ) , (A.9)where the relation s f = s f has been used to obtain the isotropic mobility, M = c s δ t ( 1 s f −
12 ) . (A.10)From the above procedure, it is clear that the axisymmetric AC equation canbe recovered correctly from the present MRT model without adopting anyapproximations. Appendix B. Chapman-Enskog analysis of axisymmetric LB modelfor the incompressible hydrodynamic equations
The present MRT-LB model for the axisymmetric NS equations is alsoanalyzed by applying the Chapman-Enskog expansions, and similarly the36article distribution function, time and space derivatives, forcing distributionfunction can be expanded as g i = g (0) i + ǫg (1) i + ǫ g (2) i + · · · , (B.1a) ∂ t = ǫ∂ t + ǫ ∂ t , ∇ = ǫ ∇ , G i = ǫG i (1) . (B.1b)Adopting Taylor expansion to Eq. (24) and substituting the expansions (B.1)into the resulting equation, one can get the zero-, first-, and second-orderequations in the parameter ǫ , ǫ : g (0) i = g ( eq ) i , (B.2a) ǫ : D i g (0) i = − δ t ( M − S g M ) ij g (1) j + G i (1) , (B.2b) ǫ : ∂ t g (0) i + D i g (1) i + δ t D i g (0) i = − δ t ( M − S g M ) ij g (2) j . (B.2c)Substituting Eq. (B.2b) into Eq. (B.2c) and multiplying matrix M on theboth sides of Eqs. (B2) separately gives ǫ : m (0) g = m ( eq ) g , (B.3a) ǫ : ˆ D m (0) g = − S g ′ m (1) g + ( I − S g MT (1) , (B.3b) ǫ : ∂ t m (0) g + ˆ D ( I − S g m (1) g + δ t D ( I − S g MT (1) = − S g ′ m (2) g , (B.3c)where S g ′ = S g /δ t . Now it is seen that the distribution functions in thediscrete-velocity space can be projected onto the macroscopic quantities inthe moment space. Based on the moment conditions, we can define m g as m g = Mg = ( − δ t r u ·∇ ρ, m g , m g , rρu z c − δ t F z c , m g , rρu r c − δ t F r c , m g , m g , m g ) T , (B.4)37nd m ( k ) g ( k ≥
1) can be denoted by m (1) g = ( − δ t r u · ∇ ρ, m (1) g , m (1) g , − δ t F z (1) c , m (1) g , − δ t F r (1) c , m (1) g , m (1) g , m (1) g ) T , (B.5a) m ( k ) g = (0 , m ( k ) g , m ( k ) g , , m ( k ) g , , m ( k ) g , m ( k ) g , m ( k ) g ) T , ( k ≥ . (B.5b)Substituting Eqs. (36-37) and Eq. (B.5a) into Eq. (B.3b), we can obtainseveral macroscopic equations at the ǫ scale, and only ones used to recoverthe hydrodynamic equations are presented here, ∂ z ( ru z ) + ∂ r ( ru r ) = 0 , (B.6a) ∂ t ( rρu z ) + ∂ z ( rp + rρu z ) + ∂ r ( rρu z u r ) = F (1) z , (B.6b) ∂ t ( rρu r ) + ∂ z rρu z u r + ∂ r ( rp + rρu r ) = F (1) r , (B.6c) ∂ t rp + rρu z + rρu r c s = − s g ′ m (1) g + T g , (B.6d) ∂ t rρu z − rρu r c + ∂ z rρu z − ∂ r rρu r − s g ′ m (1) g + T g , (B.6e) ∂ t rρu z u r c + ∂ z rρu r ∂ r rρu z − s g ′ m (1) g + T g , (B.6f)where T g = 2(1 − s g ρu r ) , (B.7a) T g = 23 (1 − s g u z ∂ z ( rρ ) − u r ∂ r ( rρ )] , (B.7b) T g = 13 (1 − s g u z ∂ r ( rρ ) + u r ∂ z ( rρ )] . (B.7c)Similarly, we substitute Eqs. (36-37) and Eqs. (B.5) into Eq. (B.3c). Withsome algebraic operations, the moment equations at the ǫ scale can be de-rived and the related ones are presented as ∂ t ( rρu z )+ c ∂ z (1 − s g m (1) g + c ∂ z (1 − s g m (1) g + c ∂ r (1 − s g m (1) g + T gz = 0 , (B.8a)38 t ( rρu r )+ c ∂ r (1 − s g m (1) g − c ∂ r (1 − s g m (1) g + c ∂ z (1 − s g m (1) g + T gr = 0 , (B.8b)where T gz and T gr are given by T gz = ∂ z (1 − s g c δ t ρu r ) + ∂ z (1 − s g c δ t u z ∂ z ( rρ ) − u r ∂ r ( rρ )]+ ∂ r (1 − s g c δ t u z ∂ r ( rρ ) + u r ∂ z ( rρ )] , (B.9)and T gr = ∂ r (1 − s g c δ t ρu r ) − ∂ r (1 − s g c δ t u z ∂ z ( rρ ) − u r ∂ r ( rρ )]+ ∂ z (1 − s g c δ t u z ∂ r ( rρ ) + u r ∂ z ( rρ )] , (B.10)In Eqs. (B.8), the variables m (1) g , m (1) g and m (1) g needs to be determined.Based on the relations (B.6), one can get − s g ′ m (1) g = ∂ t rp + rρu z + rρu r c s − − s g ρu r , (B.11a) − s g ′ m (1) g = ∂ t rρu z − rρu r c + s g u z ∂ z ( rρ ) − u r ∂ r ( rρ )]+ 23 rρ [ ∂ z u z − ∂ r u r ] , (B.11b) − s g ′ m (1) g = ∂ t rρu z u r c + s g u z ∂ r ( rρ ) + u r ∂ z ( rρ )] + 13 rρ [ ∂ r u z + ∂ z u r ] . (B.11c)Substituting Eqs. (B.11) into Eqs. (B.8) and using the incompressible limits[ ∂ t p = O ( M a ) and u = O ( M a )], we can ultimately derive the second-order equations in ǫ , ∂ t ( rρu z ) − ∂ z rνρ ( ∂ z u z − ∂ r u r − u r r ) − ∂ r rνρ ( ∂ r u z + ∂ z u r ) = 0 , (B.12a) ∂ t ( rρu r ) − ∂ z rνρ ( ∂ z u r + ∂ r u z ) − ∂ r rνρ ( ∂ r u r − ∂ z u z − u r r ) = 0 , (B.12b)39here the terms with the order of O ( δ t M a ) have been neglected, and thekinematic viscosity can be given by ν = c s δ t ( τ g −
12 ) , (B.13)where 1 /τ g = s g = s g = s g . Substituting Eq. (B.6a) into Eqs. (B.12) andfurther combining Eqs.(B.6a-B.6c) and Eqs. (B.12a-B.12b) at t and t timescales, we can obtain the hydrodynamic equations, ∂ α ( ru α ) = 0 , (B.14a) ∂ t ( rρu β ) + ∂ α ( rρu α u β ) = − ∂ β ( rp ) + ∂ α [ rνρ ( ∂ α u β + ∂ β u α )] + F β . (B.14b)From the above discussions, it is shown that our MRT-LB model can exactlyrecover the axisymmetric hydrodynamic equations under the incompressiblelimits. Appendix C. Calculation of the hydrodynamic pressure
Lastly, we give a discussion on the calculation of the hydrodynamic pres-sure p . Based on Eq. (25), we have g ( eq )0 ( x , t ) = r ( ω − c s p ( x , t ) + rρs ( u ( x , t )) , (C.1)which can be further recast as g ( x , t ) − [ g ( x , t ) − g ( eq )0 ( x , t )] = r ( ω − c s p ( x , t ) + rρs ( u ( x , t )) . (C.2)Expanding g i ( x + c i δ t , t + δ t ) in Eq. (24) about x and t , we can easily derive δ t D i g i ( x , t ) = − (cid:2) M − S g M (cid:3) ij (cid:2) g j ( x , t ) − g eqj ( x , t ) (cid:3) + δ t G i ( x , t ) (C.3)40o the order of O ( δ t ). Premultiplying the matrix M − S g − M on both sidesof Eq. (C.3), we can get g i ( x , t ) − g ( eq ) i ( x , t ) = − δ t ( M − S g − M ) ij D j g j ( x , t )+ δ t ( M − S g − M ) ij G j ( x , t ) , (C.4)which indicates g i ( x , t ) = g ( eq ) i ( x , t ) + O ( δ t ) . (C.5)With the above approximate, we can rewrite Eq. (C.4) as g i ( x , t ) − g ( eq ) i ( x , t ) = − δ t ( M − S g − M ) ij D j g ( eq ) j ( x , t )+ δ t ( M − S g − M ) ij G j ( x , t ) . (C.6)Substituting Eqs. (36) and (37) into Eq. (C.6), and taking i = 0, we canhave g ( x , t ) − g ( eq )0 ( x , t ) = ( 3 s g + 2 s g c s g s g ) δ t ∂rp∂t +( s g + s g c s g s g ) δ t ∂ ( rρ u · u ) ∂t − s g + s g − s g s g s g s g · δ t ω ρu r (C.7)to the order of O ( δ t ). Under the incompressible condition, ∂ t p = O ( M a )and | u | = O ( M a ) are satisfied, and then we can rewrite Eq. (C.7) as g ( x , t ) − g ( eq )0 ( x , t ) = − s g + s g − s g s g s g s g · δ t ω ρu r + O ( δ t + δ t M a ) . (C.8)Substituting Eq. (C.8) into Eq. (C.2) yields r ( ω − c s p = g − rρs ( u ) + s g + s g − s g s g s g s g · δ t ω ρu r (C.9)to the order of O ( δ t + δ t M a ). Considering the first component in Eq. (B.4),the zeroth-order moment of the distribution function can be defined by X i g i = − δ t r u · ∇ ρ, (C.10)41nd the hydrodynamic pressure p can be evaluated as p = c s − ω " r X i =0 g i + δ t u · ∇ ρ + ρs ( u ) − s g + s g − s g s g s g s g · δ t ω ρu r r . (C.11) References [1] C. Hirt, B. Nichols, Volume of fluid (VOF) method for the dynamics offree boundaries, J. Comput. Phys. 39 (1981) 201-225.[2] S.O. Unverdi, G. Tryggvason, A front-tracking method for viscous, in-compressible, multi-fluid flows, J. Comput. Phys. 100 (1992) 25-37.[3] M. Sussman, P. Smereka, S. Osher, A level set approach for computingsolutions to incompressible two-phase flow, J. Comput. Phys. 114 (1994)146-159.[4] D. M. Anderson, G. B. McFadden, and A. A. Wheeler, Diffuse-interfacemethods in fluid mechanics, Annu. Rev. Fluid Mech. 30 (1998) 139-165.[5] H. Ding, P. D. M. Spelt, and C. Shu, Diffuse interface model for incom-pressible two-phase flows with large density ratios, J. Comput. Phys. 226(2007) 2078-2095.[6] D. Jacqmin, Calculation of two-phase Navier-Stokes flows using phase-field modeling, J. Comput. Phys. 155 (1999) 96-127.[7] Z. L. Guo, C. Shu, Lattice Boltzmann Method and its Applications inEngineering, World Scientific, 2013.428] T. Kr¨ u ger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, and E. M.Viggen, The lattice Boltzmann Method: Principles and Practice, Springer,2016.[9] J. W. Cahn, J. E. Hilliard, Free energy of a nonuniform system. I. Inter-facial free energy, J. Chem. Phys. 28 (1958) 258-267.[10] Y. Sun and C. Beckermann, Sharp interface tracking using the phase-field equation, J. Comput. Phys. 220 (2007) 626-653.[11] P. H. Chiu and Y. T. Lin, A conservative phase field method for solvingincompressible two-phase flows, J. Comput. Phys. 230 (2011) 185-204.[12] A. K. Gunstensen, D. H. Rothman, S. Zaleski, and G. Zanetti, LatticeBoltzmann model of immiscible fluids, Phys. Rev. A 43 (1991) 4320-4327.[13] X. Shan and H. Chen, Lattice Boltzmann model for simulating flowswith multiple phases and components, Phys. Rev. E 47 (1993) 1815-1819.[14] M. Swift, W. Osborn, and J. Yeomans, Lattice Boltzmann Simulationof Nonideal Fluids, Phys. Rev. Lett. 75 (1995) 830-833.[15] X. He, S. Chen, and R. Zhang, A lattice Boltzmann scheme for incom-pressible multiphase flow and its application in simulation of Rayleigh-Taylor instability, J. Comput. Phys. 152 (1999) 642-663.[16] T. Lee and C. L. Lin, A stable discretization of the lattice Boltzmannequation for simulation of incompressible two-phase flows at high densityratio, J. Comput. Phys. 206 (2005) 16-47.4317] H. Liang, B. C. Shi, Z. L. Guo, and Z. H. Chai, Phase-field-basedmultiple-relaxation-time lattice Boltzmann model for incompressible mul-tiphase flows, Phys. Rev. E 89 (2014) 053320.[18] H. Liang, B. C. Shi, and Z. H. Chai, Lattice Boltzmann modeling ofthree-phase incompressible flows, Phys. Rev. E 93 (2016) 013308.[19] H. Liu, Q. J. Kang, C. R. Leonardi, S. Schmieschek, A. Narv´ a ez, B. D.Jones, J. R. Williams, A. J. Valocchi, and J. Harting, Multiphase latticeBoltzmann simulations for porous media applications, Comput. Geosci. 20(2016) 777-805.[20] Q. Li, H. K. Luo, Q. J. Kang, Y. L. He, Q. Chen, and Q. Liu, LatticeBoltzmann methods for multiphase flow and phasechange heat transfer,Prog. Energy Combust. Sci. 52 (2016) 62-105.[21] M. R. Nobari, Y. J. Jan, G. Tryggvason, Head-on collision of drops-Anumerical investigation, Phys. Fluids 8 (1996) 29-42.[22] K. N. Premnath, J. Abraham, Simulations of binary drop collisionswith a multiple-relaxation-time lattice-Boltzmann model, Phys. Fluids 17(2005) 122105.[23] J. Hua, J. Lou, Numerical simulation of bubble rising in viscous liquid,J. Comput. Phys. 222 (2007) 769-795.[24] H. B. Huang, J. J. Huang, X. Y. Lu, A mass-conserving axisymmetricmultiphase lattice Boltzmann method and its application in simulation ofbubble rising, J. Comput. Phys. 269 (2014) 386-402.4425] D. R. Link, S. L. Anna, D. A. Weitz, and H. A. Stone, Geometricallymediated breakup of drops in microfluidic devices, Phys. Rev. Lett. 92(2004) 054503.[26] A. S. Utada, A. Fernandez-Nieves, H. A. Stone, and D. A. Weitz, Drip-ping to jetting transitions in coflowing liquid streams, Phys. Rev. Lett. 99(2007) 094502.[27] M. Bussmann, S. Chandra, and J. Mostaghimi, Modeling the splash ofa droplet impacting a solid surface, Phys. Fluids 12 (2000) 3121-3132.[28] C. Josseranda and S. Zaleskib, Droplet splashing on a thin liquid film,Phys. Fluids 15 (2003) 1650-1657.[29] I. Halliday, L. A. Hammond, C. M. Care, K. Good, and A. Stevens, Lat-tice Boltzmann equation hydrodynamics, Phys. Rev. E 64 (2001) 011208.[30] Y. Peng, C. Shu, Y.T. Chew, J. Qiu, Numerical investigation of flows inCzochralski crystal growth by an axisymmetric lattice Boltzmann method,J. Comput. Phys. 186 (2003) 295-307.[31] Z. L. Guo, H. F. Han, B. C. Shi, and C. G. Zheng, Theory of the latticeBoltzmann equation: Lattice Boltzmann model for axisymmetric flows,Phys. Rev. E 79 (2009) 046708.[32] Q. Li, Y. L. He, G. H. Tang, and W. Q. Tao, Improved axisymmetriclattice Boltzmann scheme, Phys. Rev. E 81 (2010) 056707.[33] L. Zheng, Z. L. Guo, B. C. Shi, C. G. Zheng, Kinetic theory based45attice Boltzmann equation with viscous dissipation and pressure work foraxisymmetric thermal flows, J. Comput. Phys. 229 (2010) 5843-5856.[34] L. Zhang, S. Yang, Z. Zeng, J. Chen, L. Wang, J.W. Chew, Alternativeextrapolation-based symmetry boundary implementations for the axisym-metric lattice Boltzmann method, Phys. Rev. E 95 (2017) 043312.[35] K. N. Premnath, J. Abraham, Lattice Boltzmann model for axisymmet-ric multiphase flows, Phys. Rev. E 71 (2005) 056706.[36] S. Mukherjee, J. Abraham, Lattice Boltzmann simulations of two-phaseflow with high density ratio in axially symmetric geometry, Phys. Rev. E75 (2007) 026701.[37] Q. Lou, Z. L. Guo, and B. C. Shi, Effects of force discretization on massconservation in lattice Boltzmann equation for two-phase flows, Europhys.Lett. 99 (2012) 64005.[38] J. J. Huang, H. B. Huang, C. Shu, Y. T. Chew, and S. L. Wang, Hybridmultiple-relaxation-time lattice-Boltzmann finite-difference method for ax-isymmetric multiphase flows, J. Phys. A: Math. Theor. 46 (2003) 055501.[39] S. Srivastava, P. Perlekar, J. H. M. ten Thije Boonkkamp, N. Verma,and F. Toschi, Axisymmetric multiphase lattice Boltzmann method, Phys.Rev. E 88 (2013) 013309.[40] H. Liang, Z. H. Chai, B. C. Shi, Z. L. Guo, and T. Zhang, Phase-field-based lattice Boltzmann model for axisymmetric multiphase flows, Phys.Rev. E 90 (2014) 063311. 4641] H. Liu, L. Wu, Y. Ba, G. Xi, Y. Zhang. A lattice Boltzmann method foraxisymmetric multicomponent flows with high viscosity ratio, J. Comput.Phys. 327 (2016) 873-893.[42] H. L. Wang, Z. H. Chai, B. C. Shi, H. Liang, Comparative study of thelattice Boltzmann models for Allen-Cahn and Cahn-Hilliard equations,Phys. Rev. E 94 (2016) 033304.[43] A. Fakhari and D. Bolster, Diffuse interface modeling of three-phasecontact line dynamics on curved boundaries: A lattice Boltzmann modelfor large density and viscosity ratios, J. Comput. Phys. 334 (2017) 620-638.[44] Z. H. Chai, D. K Sun, H. L. Wang, B. C. Shi, A comparative study oflocal and nonlocal Allen-Cahn equations with mass conservation, Int. J.Heat Mass Tran. 122 (2018) 631-642.[45] H. Liang, J. R. Xu, J. X. Chen, H. L. Wang, Z. H. Chai, and B. C.Shi, Phase-field-based lattice Boltzmann modeling of large-density-ratiotwo-phase flows, Phys. Rev. E 97 (2018) 033309.[46] F. Ren, B. Song, M. C. Sukop, and H. Hu, Improved lattice Boltzmannmodeling of binary flow based on the conservative Allen-Cahn equation,Phys. Rev. E 94 (2016) 023311.[47] J. Kim, A continuous surface tension force formulation for diffuse-interface models, J. Comput. Phys. 204 (2005) 784-804.[48] P. Lallemand and L. S. Luo, Theory of the lattice Boltzmann method:Dispersion, dissipation, isotropy, Galilean invariance, and stability, Phys.Rev. E 61 (2000) 6546-6562. 4749] I. Ginzburg, F. Verhaeghe, and D. d’Humi` e res, Two-relaxation-timelattice Boltzmann scheme: About parametrization, velocity, pressure andmixed boundary conditions, Commun. Comput. Phys. 3 (2008) 427-478.[50] Y. H. Qian, D. d’Humi` ee