On the formulations of interfacial force in the phase-field-based lattice Boltzmann method
AAPS/123-QED
On the formulations of interfacial force in the phase-field-basedlattice Boltzmann method
Chunhua Zhang and Zhaoli Guo ∗ State Key Laboratory of Coal Combustion,Huazhong University of Science and Technology, Wuhan 430074, China
Hong Liang
Department of Physics, Hangzhou Dianzi University, Hangzhou 310018, China (Dated: January 1, 2021)
Abstract
Different formulations of interfacial force have been adopted in phase-field-based lattice Boltz-mann method for two-phase flows. Although they are identical mathematically, their numericalperformances may be different due to truncation errors in the discretization. In this paper, four-type formulations of interfacial force available in the literature, namely stress tensor form (STF),chemical potential form (CPF), pressure form (PF) and continuum surface force (CSF) form, arecompared and discussed. A series of benchmark problems, including stationary droplet, two merg-ing droplets, Capillary wave, rising bubble and drop deformation in shear flow, are simulated.Numerical results show that CPF is a good choice for small surface deformation problems whileSTF is preferred for dynamical problems, both STF and CSF demonstrate good numerical stability. ∗ [email protected] a r X i v : . [ phy s i c s . f l u - dyn ] D ec . INTRODUCTION Multiphase flows are ubiquitous in both natural processes and industrial applications,such as droplet dynamics [1], lab-on-chip devices [2], surfactant behavior [3], undergroundwater flows [4] and enhanced oil recovery [5]. A number of numerical methods have beendeveloped for simulating such flows, which can be divided into two categories, i.e, interfacetracking approach and interface capturing approach. In the former, interfaces are explicitlytracked, such as the marker and cell method [6] and front-tracking method [7]. In the latter,interfaces are implicitly tracked and an interface function that marks the location of theinterface is governed by the advection (diffusion) equations, such as volume of fluid (VOF)method [8], level set (LS) method [9] and phase field method [10].Among these methods, the phase field method is an increasingly popular choice for mul-tiphase fluids simulations. The basic idea is to introduce a so-called order parameter thathas distinct values in the bulk phases but varies smoothly over the interfacial region. Theorder parameter defined as the volume fraction or mass fraction is usually governed by thephase field equations, such as the Cahn-Hilliard equation or the Allen-Cahn equation, whichleads to the Navier-Stokes-Cahn-Hilliard (NSCH) system or the Navier-Stokes-Allen-Cahn(NSAC) system. If the fluid density is taken as an order parameter, the flow can be de-scribed by the Navier-Stokes-Kortweg (NSK) system [11, 12]. Although the momentumequations with interfacial force in NSCH, NSAC and NSK are very similar, the propertiesof these equations are different. In the NSCH and NSAC equations, the pressure servesas an auxiliary variable associated with the incompressibility (or quasi-incompressibility)condition. In the NSK equations, the pressure is connected to the density via an equa-tion of state. In the traditional computational fluid dynamics (CFD), many discretizationmethods have been developed to numerically solve the above governing equations. Recently,the lattice Boltzmann method (LBM) has grown as an alternative tool for multiphase flowsimulations [13–15]. The LBM is a mesoscopic method based on certain kinetic models. InLBM, the fluid is represented by a discrete set of particle distribution functions which onlyperform propagation and collision processes on a fixed lattice. The macroscopic quantitiesof the flow are calculated by taking the moments of the particle distribution functions. LBMis simple and easy to be implemented compared with the traditional CFD to discretize the2acroscopic governing equations. However, it can be shown that the corresponding phasefield equation and hydrodynamic equations can be recovered from the lattice Boltzmannequations through the Chapman-Enskog (CE) analysis.In computational methods for multiphase flows, approximating the surface tension forceaccurately is critical to capture correct flow behaviors. A number of mathematica modelsfor the interfacial force are available in phase-field-based lattice Boltzmann methods so far.In fact, the interfacial force can be strictly derived based on the entropy principle of rationalthermodynamics [10, 16–18]. The resulting interfacial force appears as a gradient of the stresstensor of the order parameter in the modified momentum equation. These formulations canbe called stress tensor form (STF). The stress form can be further simplified by redefiningthe pressure. Then, the interfacial force can be expressed as the forms dependent on thegradients of the order parameter [19–23]. These formulations can be called pressure form(PF). If the chemical potential related to the order parameter is employed, the interfacialforce can also be expressed as the forms dependent on the chemical potential [24–27], whichcan be named as chemical potential form (CPF). Mathematically, the STF, PF and CPF areequivalent. In addition, in the continuum surface force (CSF) model of Brackbill et al. [28],the interfacial force is treated as a volumetric force proportional to the normal vector andcurvature of the interface and a surface Dirac function localizing the interfacial force to theinterface, which has been widely used in the VOF and LS methods. Based on the CSFmodel, Kim et al. [29] proposed a CSF type interfacial force for phase field methods. Thebasic idea is to replace the level set by the order parameter and take the square of gradient ofthe order parameter as the surface Dirac function. An advantage of the CSF formulation isthat the pressure field can be calculated explicitly while the calculated pressure field with theprevious interfacial forces includes some gradient terms of the order parameter except thetrue pressure. The surface Dirac function in CSF model can also be defined in other ways.For instance, Lee and Kim et al. [30] compared various types of surface Dirac functions inthe CSF model. They argued that the absolute value of the gradient of the order parameterhas the best performances in their considered numerical experiments. These formulationsare called CSF form of the interfacial force in the present work. It’s worth noting that thecalculation of the normal vectors and the curvature at the interface is critical in the CSFmodels.Although most of the above interfacial force formulations are mathematically equivalent,the performance of each formulation may be different in practical computations. For ex-3mple, Lee and Fischer et al. [19] compared the parasitic currents between the pressureform and potential form in LBM, and the results showed that potential form yielded muchsmaller parasitic currents. Chao and Mei et al. [31] compared the interface force distribu-tion between the pressure form and the CSF form, and the results showed that the pressureform could generate wiggles over the interface region while the CSF form produced no suchunphysical results. However, there is a lack of systematic study of the performance of thesefour interfacial force formulatiions widely used in LBM, and this paper will focus on thistopic.The paper is organized as follows. In section 2, the governing equations of the phasefield model for binary fluids are presented and the formulas of surface tension force aresummarized. The phase-field-based lattice Boltzmann method is briefly introduced in section3. In section 4, several benchmark problems are investigated and the results are compared.Finally, conclusions are drawn in Section 5.
II. MATHEMATICAL FORMULATIONA. Governing equations
In this study, we consider the NSCH equations for multiphase flows. The Cahn-Hillardequation is expressed as [10, 32] ∂φ∂t + ∇ · ( φ u ) = ∇ · M ∇ µ φ , (1)where φ is the order parameter to identify different phases, M is the mobility, µ φ is thechemical potential that is defined as µ φ = δψδφ = ∂f ∂φ − κ ∇ φ, (2)where ψ is the system free energy, ψ = (cid:90) V (cid:104) f ( φ ) + κ |∇ φ | (cid:105) dV, (3)where f = β (1 − φ ) is the bulk energy density, the second term is the interface energydensity, β and κ are determined by the surface tension σ and the interface width W .For a plane interface at equilibrium, the equilibrium profile for the order parameter canbe obtained by solving µ φ = 0, φ ( r ) = tanh (cid:32)(cid:114) βκ r (cid:33) , (4)4here r is the signed distance function which is the coordinate normal to the interface. (cid:112) κ/ β has a length scale of interface thickness. As the surface tension is interpreted asenergy per unit surface area, the surface tension for a flat interface with equilibrium profilecan be calculated by σ = (cid:90) + ∞−∞ (cid:16) f ( φ ) + κ |∇ φ | (cid:17) dr = κ (cid:90) + ∞−∞ |∇ φ | dx = 43 (cid:112) βκ, (5)In Ref.[10], (cid:112) κ/ β is defined as W/
2, which leads to β = 34 σW , κ = 38 W σ. (6)The dynamics of a fluid mixture of two incompressible viscous fluids can be described bythe Navier-Stokes equations with interfacial force [10, 26] ∇ · u = 0 , (7) ∂ ( ρ u ) ∂t + ∇ · ( ρ uu ) = −∇ P sf + ∇ · µ ( ∇ u + ∇ u T ) + F g + F sf , (8)where ρ is the fluid density, u is the flow velocity, P sf is the generalized pressure dependenton the definition of the interfacial force, µ is the dynamic viscosity, F g = ( ρ − ρ ) g is thegravitational force with g being the gravitational acceleration and ρ being the backgrounddensity, F sf is the interfacial force. The subscript sf (= stf, cpf, pf, csf ) denotes differentformulations of interfacial force.The mixture density ρ and viscosity µ can be given by ρ = ρ φ ρ − φ µ = µ φ µ − φ u (cid:48) = u U c , x (cid:48) = x L c , t (cid:48) = tT c , p (cid:48) = P sf p c , µ (cid:48) φ = µ φ µ φ,c , F (cid:48) sf = F sf L c σ , (11)where U c , L c , T c (= L c /U c ) , p c (= ρ c U c ) , µ φ,c (= 4 β ) are respectively the reference velocity,length, time, pressure and chemical potential. In this paper, the density and dynamicalviscosity of fluid 1 are chosen as the reference quantities, i.e, ρ c = ρ , µ c = µ . With the5bove variables and dropping the primes, the dimensionless governing equations can bewritten as ∂ t φ + ∇ · ( φ u ) = 1Pe ∇ · ( M ∇ µ φ ) , (12) ∂ t ( ρ u ) + ∇ · ( ρ uu ) = −∇ P sf + 1Re ∇ · µ ( ∇ u + ∇ u T ) + 1We F sf + 1Fr F g , (13) ∇ · u = 0 , (14)with µ φ = φ ( φ − − Cn ∇ φ,ρ = 1 + φ − φ ρ ρ ,µ = 1 + φ − φ µ µ . (15)The dimensionless groups used above are the Reynolds number Re, Peclet number Pe, Webernumber We, Frounde number Fr and Cahn number Cn, which are respectively defined byRe = ρ c U c L c µ c , Pe = U c L c M β ,
We = ρ c L c U c σ , Fr = U c √ gL c , Cn = WL c , (16) B. Interfacial force formulations
Based on the energetic variational approach or the free energy inequality, the surfacetension force in the momentum equation can be defined as [17, 33, 34] F stf − = −∇ · κ ( ∇ φ ⊗ ∇ φ ) , (17)where ∇ φ ⊗ ∇ φ is the usual tensor product and denotes the induced elastic stress due to themixing of the different species. In this case, the generalized pressure P sf in Eq. (8) includesboth the hydrostatic pressure p h due to the incompressibility and the contributions fromthe induced stress, P stf − = p h + κ |∇ φ | . In Ref. [35, 36], the surface tension force term isdefined as F stf − = ∇ · κ ( |∇ φ | I − ∇ φ ⊗ ∇ φ ) , (18)which implies that the principle axes of the tensor are perpendicular to the tangent planeof the interface. The normal stress perpendicular to the tangent plane of the interface iszero and the two tangent normal stresses are equal. In this case, the generalized pressure inEq. (8) becomes the true pressure, namely, P stf − = p h [10, 36].6or simplicity, we assume that the surface tension σ is constant. By using the followingidentity κ ∇ · ( ∇ φ ⊗ ∇ φ ) = κ ∇|∇ φ | + κ ∇ φ ∆ φ = ∇ (cid:16) κ |∇ φ | + κφ ∆ φ (cid:17) − κφ ∇ ∆ φ = ∇ (cid:16) κ |∇ φ | + f (cid:17) − µ φ ∇ φ = ∇ (cid:16) κ |∇ φ | + f − φµ φ (cid:17) + φ ∇ µ φ , (19)and absorbing the gradient terms into pressure p h , the surface tension force can be expressedas F cpf − = − φ ∇ µ φ , F cpf − = µ φ ∇ φ, F pf − = − κ ∇ φ ∆ φ, F pf − = κφ ∇ ∆ φ. (20)The corresponding generalized pressure is redefined as P cpf − = p h + f − φµ φ − κ |∇ φ | ,P cpf − = p h + f − κ |∇ φ | ,P pf − = p h − κ |∇ φ | ,P pf − = p h + κφ ∆ φ − κ |∇ φ | . (21) F cpf − and F cpf − are termed as chemical potential form. F pf − and F pf − are the pressureform. It is noted that F stf − is used in [37, 38] and F stf − is used in [35, 39–41], F cpf − isused in [24, 42, 43] and F cpf − is used in [13, 25–27], F pf − is used in [22, 23] and F pf − isused in [20, 21].Based on the CSF model, the surface tension force can be given by [28, 44] F csf = σ (cid:101) κδ s n , (22)where n is the unit normal vector, (cid:101) κ = −∇ · n is the local mean curvature, δ s is the surfaceDirac function used to ensure the force acting on the interfacial region. To match the surfacetension of the sharp interface model, the Dirac function should satisfy (cid:90) ∞−∞ δ s dr = 1 . (23)There are many possible choices for δ s . Kim [29] proposed to use α |∇ φ | as the Diracfunction with α = 3 W/ F csf − = − κ ∇ φ |∇ φ |∇ · n . (24)7ee and Kim et.al [30] proposed α |∇ φ | as the Dirac function with α = 0 . F csf − = − σ ∇ φ ( ∇ · n ) . (25)The derivation of α is referred to Appendix.B. In above interfacial force formulations,Eqs. (17), (18) and (20) are identical mathematically. In fact, these formulations can berewritten as F stf − = F csf − − (cid:20) ∇ κ |∇ φ | κ ∇ φ ( ∇ φ · ∇|∇ φ | ) |∇ φ | (cid:21) , F stf − = F csf − − (cid:20) −∇ κ |∇ φ | κ ∇ φ ( ∇ φ · ∇|∇ φ | ) |∇ φ | (cid:21) , F cpf − = F csf − − (cid:20) ∇ ( φµ φ ) − ∇ f + κ ∇ φ ( ∇ φ · ∇|∇ φ | ) |∇ φ | (cid:21) , F cpf − = F csf − − (cid:20) −∇ f + κ ∇ φ ( ∇ φ · ∇|∇ φ | ) |∇ φ | (cid:21) , F pf − = F csf − − κ ∇ φ ( ∇ φ · ∇|∇ φ | ) |∇ φ | , F pf − = F csf − − (cid:20) −∇ ( κφ ∆ φ ) + κ ∇ φ ( ∇ φ · ∇|∇ φ | ) |∇ φ | (cid:21) . (26)It is clear that there are some gradient terms in F stf , F cpf and F pf . This is why the previousformulations cannot be used to calculate the pressure field explicitly [29].By using Eq. (4), the following relations can be obtained |∇ φ | = 2 W (1 − φ ) , ∇ φ ( ∇ φ · ∇|∇ φ | ) |∇ φ | = 12 ∇|∇ φ | . (27)Inserting Eq. (27) into Eq. (26) leads to F stf − = F csf − − ∇ κ |∇ φ | , F cpf − = F csf − − ∇ ( φµ φ ) , F pf − = F csf − − ∇ κ |∇ φ | , F pf − = F csf − − ∇ ( − κφ ∆ φ + κ |∇ φ | ) , F stf − = F cpf − = F csf − . (28)Therefore, F csf − , F stf − and F cpf − are identical when the system is at equilibrium. Themain difference between F csf − and F csf − is the definition of the Dirac delta function. Allabove formulations have been used to mimic the interfacial force in the phase-field-basedLBM. In Sec. V, the performance of the LBM models with the above eight formulations ofsurface tension force will be compared. 8 II. PHASE-FIELD-BASED LATTICE BOLTZMANN METHOD
We adopted the multiphase LBM of He et al [20] for the hydrodynamic equations andthe improved LBM of Zhang et al [45] for Cahn-Hilliard equation. The evolutions of thedistribution functions f i and h i are respectively expressed as f i ( x + c i δt, t + δt ) − f i ( x , t ) = − τ f [ f i ( x , t ) − f eqi ( x , t )] + δt (cid:18) − τ f (cid:19) F i , (29) h i ( x + c i δt, t + δt ) − h i ( x , t ) = − τ h [ h i ( x , t ) − h eqi ( x , t )] + δt (cid:18) − τ h (cid:19) H i , (30)where f i ( x , t ) and h i ( x , t ) are the distribution functions for the hydrodynamics and orderparameter fields respectively, c i is the discrete velocity in the i-th direction, δt is the timestep, τ f and τ h are the dimensionless relaxation times related to the shear viscosity andmobility respectively, F i and H i are the discrete force terms. To recover the correct governingequations, the equilibrium distributions f eqi and h eqi are defined as f eqi = ω i [ P sf + c s ρs i ( u )] (31) h eqi = φ + ( ω − ηµ φ , i = 0 ω i ηµ φ + ω i φ c i · u c s , i (cid:54) = 0 (32)with s i ( u ) = c i · u c s + uu : ( c i c i − c s I )2 c s , (33)where ω i is the weighting coefficient corresponding to the discrete velocity c i , c s = c/ √ c = δx/δt is the lattice speed with δx being the lattice lengthscale, and η is an adjustable parameter for the mobility. In this work, the two-dimensionalnine-velocity (D2Q9) model is used in which the discrete velocity is( c , c , c , c , c , c , c , c , c )= c − − − − − − (34)and the corresponding weighting coefficients are ω = 4 / ω − = 1 / ω − = 1 / F i and H i are given by F i = ( c i − u ) · (cid:2) Γ i ( u )( F sf + F g ) + s i ( u ) ∇ c s ρ (cid:3) (35)9 i = ¯ ω i τ τ h τ δt ∇ · ( u φ ) + ω i c i · ∂ t ( φ u ) c s . (36)where Γ( u ) = ω i + s i ( u ), τ = ( τ h − τ h + ) δt , τ = ( τ h − . δt and ¯ ω i is a new weightcoefficient and satisfies (cid:80) i ¯ ω i = (cid:80) i ¯ ω i c i = 0 , (cid:80) i ¯ ω i c i c i = c s I . In particular, if τ g = 0 . √ /
6, the above improved LBM for CH can be simplified to the one of Liang et al [13].However, the relationship between the Peclet number and Cahn number should be Pe ∼ Cn − to achieve the sharp-interface limit with continuous mesh refinement [46]. Then,the relaxation time may have a value except the optimum one in some situations and theimproved LBM should be considered.The macroscopic quantities are calculated by P sf = (cid:88) i f i + δt c s u · ∇ ρ, u = 1 c s ρ (cid:34)(cid:88) i c i f i + δt c s ( F sf + F g ) (cid:35) φ = (cid:88) i h i , (37)Through the Chapman-Enskog expansion, the macroscopic governing equations recoveredfrom the above LBM are1 c s ρ ∂ t P sf + ∇ · u = 0 , (38) ∂ t ( ρ u ) + ∇ · ( ρ uu ) = −∇ P sf + ∇ · µ ( ∇ u + ∇ u T ) + F sf + F g , (39) ∂ t φ + ∇ · ( φ u ) = ∇ · M ∇ µ φ , (40)where the viscosity µ and the mobility M are defined as µ = ρc s ( τ f − . δt and M = c s η ( τ h − . δt , respectively.The gradient terms in each formulation of interfacial force can be calculated with differentschemes. In the present work, we will use the isotropic central scheme [47], ∇ Ψ = 1 c s δt (cid:88) i =1 ω i c i Ψ( x + c i δt ) , (41) ∇ Ψ = 2 c s δt (cid:88) i =1 ω i [Ψ( x + c i δt ) − Ψ( x )] , (42)where Ψ denotes arbitrary quantity. For a node located at wall boundary, a second-orderone-side finite difference is employed. 10 V. BOUNDARY CONDITIONS
Boundary treatment is one of the most important tasks in numerical methods. In LBM,the classical boundary condition to model walls is the bounce-back method, which can berealised by both the full-way bounce-back and the half-way bounce-back [48]. As the half-way bounce-back can be implemented without solid nodes and is more accurate for unsteadyflows, we will only consider the half-way bounce-back in the practical calculation. As shownin Fig 1, following Ladd’s half-way bounce-back scheme, the unknown distribution functionis determined by [14, 49] f ¯ i ( x f , t + δt ) = f + i ( x f , t ) − ω i ρ ( x w , t ) c i · u w ,g ¯ i ( x f , t + δt ) = g + i ( x f , t ) − ω i φ ( x w , t ) c i · u w c s , (43)where f ¯ i and g ¯ i are the distribution function with the velocity c ¯ i = − c i , the superscript’+’ denotes the post-collision value of the corresponding distribution function and u w is theprescribed wall velocity. For a stationary boundary with u w = 0, the above equations canbe used for the non-slip boundary.For the order parameter, the following boundary conditiions are employed, n w · ∇ φ = 0 , n w · ∇ µ φ = 0 , (44)where n w is the unit outward normal defined at the solid boundary. Eq.(44) means that theorder parameter conserves mass over the entire domain. In addition, the density ρ ( x w , t )can be approximated by ρ ( x f , t ). Here we use ∇ φ · n w = 0 to interpolate the density at thewall. V. NUMERICAL RESULTS AND DISCUSSION
In this section, the performance of each interfacial force formulation is validated by a seriesof benchmark tests, including stationary droplet, two merging droplets, capillary wave, risingbubble and the deformation droplet in a shear flow. For each test, the results obtained bythe lattice Boltzmann equation (LBE) model with different interfacial force formulations arecompared with the theoretical solutions or the available reference solutions in the literature.In Eq. (36), the time derivative is calculated by explicit Euler scheme, and ¯ ω = ω − ω i = ω i for i >
0. The Peclet number is set to be 1 . / Cn and the interface width is set tobe four grids unless otherwise stated. 11ig. 1: Illustration for the half-way bounce-back. The thin solid straight line is the gridline and the dashed line corresponds to the computational boundary. The black circles arethe fluid nodes and the black square is the solid node. The arrow represents the particle’sdirection, the rightmost grey shaded domain is the solid region.
A. Stationary droplet
We first make a comparison among different interfacial force formulations by simulating astationary droplet. Theoretically, the exact solution is zero velocity for all time. Initially, acircle droplet with radius R is placed at the center of the domain L × L . The order parameteris set to be φ ( x, y ) = tanh (cid:32) R − (cid:112) ( x − x c ) + ( y − y c ) W (cid:33) , (45)where ( x c , y c ) is the center coordinate of the droplet. Periodic boundary conditions areapplied to all the boundaries. The initial velocity field is set to be zero. The physicalparameters are set to be L = 1m, R = 0 . ρ =4kg/ m , ρ =1 kg/ m , ν = ν = 0 . m /s and σ = 0 . × , × , ×
240 are used. Thecharacteristic velocity is U c = σ/µ .We first examine the shape of the droplet at equilibrium. The interface profile of thedroplet obtained by all interfacial force formulations are similar and agree well with theinitial interface profile, and the results are not shown here. It is also found that the deviationbetween the numerical results given by all formulations and the analytical interface profilebecomes small as the value of mobility decreases, which is also consistent with the results in[13]. Since the definition of characteristic velocity is artificial to some extent, the relationshipof Pe ∼ / Cn may be unable to produce the closest results to the exact one.From the Laplace law, the numerical surface tension can be calculated by σ num = R num × ( p in − p out ). The relative error, Err = | σ num − σ exact | /σ exact × F csf − gives the smallest error while F pf − gives the largest one .TABLE I: Comparison of numerical surface tension based on Laplace law ( σ = 0 . F sf ×
60 120 ×
120 240 × σ num Err (%) σ num Err (%) σ num Err (%) F stf − F stf − F cpf − F cpf − F pf − F pf − F csf − F csf − The pressure field p h on 240 ×
240 meshes is presented in Figure 2. It can be seen thatthe pressure inside the droplet is higher than that in the surrounding fluid. However, F stf − , F csf − and F csf − give smooth pressure field across the interface while the others give obviousoscillation near the interface.The spurious velocity for each formulation is also examined. The magnitudes of spuriousvelocities denoted by Ca = µ | u max | /σ are presented in Table II. It can be seen that both F cpf − and F cpf − give small spurious velocities while the others give larger ones.Finally, the absolute values of interfacial force across the drop center with different for-mulations are compared. The results are shown in Figure 3. Theoretically, the interfacialforce should be zero everywhere except in the vicinity of the interface. However, the absolutevalues of F cpf − have non-zero values in the whole domain. This may cause earlier motion ofthe interface although the amplitude of interfacial force is small. Since the interface widthis fixed, the range of nonzero interfacial force decreases with increasing grid resolution. Inaddition, based on the definition of each formulation and the equilibrium state, the profileof the interfacial force should be symmetric with respect to the phase interface ( φ = 0).However, the interfacial force profiles of F stf − , F cpf − , F csf − , F csf − are symmetric whilethe profiles of the others are asymmetrical. 13 a) (b) (c) (d)(e) (f) (g) (h) Fig. 2: The pressure field on 240 ×
240 grid for (a) F stf − ,(b) F stf − ,(c) F cpf − ,(d) F cpf − ,(e) F pf − ,(f) F pf − ,(g) F csf − and (h) F csf − .TABLE II: The maximum spurious velocities of stationary droplet at equilibrium. F sf ×
60 120 ×
120 240 × | u | max Ca | u | max Ca | u | max Ca F stf − . × − . × − . × − . × − . × − . × − F stf − . × − . × − . × − . × − . × − . × − F cpf − . × − . × − . × − . × − . × − . × − F cpf − . × − . × − . × − . × − . × − . × − F pf − . × − . × − . × − . × − . × − . × − F pf − . × − . × − . × − . × − . × − . × − F csf − . × − . × − . × − . × − . × − . × − F csf − . × − . × − . × − . × − . × − . × − B. Droplets merging
To test the performance of the LBM with different interfacial force formulations, themerging of two droplets is simulated in this section. Initially, two circular droplets (density ρ d and viscosity ν d ) are placed in another fluid (density ρ s and viscosity ν s ) in a rectangledomain of L x × L y . When the initial gap d between two droplets is smaller than 2 W , merging14 x/R | F s tf - |
60 60120 120240 240 (a) x/R | F s tf - | -3
60 60120 120240 240 (b) x/R | F c p f - | -4
60 60120 120240 240 (c) x/R | F c p f - | -3
60 60120 120240 240 (d) x/R | F p f - |
60 60120 120240 240 (e) x/R | F p f - |
60 60120 120240 240 (f) x/R | F cs f - | -3
60 60120 120240 240 (g) x/R | F cs f - | -3
60 60120 120240 240 (h)
Fig. 3: The interfacial force profiles along the midline of the drop for (a) F stf − ,(b) F stf − ,(c) F cpf − ,(d) F cpf − ,(e) F pf − ,(f) F pf − ,(g) F csf − and (h) F csf − .will occur due to the surface tension effect. The order parameter is initialized to be φ ( x, y ) = 1+tanh (cid:32) R − (cid:112) ( x − x ) + ( y − y ) W (cid:33) +tanh (cid:32) R − (cid:112) ( x − x ) + ( y − y ) W (cid:33) , (46)where ( x , y ) = ( L x / − R − d/ , L y /
2) and ( x , y ) = ( L x / R + d/ , L y /
2) are thecenters of the two droplets, respectively. The initial velocity field is zero in the wholedomain. In simulations, the computational domain of L x × L y = 1 . ×
1m is discretizedby a uniform mesh 240 × R = R = 0 . d = 1 . W and W = 0 . ρ d = 5kg / m , ρ s = 1kg / m and the viscosities are ν d = ν s = 0 . / s. The surface tensioncoefficient is σ = 0 . / m, and the characteristic velocity is given by U c = (cid:112) σρ /R . ThePeclet number is set as Pe = 0 . / Cn. Periodic boundary conditions are implemented at allboundaries. With these parameters, merging will take place. Figure 4 shows the interfacialshapes of the droplets at t = 2T and 30T with T = (cid:112) ρ R /σ . The interfacial shapes at t = 30T are compared with analytical results. From Fig. 4, it is observed that the twodroplets gradually merge, oscillate and finally form a larger stationary droplet. Especially,the final interface shapes predicted by all formulations are in good agreement with theanalytical solutions. However, the interface positions predicted by the LBE models with15 cpf − , F pf − , F csf − at t = 2T are different from those of the other formulations. Thedroplets of the LBE models with F cpf − , F pf − and F csf − have started to merge while thedroplets with the other interfacial force formulations remain at a distinct distance. As noexternal forces are presented in the system, the mass centre of the droplets should not changeduring coalescence. Figure 5 shows the time development of the position of the mass centreof the droplets. All interfacial forces present similar accuracy. It’s worth pointing out thatthe computations with F cpf − , F cpf − , F pf − and F pf − become unstable when Pe = 1 / Cn.This implies that both F stf and F csf have a better numerical stability for this problem. x y (a) x y (b) x y (c) x y (d) x y (e) x y (f) x y (g) x y (h) Fig. 4: Interfaces of two droplets of equal sizes at t = 2 T (dotted line) and 30T(dashedline): (a) F stf − ,(b) F stf − ,(c) F cpf − ,(d) F cpf − ,(e) F pf − ,(f) F pf − ,(g) F csf − and (h) F csf − . Solid line represents the analytical solutions.We further simulate the above system but with two droplets of unequal sizes ( R =0 . , R = 0 . t = 2Tand 30T. The interface positions are different for each interfacial force formulation. Inparticular, the merged droplets predicted by the BE models with F cpf − , F cpf − and F csf − have a distinct movement. Figure 7 shows the time development of the position of the masscentre of the droplets, which shows that the positions predicted by F cpf ( F cpf − , F cpf − )and F csf ( F csf − , F csf − ) display significant deviations from their initial positions as timeincreases. 16 t/T x c F stf-1 F stf-2 F cpf-1 F cpf-2 F pf-1 F pf-2 F csf-1 F csf-2 Fig. 5: Time history of mass center x c for droplets of equal size. x y (a) x y (b) x y (c) x y (d) x y (e) x y (f) x y (g) x y (h) Fig. 6: Interfaces of two droplets of unequal sizes at t = 2T(dotted line) and 30T(dashedline): (a) F stf − ,(b) F stf − ,(c) F cpf − ,(d) F cpf − ,(e) F pf − ,(f) F pf − ,(g) F csf − and (h) F csf − . Solid line represents the analytical solutions. C. Capillary wave
We further test the numerical accuracy of the interfacial force formulations by a two-dimensional capillary wave problem. Initially, a heavier fluid is placed under a lighter fluidwith a small perturbation y = 1 . H + h cos( kx ) on the interface in a rectangle domain of H × H , where h is the initial perturbation amplitude and k = 2 π/H is the wave number.17 t/T x c F stf-1 F stf-2 F cpf-1 F cpf-2 F pf-1 F pf-2 F csf-1 F csf-2 Fig. 7: Time history of mass center x c for droplets of unequal size.The evolution of the interface wave amplitude h ( t ) is given by [50] h ( t ) h = 4(1 − β ) ν k − β ) ν k + ω erfc( √ νk t ) + (cid:88) i =1 z i Z i ω z i − νk e ( z i − νk ) t erfc( z i √ t ) , (47)where β = ρ ρ / ( ρ + ρ ) , ω = ( σk ) / ( ρ + ρ ), erfc( z i ) is the complementary error functionof a complex variable z i , z i ( i = 1 , . . . ,
4) are the four roots of the following algebraic equation z − β √ νk z + 2(1 − β ) νk z + 4(1 − β )( νk ) / z + (1 − β ) ν k + ω = 0 , (48)and Z i is defined as Z i = (cid:89) j (cid:54) = i ( z j − z i ) , i, j = 1 , · · · , . (49)In simulations, periodic boundaries are applied to the left and right sides and no-slip bound-aries are imposed on the top and bottom walls [51]. The physical parameters are set as H = 1m , ρ = ρ = 1kg/m , ν = ν = 0 . /s , σ = 0 . U c = (cid:112) σ/L c /ρ . Hence, the Reynolds number is Re = 50 and the Webernumber is We = 1. Two uniform grids of H = 80 and 160 are used. Figure 8 shows theevolution of the capillary amplitude for each grid. All the numerical results agree well withthe theoretical solutions in the initial stage. However, the decaying amplitudes with F cpf − and F pf − on 80 ×
240 meshes reach the steady state faster than the other forms as timeincreases. We found that this behavior can be improved by increasing the Peclet number.We further repeated the above simulations with ρ /ρ = 10. The results are shown inFig. 9. In this case, all the results give a good agreement with the theoretical solutions.For quantitative comparison, the time averaged L -norm error for the wave amplitude is18easured, which is defined as E ( h ) = (cid:115) ω (cid:90) ω | ¯ h ( t ) − ¯ h exact ( t ) | dt. (50)Table III presents the time averaged L -norm error of wave amplitude, from which we canobserve that all the averaged errors monotonically decrease as the numerical grid increases.Among the results, it can be found that the results given by F pf − and F csf − are closer tothe analytical solutions. t -0.8-0.6-0.4-0.200.20.40.60.81 h ( t ) / h TheoreticalH=80H=160 (a) t -0.8-0.6-0.4-0.200.20.40.60.81 h ( t ) / h TheoreticalH=80H=160 (b) t -0.8-0.6-0.4-0.200.20.40.60.81 h ( t ) / h TheoreticalH=80H=160 (c) t -0.8-0.6-0.4-0.200.20.40.60.81 h ( t ) / h TheoreticalH=80H=160 (d) t -0.8-0.6-0.4-0.200.20.40.60.81 h ( t ) / h TheoreticalH=80H=160 (e) t -0.8-0.6-0.4-0.200.20.40.60.81 h ( t ) / h TheoreticalH=80H=160 (f) t -0.8-0.6-0.4-0.200.20.40.60.81 h ( t ) / h TheoreticalH=80H=160 (g) t -0.8-0.6-0.4-0.200.20.40.60.81 h ( t ) / h TheoreticalH=80H=160 (h)
Fig. 8: Time evolution of capillary wave amplitude with ρ /ρ = 1 for (a) F stf − ,(b) F stf − ,(c) F cpf − ,(d) F cpf − ,(e) F pf − ,(f) F pf − ,(g) F csf − and (h) F csf − .TABLE III: The time averaged L -norm error for capillary wave time evolution. ρ ρ grid F stf − F stf − F cpf − F cpf − F pf − F pf − F csf − F csf − t -0.6-0.4-0.200.20.40.60.811.2 h ( t ) / h TheoreticalH=80H=160 (a) t -0.6-0.4-0.200.20.40.60.811.2 h ( t ) / h TheoreticalH=80H=160 (b) t -0.6-0.4-0.200.20.40.60.811.2 h ( t ) / h TheoreticalH=80H=160 (c) t -0.6-0.4-0.200.20.40.60.811.2 h ( t ) / h TheoreticalH=80H=160 (d) t -0.6-0.4-0.200.20.40.60.811.2 h ( t ) / h TheoreticalH=80H=160 (e) t -0.6-0.4-0.200.20.40.60.811.2 h ( t ) / h TheoreticalH=80H=160 (f) t -0.6-0.4-0.200.20.40.60.811.2 h ( t ) / h TheoreticalH=80H=160 (g) t -0.6-0.4-0.200.20.40.60.811.2 h ( t ) / h TheoreticalH=80H=160 (h)
Fig. 9: Time evolution of capillary wave amplitude with ρ /ρ = 10 for (a) F stf − ,(b) F stf − ,(c) F cpf − ,(d) F cpf − ,(e) F pf − ,(f) F pf − ,(g) F csf − and (h) F csf − . D. Rising bubble
We now examine the performance of different interfacial force formulations by simulatinga bubble rising in a two-dimensional domain, which was also simulated by Hysing et al. [52]and S. Aland et al . [53]. Although no analytical solution is available for this problem, somenumerical results were presented in [52]. The results from group 3 on the finest grids in[52] are taken as the reference solutions. The schematic of the domain is shown in Fig 10.Initially, a bubble with radius R = 0 .
25m is placed at (0 . , . × ×
240 and 240 ×
480 are employed. The fluid parameters are listed in Table IV.The related non-dimensional numbers are given byRe = ρ U g Lµ , Eo = ρ U g Lσ ,
Mo = Eo Re = U g µ ρ σ R (51)where U g = √ Rg and L are the reference velocity and length, respectively.For comparison, the benchmark quantities, including bubble shape at t = 3 s , rising20ig. 10: Initial configuration for the rising bubble.velocity, center of mass and circularity are measured by y c = (cid:82) Ω (1 − φ ) yd x (cid:82) Ω (1 − φ ) d x , (52) v c = (cid:82) Ω (1 − φ ) vd x (cid:82) Ω (1 − φ ) d x , (53) C = perimenter of area-equivalent circleperimeter of bubble = 2 (cid:113)(cid:82) φ< πdxP b (54)where v is the velocity component in the vertical direction and P b is obtained by integrationover the contour line at φ = 0 in Matlab.Figure 11 shows the bubble shapes predicted by various interfacial force formulationsat t = 3s. It can be seen that all the results agree well with the benchmark solutions.However, the shapes of the bubble obtained by F cpf − , F cpf − and F pf − are clearly lowerthan the reference solutions for both grids. Figure 12 shows the time histories of the centerof mass. At the initial stage, all the results are in good agreement with those reported in [53].However, the discrepancy between the results with F cpf − , F cpf − and F pf − and the referencesolutions becomes larger after t = 1 . ρ (kg / m ) ρ (kg / m ) µ (Pa · s) µ (Pa · s) g (m/ s ) σ (N/m) Eo Re1000 100 10 1 0.98 24.5 10 35 well with the reference values. However, the results with F pf − and F pf − on the coarse meshdeviate slightly from the reference solutions. The minimum circularity on the finer mesh issignificantly lower than that of the reference solution except for F pf − .For quantitative comparison, the maximum mass center position, the maximum risingvelocity and minimum circularity with each force formulation are calculated and comparedwith the reference results. The results are presented in Table V. Overall, the values obtainedby F stf − , F stf − , F csf − and F csf − are similar and in better agreement with the referencedata. x y S. Hysing (2009)S. Aland and A.Voigt (2012)120 240240 480 (a) x y S. Hysing (2009)S. Aland and A.Voigt (2012)120 240240 480 (b) x y S. Hysing (2009)S. Aland and A.Voigt (2012)120 240240 480 (c) x y S. Hysing (2009)S. Aland and A.Voigt (2012)120 240240 480 (d) x y S. Hysing (2009)S. Aland and A.Voigt (2012)120 240240 480 (e) x y S. Hysing (2009)S. Aland and A.Voigt (2012)120 240240 480 (f) x y S. Hysing (2009)S. Aland and A.Voigt (2012)120 240240 480 (g) x y S. Hysing (2009)S. Aland and A.Voigt (2012)120 240240 480 (h)
Fig. 11: Bubble shapes at time t = 3 s for (a) F stf − ,(b) F stf − ,(c) F cpf − ,(d) F cpf − ,(e) F pf − ,(f) F pf − ,(g) F csf − and (h) F csf − . E. Droplet deformation in shear flow
Finally, we consider a circle drop deformation in a shear flow. The schematic of theflow field is shown in Fig 15. Initially, a circle drop is located at the center of a rectangle22
T (s) Y c ( m ) S. Aland and A.Voigt (2012)120 240240 480 (a)
T (s) Y c ( m ) S. Aland and A.Voigt (2012)120 240240 480 (b)
T (s) Y c ( m ) S. Aland and A.Voigt (2012)120 240240 480 (c)
T (s) Y c ( m ) S. Aland and A.Voigt (2012)120 240240 480 (d)
T (s) Y c ( m ) S. Aland and A.Voigt (2012)120 240240 480 (e)
T (s) Y c ( m ) S. Aland and A.Voigt (2012)120 240240 480 (f)
T (s) Y c ( m ) S. Aland and A.Voigt (2012)120 240240 480 (g)
T (s) Y c ( m ) S. Aland and A.Voigt (2012)120 240240 480 (h)
Fig. 12: The evolution of the center of mass for (a) F stf − ,(b) F stf − ,(c) F cpf − ,(d) F cpf − ,(e) F pf − ,(f) F pf − ,(g) F csf − and (h) F csf − . T (s) v c ( m / s ) S. Aland and A.Voigt (2012)120 240240 480 (a)
T (s) v c ( m / s ) S. Aland and A.Voigt (2012)120 240240 480 (b)
T (s) v c ( m / s ) S. Aland and A.Voigt (2012)120 240240 480 (c)
T (s) v c ( m / s ) S. Aland and A.Voigt (2012)120 240240 480 (d)
T (s) v c ( m / s ) S. Aland and A.Voigt (2012)120 240240 480 (e)
T (s) v c ( m / s ) S. Aland and A.Voigt (2012)120 240240 480 (f)
T (s) v c ( m / s ) S. Aland and A.Voigt (2012)120 240240 480 (g)
T (s) v c ( m / s ) S. Aland and A.Voigt (2012)120 240240 480 (h)
Fig. 13: The evolution of the rising velocity for (a) F stf − ,(b) F stf − ,(c) F cpf − ,(d) F cpf − ,(e) F pf − ,(f) F pf − ,(g) F csf − and (h) F csf − .domain of 2 H × H . The effect of gravity is ignored. The top and bottom walls maintainvelocities U and − U , respectively, leading to a shear rate E = 2 U/H . The periodic boundaryconditions are applied to the left and right boundaries. The same density and viscosity23
T(s) C S. Aland and A. Voigt (2012)120 120240 480 (a)
T(s) C S. Aland and A. Voigt (2012)120 120240 480 (b)
T(s) C S. Aland and A. Voigt (2012)120 120240 480 (c)
T(s) C S. Aland and A. Voigt (2012)120 120240 480 (d)
T(s) C S. Aland and A. Voigt (2012)120 120240 480 (e)
T(s) C S. Aland and A. Voigt (2012)120 120240 480 (f)
T(s) C S. Aland and A. Voigt (2012)120 120240 480 (g)
T(s) C S. Aland and A. Voigt (2012)120 120240 480 (h)
Fig. 14: The evolution of circularity for (a) F stf − ,(b) F stf − ,(c) F cpf − ,(d) F cpf − ,(e) F pf − ,(f) F pf − ,(g) F csf − and (h) F csf − .TABLE V: Benchmark quantities for rising bubble on 240 ×
480 meshes.
Parameter Ref.[52] F stf − F stf − F cpf − F cpf − F pf − F pf − F csf − F csf − y max v max C min are specified for both the drop and surrounding fluid. In the simulation, we set H =8m , R = 1m , U w = 4m/s , ρ d = ρ s = 1kg/m . The Reynolds number Re = Eρ d R /µ d = 0 . µER/σ is varied from 0 . . σ . The uniformgrid size of 200 ×
200 is employed. The shapes of the deformed drop at steady state areillustrated in Fig. 16. It can be seen that the shapes of the drop given by all interfacial forceformulations deform into an ellipsoidal one and are elongated as Ca increases. In particular,the shapes of the drop obtained by F csf − are overstretched compared to other results.The shape of the drop can be characterized by a Taylor deformation parameter defined as D = ( L − B ) / ( L + B ), where L and B are the lengths along the major axis and the minoraxis of the droplet, respectively. A theoretical solution derived on the assumptions of theStokes flow and small deformation shows that the Taylor deformation parameter is related24ig. 15: Drop deformation in a shear flow. L is the major axis and B is the minor axis.to the capillary number and the viscosity ratio [56, 57] D = L − BL + B = Ca 19 λ + 1616 λ + 16 , (55)where λ = µ d /µ f is the viscosity ratio between the drop fluid and the surrounding fluid.Table VI shows the Taylor deformation parameters with different force formulations. Itcan be seen that the values predicted by F csf − are significantly higher than the theoreticalvalues. Overall, the values with F stf − and F stf − are close to the theoretical ones.TABLE VI: Comparison of Taylor deformation number D with linear theory F sf ρ d /ρ s = 1 ρ d /ρ s = 0 . F stf − F stf − F cpf − F cpf − F pf − F pf − F csf − F csf − x y Ca=0.1Ca=0.2Ca=0.3Ca=0.4 (a) x y Ca=0.1Ca=0.2Ca=0.3Ca=0.4 (b) x y Ca=0.1Ca=0.2Ca=0.3Ca=0.4 (c) x y Ca=0.1Ca=0.2Ca=0.3Ca=0.4 (d) x y Ca=0.1Ca=0.2Ca=0.3Ca=0.4 (e) x y Ca=0.1Ca=0.2Ca=0.3Ca=0.4 (f) x y Ca=0.1Ca=0.2Ca=0.3Ca=0.4 (g) x y Ca=0.1Ca=0.2Ca=0.3Ca=0.4 (h)
Fig. 16: The φ = 0 contours of the order parameter at Ca = 0 . , . , . ×
200 meshwith (a) F stf − ,(b) F stf − ,(c) F cpf − ,(d) F cpf − ,(e) F pf − ,(f) F pf − ,(g) F csf − and (h) F csf − . VI. CONCLUSIONS
In this paper, we successfully implemented the phase-field-based lattice Bolzmann methodwith different interfacial force formulations for two phases flow. The performance of eachsurface tension formulation has been validated and compared. For a stationary drop, F csf − provides the most accurate prediction in terms of the surface tension coefficient. The poten-tial form can generate small spurious currents. F stf − , F csf − and F csf − produce a smoothpressure field across the interface and symmetric distribution of the interfacial force. Thedistributions of F stf − , F cpf − F pf − and F pf − become symmetric with respect to the phaseinterface. For the droplets merging problems, there are obvious differences for the interfaceshapes of the droplets during coalescence. The droplets is more prone to merge due to thesurface tension effects when F cpf − , F pf − and F csf − are used. In particular, the unexpectedmovement of droplets with unequal sizes occurs when F cpf − , F cpf − , F csf − and F csf − areused. It is also found that F stf and F csf show better numerical stability than F cpf and F csf .For the test of capillary wave, the evolution processes of the interface amplitude from F pf − and F csf − are closer to the analytical solutions in all formulations. It is worth noting that F cpf − can yield good results but the Peclet number should be carefully chosen. For the26imulation of a rising bubble, both the stress form and CSF form give a good results in termsof the mass center. F cpf − , F cpf − and F pf − clearly underestimate the center at the latestage. For the rising velocity, all formulations underestimate the maximum rise velocity. Interms of the circularity, only F cpf − and F csf − give the predictions closer to the referencesolutions. For shear flow, all formulations give accurate predictions in comparison with thelinear theory at Ca = 0 .
1. With the increase of capillary number, F csf − produces a largerdeformation than the theoretical predictions. For all the considered capillary number, F stf − and F stf − can give a satisfactory prediction.In summary, it seems that no surface tension force formulation can give satisfactory resultsin all tests. Different forms may be considered for different problems. Overall, F cpf is goodfor calculating multiphase flows with small interface deformation. Both F stf and F csf aregood for dynamical situations. We hope the present comparison can provide insights intothe advantages and limitations of each formulation. DATA AVAILABILITY
The data that support the findings of this study are available from the correspondingauthor upon reasonable request.
ACKNOWLEDGEMENTS
This study was supported by the National Science Foundation of China(51836003).
Appendix A: Relations among different interfacial force formulations
This appendix presents the relations among different interfacial force formulations. InEq.(24), the curvature term can be written as ∇ · n = ∇ · (cid:18) ∇ φ |∇ φ | (cid:19) = 1 |∇ φ | (cid:18) ∇ φ − ∇ φ · ∇|∇ φ ||∇ φ | (cid:19) . (A1)Substituting the above equation into Eq.(24) yields F csf − = − κ ∇ φ |∇ φ |∇ · n = − κ ∇ φ (cid:18) ∇ φ − ∇ φ · ∇|∇ φ ||∇ φ | (cid:19) , = F pf − + κ ∇ φ ( ∇ φ · ∇|∇ φ | ) |∇ φ | , (A2)27y using the equality ∇ φ ∇ φ = ∇ ( φ ∇ φ ) − φ ∇∇ φ , one can obtain the following relation-ship, F csf − = F pf − − κ ∇ ( φ ∇ φ ) + κ ∇ φ ( ∇ φ · ∇|∇ φ | ) |∇ φ | , (A3)Based on Eq.(2), F csf − can be rewritten as F csf − = ∇ φ (cid:18) µ φ − ∂f ∂φ (cid:19) + κ ∇ φ ( ∇ φ · ∇|∇ φ | ) |∇ φ | , = F cpf − − ∇ f + κ ∇ φ ( ∇ φ · ∇|∇ φ | ) |∇ φ | , = F cpf − + ∇ ( φµ φ ) − ∇ f + κ ∇ φ ( ∇ φ · ∇|∇ φ | ) |∇ φ | , (A4)where we have used the equality ∇ ( φµ φ ) = µ φ ∇ φ + φ ∇ µ φ . By using the following equality, − κ ∇ φ ∆ φ = κ ∇|∇ φ | − ∇ · κ ( ∇ φ ⊗ ∇ φ ) , (A5)one can obtain F csf − = − κ ∇ φ ∇ φ + κ ∇ φ ( ∇ φ · ∇|∇ φ | ) |∇ φ | , = F stf − + κ ∇|∇ φ | + κ ∇ φ ( ∇ φ · ∇|∇ φ | ) |∇ φ | , (A6)By virtue of Eq.(4), we have |∇ φ | = ∂φ∂r = 2 W (1 − φ ) , ∇ φ ( ∇ φ · ∇|∇ φ | ) |∇ φ | = 12 ∇|∇ φ | = 8 W φ ( φ − ∇ φ, (A7)Eq.(A6) is then derived as F csf − = −∇ · κ ( ∇ φ ⊗ ∇ φ ) + κ ∇|∇ φ | = F stf − . (A8) Appendix B: The Dirac function in F csf − and F csf − In F csf − , the surface Dirac function is chosen as α |∇ φ | with α being an undeterminedparameter. Based on Eq.(4), |∇ φ | = ∂φ∂r = 2 W (1 − φ ) , (B1)28nserting the above equation into Eq.(23) yields, (cid:90) ∞−∞ α |∇ φ | dr = (cid:90) ∞−∞ α W (1 − φ ) ∂φ∂r dr = (cid:90) ∞−∞ α W (1 − φ ) dφ = α W (cid:90) ∞−∞ d (cid:18) φ − φ (cid:19) = 8 α W = 1 , (B2)where φ | r = ∞ = 1 and φ | r = −∞ = − α = W .In F csf − , the surface Dirac function is chosen as α |∇ φ | . Analogously, one can have (cid:90) ∞−∞ α |∇ φ | dr = (cid:90) ∞−∞ ∂φ∂r dr = (cid:90) ∞−∞ αdφ = 2 α = 1 . (B3)This leads to α = . REFERENCES [1] V. Cristini and Y.-C. Tan, Lab on a Chip , 257 (2004).[2] L. Clime, D. Brassard, and T. Veres, Journal of Applied Physics , 07B517 (2009).[3] H. Liu, Y. Ba, L. Wu, Z. Li, G. Xi, and Y. Zhang, Journal of Fluid Mechanics , 381(2018).[4] S. Shad, I. Gates, et al. , Journal of Canadian Petroleum Technology , 48 (2010).[5] A. Jafari, M. Hasani, M. Hosseini, and R. Gharibshahi, Petroleum Science , 1 (2019).[6] S. McKee, M. F. Tom´e, V. G. Ferreira, J. A. Cuminato, A. Castelo, F. Sousa, and N. Man-giavacchi, Computers & Fluids , 907 (2008).[7] G. Tryggvason, R. Scardovelli, and S. Zaleski, Direct numerical simulations of gas–liquidmultiphase flows (Cambridge University Press, 2011).[8] C. W. Hirt and B. D. Nichols, Journal of computational physics , 201 (1981).[9] Y.-C. Chang, T. Hou, B. Merriman, and S. Osher, Journal of computational Physics ,449 (1996).[10] D. Jacqmin, Journal of Computational Physics , 96 (1999).[11] D. J. Korteweg, Archives N´eerlandaises des Sciences exactes et naturelles , 1 (1901).
12] D. M. Anderson, G. B. McFadden, and A. A. Wheeler, Annual review of fluid mechanics ,139 (1998).[13] H. Liang, B. Shi, Z. Guo, and Z. Chai, Physical Review E , 053320 (2014).[14] H. Liu, A. J. Valocchi, Y. Zhang, and Q. Kang, Journal of Computational Physics , 334(2014).[15] A. Fakhari, Y. Li, D. Bolster, and K. T. Christensen, Advances in water resources , 119(2018).[16] J. Lowengrub and L. Truskinovsky, Proceedings of the Royal Society of London. Series A:Mathematical, Physical and Engineering Sciences , 2617 (1998).[17] H. Abels, H. Garcke, and G. Gr¨un, Mathematical Models and Methods in Applied Sciences , 1150013 (2012).[18] K. F. Lam and H. Wu, European Journal of Applied Mathematics , 595 (2018).[19] T. Lee and P. F. Fischer, Physical Review E , 046709 (2006).[20] X. He, S. Chen, and R. Zhang, Journal of Computational Physics , 642 (1999).[21] A. Fakhari and M. H. Rahimian, International journal for numerical methods in fluids , 827(2010).[22] A. Shah, S. Saeed, and S. A. Khan, Heliyon , e01024 (2018).[23] S. A. Khan and A. Shah, AIP Advances , 085312 (2019).[24] Y. Zu and S. He, Physical Review E , 043301 (2013).[25] D. Jacqmin, in (1996) p. 858.[26] H. Ding, P. D. Spelt, and C. Shu, Journal of Computational Physics , 2078 (2007).[27] A. Fakhari and M. H. Rahimian, Physical Review E , 036707 (2010).[28] J. U. Brackbill, D. B. Kothe, and C. Zemach, Journal of computational physics , 335(1992).[29] J. Kim, Journal of Computational Physics , 784 (2005).[30] H. G. Lee and J. Kim, International Journal for Numerical Methods in Engineering , 269(2012).[31] J. Chao, R. Mei, R. Singh, and W. Shyy, International journal for numerical methods in fluids , 622 (2011).[32] J. W. Cahn and J. E. Hilliard, The Journal of chemical physics , 258 (1958).[33] P. Yue, J. J. Feng, C. Liu, and J. Shen, Journal of Fluid Mechanics , 293 (2004).
34] H. Abels, H. Garcke, G. Gr¨un, and S. Metzger, in
Transport Processes at Fluidic Interfaces (Springer, 2017) pp. 203–229.[35] V. Starovoitov, Journal of applied mechanics and technical physics , 891 (1994).[36] D. Jacqmin, Journal of Fluid Mechanics , 57 (2000).[37] X. Yang, J. J. Feng, C. Liu, and J. Shen, Journal of Computational Physics , 417 (2006).[38] J. Shen and X. Yang, Journal of computational physics , 2978 (2009).[39] H.-G. Lee, J. Lowengrub, and J. Goodman, Physics of Fluids , 492 (2002).[40] J. Kim, Applied mathematics and computation , 589 (2005).[41] C. Zhang, Z. Guo, and Y. Li, International Journal of Heat and Mass Transfer , 1128(2019).[42] Y. Wang, C. Shu, H. Huang, and C. Teo, Journal of Computational Physics , 404 (2015).[43] Z. Chen, C. Shu, D. Tan, X. Niu, and Q. Li, Physical Review E , 063314 (2018).[44] S. Popinet, Annual Review of Fluid Mechanics , 49 (2018).[45] C. Zhang, Z. Guo, and H. Liang, Physical Review E , 043310 (2019).[46] F. Magaletti, F. Picano, M. Chinappi, L. Marino, and C. M. Casciola, Journal of FluidMechanics , 95 (2013).[47] Z. Guo, C. Zheng, and B. Shi, Physical Review E , 036707 (2011).[48] T. Kr¨uger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, and E. M. Viggen, SpringerInternational Publishing , 4 (2017).[49] A. J. Ladd, Journal of fluid mechanics , 285 (1994).[50] A. Prosperetti, The Physics of Fluids , 1217 (1981).[51] P. Lallemand, L.-S. Luo, and Y. Peng, Journal of Computational Physics , 1367 (2007).[52] S.-R. Hysing, S. Turek, D. Kuzmin, N. Parolini, E. Burman, S. Ganesan, and L. Tobiska,International Journal for Numerical Methods in Fluids , 1259 (2009).[53] S. Aland and A. Voigt, International Journal for Numerical Methods in Fluids , 747 (2012).[54] S. H. Kim and H. Pitsch, Journal of Computational Physics , 19 (2015).[55] L. Amaya-Bower and T. Lee, Computers & fluids , 1191 (2010).[56] G. I. Taylor, Proceedings of the Royal Society of London. Series A, containing papers of amathematical and physical character , 501 (1934).[57] G. I. Taylor, Proceedings of the Royal Society of London. Series A, Containing Papers of aMathematical and Physical Character , 41 (1932)., 41 (1932).