A versatile lattice Boltzmann model for immiscible ternary fluid flows
AA versatile lattice Boltzmann model for immiscible ternaryfluid flows
Yuan Yu a , Haihu Liu b, ∗ , Dong Liang a,c , Yonghao Zhang d a School of Engineering, Sun Yat-Sen University, Guangzhou 510006, China b School of Energy and Power Engineering, Xi’an Jiaotong University, Xi’an 710049, China c Guangdong Provincial Key Laboratory of Fire Science and Technology, Guangzhou 51006, China d Department of Mechanical and Aerospace Engineering, University of Strathclyde, Glasgow G1 1XJ,United Kingdom
Abstract
We propose a lattice Boltzmann color-gradient model for immiscible ternary fluidflows, which is applicable to the fluids with a full range of interfacial tensions, es-pecially in near-critical and critical states. An interfacial force for N-phase systems isderived based on the previously developed perturbation operator and is then introducedinto the model using a body force scheme, which helps reduce spurious velocities. Ageneralized recoloring algorithm is applied to produce phase segregation and ensureimmiscibility of three di ff erent fluids, where a novel form of segregation parameters isproposed by considering the existence of Neumann’s triangle and the e ff ect of equilib-rium contact angle in three-phase junction. The proposed model is first validated withthree typical examples, namely the interface capturing for two separate static droplets,the Young-Laplace test for a compound droplet, and the spreading of a droplet betweentwo stratified fluids. This model is then used to study the structure and stability of dou-ble droplets in a static matrix. Consistent with the theoretical stability diagram, sevenpossible equilibrium morphologies are successfully reproduced by adjusting two ratiosof the interfacial tensions. By simulating Janus droplets in various geometric config-urations, the model is shown to be accurate when three interfacial tensions satisfy aNeumann’s triangle. In addition, we also simulate the near-critical and critical statesof double droplets where the outcomes are very sensitive to the model accuracy. Ourresults show that the present model is advantageous to three-phase flow simulations,and allows for accurate simulation of near-critical and critical states. Keywords:
Immiscible multiphase flows, full range of interfacial tensions,color-gradient model, double emulsion droplets, near-critical and critical states ∗ Corresponding author: Tel.: +
86 (0) 298 266 5700;
Email address: [email protected] (Haihu Liu)
Preprint submitted to Elsevier August 28, 2018 a r X i v : . [ phy s i c s . c o m p - ph ] A ug . Introduction An emulsion is a mixture of a dispersed phase as droplets in another immisciblefluid that forms a continuous phase. Two basic types of emulsions are the oil-in-water(O / W) and water-in-oil (W / O) emulsions[1]. Recently, more complex systems referredto as double emulsions and Janus emulsions have received a rapidly growing interestdue to their unique properties[2, 3, 4] and potential applications[5, 6, 7, 8, 9, 10]. Dou-ble emulsions, also known as ‘emulsion of emulsion’ or ‘emulsion within emulsion’,are emulsions with smaller droplets encapsulated in larger droplets. The shell fluid canserve as a barrier between the core droplets and the outer environment, which makesdouble emulsions highly desirable for applications in controlled release, separation,and encapsulation[1, 2, 3, 4]. Janus emulsions, which are named after the two-facedRoman god Janus, are highly structured fluids consisting of emulsion droplets thathave two distinct physical properties[11]. Because of their natural asymmetric abil-ity in the compositions and the shapes, Janus emulsions are often used in the fieldsthat need asymmetry in the shape and the materials. In the applications of emulsions,morphology is one of the most important properties and closely related to other emul-sion properties such as rheology, droplet size, relative stability, electrical conductivityand zeta potential[2, 3, 4, 12]. A number of theoretical and experimental studies havebeen devoted to identifying di ff erent equilibrium morphologies and their transforma-tion. For example, Torza and Mason[13] studied the droplet morphology in termsof spreading coe ffi cients and obtained the theoretical relationship between the dropletmorphology and spreading coe ffi cients. They experimentally observed three equilib-rium morphologies of double droplets, i.e. complete engulfing, partial engulfing andnon-engulfing, which correspond to three sets of spreading coe ffi cients. Beyond thesethree equilibrium states, Pannacci et al.[14] identified several new morphologies ofdouble droplets, and found the non-equilibrium morphologies can have long lifetimescontrolled by hydrodynamics, which facilitates the use of double droplets to produceencapsulated particles at early times and Janus particles at longer times. Guzowski etal.[15] presented a detailed theoretical analysis on the possible equilibrium morpholo-gies of double droplets and designed the structure of double emulsions by tuning thevolumes of the constituent segments experimentally. As a supplement to theoreticaland experimental studies, numerical modelling and simulations are becoming increas-ingly popular in investigation of the behavior of Janus / double emulsions, which aretypical of three-phase flow problems.Traditionally, three-phase flows are simulated by solving the macroscopic Navier-Stokes equations together with various approaches to capturing or tracking the inter-faces between fluids. Among these approaches, the front-tracking[16], volume-of-fluid(VOF)[17], level-set[18, 19, 20] and phase-field[21, 22, 23, 24, 25, 26] methods arecommonly used. However, the front-tracking method is not suitable for simulatinginterface breakup and coalescence; the VOF and level-set methods require either so-phisticated interface reconstruction algorithms or unphysical re-initialization processesto represent the interfaces; and the phase-field method yields an interface thicknessfar greater than its actual value, which may lead to unphysical dissolution of smalldroplets and mobility-dependent numerical results[27]. It still remains an open ques-tion for the phase-field method to choose an optimal mobility, even for a two-phase2ow problem[28].In the past decades, the lattice Boltzmann (LB) method has developed into a promis-ing alternative to the traditional Navier-Stokes-based solvers, for simulating complexflow problems. It is a pseudo-molecular method tracking evolution of the distribu-tion function of an assembly of molecules, built upon microscopic models and meso-scopic kinetic equations[29]. The LB method has several advantages over the tradi-tional Navier-Stokes-based solvers, e.g. the algorithm simplicity and parallelizability,and the ease of handling complex boundaries[30]. In addition, its kinetic nature allowsa simple incorporation of microscopic physics without su ff ering from the limitations interms of length and time scales typical of molecular dynamics simulations[31]. Thus,the LB method is particularly useful in the simulation of multiphase flows. The existingLB models for multiphase flows can be generally classified into four categories: color-gradient model[32, 33, 28], interparticle-potential model[34, 35, 36, 37, 38], phase-field-based model[39, 40, 41], and mean-field theory model[42]. These models haveshown great success as in dealing with two-phase flow problems, and all of them exceptthe mean-field theory model have been extended to the modeling of immiscible ternaryfluids, see, e.g. Refs[43, 44, 45, 46, 47, 48, 49, 50, 51]. The ternary color-gradientmodels[49, 50] inherit a series of advantages of its two-phase counterpart, such asstrict mass conservation for each fluid, flexibly tunable interfacial tensions, and thestability for a broad range of viscosity ratios, and they are well suited to exploring thedynamic processes occurring in ternary fluid systems as previously demonstrated by Fuet al.[51] and Jiang et al.[52]. The existing color-gradient models, however, commonlysu ff er from a problem, i.e. three interfacial tensions should satisfy a Neumann’s trian-gle. In industrial processes, surfactants are often added to emulsions to stabilize themagainst droplet coalescence. The presence of surfactants could significantly modify theinterfacial tensions so that the interfacial tensions do not always yield a Neumann’striangle. To correctly predict the dynamical behavior of emulsions, thereby allowingprecise control over the droplet geometry and composition, it is necessary for a nu-merical model to be capable of simulating ternary fluids with a full range of interfacialtensions. On the other hand, it is challenging to simulate the near-critical and criticalstates of a ternary fluid system where the largest interfacial tension is close to the sumof the other two, as the outcomes are very sensitive to the model accuracy.In this paper, we develop a LB color-gradient model for simulating immiscibleternary fluids with a full range of interfacial tensions. Based on the perturbation op-erator developed by Leclaire et al.[50], an interfacial force formulation is derived todescribe the interactions among di ff erent fluids and is then introduced into the modelusing a body force scheme, which is found to e ff ectively reduce spurious velocities. Inaddition, the recoloring algorithm proposed by Spencer et al.[49] is applied to maintainthe interfaces and ensure immiscibility of three di ff erent fluids, where a new form ofsegregation parameters is proposed by considering both the existence of Neumann’striangle and the e ff ect of equilibrium contact angle in three-phase junction. The capa-bility and accuracy of this model are first assessed by simulating the interface capturingfor two separate static droplets, the Young-Laplace test for a compound droplet, and thespreading of a droplet between two stratified fluids. It is then used to study the struc-ture and stability of double droplets in a static matrix fluid, where we emphasize themodel’s capability for simulating ternary fluid flows in near-critical and critical states.3 . Numerical method The two-phase color-gradient LB model of Liu et al.[28, 53] is extended to thesimulation of immiscible ternary fluids. The ternary color-gradient model consists ofthree steps, i.e. the collision step, the recoloring step and the streaming step. In thecollision step, an interfacial force that describes the interactions among di ff erent fluidsis derived from the perturbation operator presented in Leclaire et al.[50], and is thenintroduced by the body force scheme of Guo et al.[54] In the recoloring step, a novelform of segregation parameters is proposed to ensure accurate phase segregation inthree-phase junction and allow for the states where three interfacial tensions betweenthe fluids cannot form a triangle, known as the Neumann’s triangle. The distributionfunctions f i , r , f i , g and f i , b are introduced to represent three immiscible fluids, i.e. redfluid, green fluid and blue fluid, where the subscript i is the lattice velocity directionand ranges from 0 to ( n -1) for a given m -dimensional D m Q n lattice model. The totaldistribution function is defined as f i = (cid:80) k f i , k ( k = r , g or b ), which undergoes acollision step as f † i ( x , t ) = f i ( x , t ) + Ω i ( x , t ) + Φ i ( x , t ) , (1)where f i ( x , t ) is the total distribution function in the i -th velocity direction at the posi-tion x and the time t , f † i is the post-collision distribution function, Ω i is the Bhatnagar-Gross-Krook (BGK) collision operator, and Φ i is the forcing term (also known as per-turbation operator), which contributes to the mixed interfacial regions and creates theinterfacial tensions between di ff erent fluids.In the BGK collision operator, the total distribution functions are relaxed toward alocal equilibrium with a single relaxation time: Ω i ( x , t ) = − τ f (cid:104) f i ( x , t ) − f eqi ( x , t ) (cid:105) , (2)where τ f is the dimensionless relaxation time, and f eqi is the equilibrium distributionfunction of f i . The equilibrium distribution function is obtained by a second orderTaylor expansion of Maxwell-Boltzmann distribution with respect to the local fluidvelocity u : f eqi = w i ρ (cid:34) + e i · u c s + ( e i · u ) c s − u · u c s (cid:35) , (3)where ρ = (cid:80) k ρ k is the total density and ρ k is the density of the fluid k ; c s is the speed ofsound; e i is the lattice velocity in the i -th direction; and w i is the weighting factor. Forthe two-dimensional nine-velocity (D2Q9) model, e i is defined as e = (0 , e , = ( ± c , e , = (0 , ± c ), e , = ( ± c , ± c ), and e , = ( ∓ c , ± c ), where c = δ x /δ t = √ c s with δ x and δ t being the lattice spacing and time step, respectively (for the sake ofsimplicity, δ x = δ t = w i is given by w = / w − = / w − = / ff erent fluids in three-phase simulations. FollowingLeclaire et al.[50], the perturbation operator is given by Φ i = (cid:88) k Φ i , k , (4) Φ i , k = (cid:88) l , l (cid:44) k A kl C kl | G kl | (cid:34) w i ( e i · G kl ) | G kl | − B i (cid:35) , (5)where G kl = ρ l ρ ∇ ρ k ρ − ρ k ρ ∇ ρ l ρ is the color gradient [50] and is introduced to identify thelocation of the k - l interface, i.e. the interface between the fluid k and the fluid l . C kl is a concentration factor that controls the activation of the interfacial tension at the k - l interface, and is given by [50] C kl = min ρ k ρ l ρ k ρ l , , (6)where ρ k is the density of the pure fluid k , and A kl is a parameter related to the interfacialtension between the fluids k and l , i.e. σ kl = ( A kl + A lk ) τ f . The generalized expres-sion for B i was given by Liu et al.[28] and it was in particular taken as B = − / B − = /
27 and B − = /
108 in the work of Leclaire et al.[50]. It is worth notingthat Eq.(5) is not limited to the case with ternary fluids, and can be also applicable to N -phase ( N >
3) systems.Using the Chapman-Enskog multiscale analysis, it is shown that the perturbationoperator, given by Eqs.(4) and (5), can lead to the following interfacial force: F s = −∇ · τ f δ t (cid:88) i Φ i e i e i = (cid:88) k (cid:88) l , l (cid:44) k ∇ · (cid:20) σ kl C kl | G kl | ( I − n kl n kl ) (cid:21) , (7)where n kl is the unit normal vector of the k - l interface and is defined by n kl = G kl / | G kl | .Instead of using Eqs.(4) and (5), the e ff ect of interfacial tension is realized throughthe body force scheme of Guo et al.[54], which is able to reduce e ff ectively spuriousvelocities while keeping high numerical accuracy [53, 55]. According to Guo et al.[54],the forcing term Φ i in Eq. (1) is written as Φ i ( x , t ) = w i (cid:32) − τ f (cid:33) (cid:32) e i − u c s + e i · u c s e i (cid:33) · F s ( x , t ) δ t , (8)where the local fluid velocity is defined by the averaged momentum before and afterthe collision, i.e., ρ u ( x , t ) = (cid:88) i f i ( x , t ) e i + F s ( x , t ) δ t . (9)In this work, we assume equal densities for the red, green and blue fluids. To allowfor unequal viscosities of the three fluids, we determine the local kinematic viscosity ν by a harmonic mean ρν = (cid:88) k ρ k ν k , (10)5here ν k ( k = R , G or B ) is the kinematic viscosity of the fluid k . The local relaxationtime τ f can be calculated from the local viscosity using the following equation: ν = (cid:32) τ f − (cid:33) c s δ t . (11)The partial derivatives in the interfacial force F s should be evaluated through suit-able di ff erence schemes. To minimize the discretization errors, the fourth-order isotropicfinite di ff erence scheme ∂ α ϕ ( x , t ) = c s (cid:88) i w i ϕ ( x + e i δ t , t ) e i α , (12)is used to evaluate the derivatives of a variable ϕ .Although the forcing term generates the interfacial tensions, it does not guaranteethe immiscibility of di ff erent fluids. In order to minimize the mixing of the fluids, arecoloring step is applied. Based on the pioneering work of D’Ortona et al.[56], Latva-Kokko and Rothman[57] developed a recoloring algorithm to demix two immisciblefluids, which can overcome the lattice pinning problem and creates a symmetric dis-tribution of particles around the interface so that unphysical spurious velocities canbe e ff ectively reduced. This recoloring algorithm was later generalized by Spencer etal.[49] to three-phase fluid flows. Following Spencer et al.[49], the recolored distribu-tion functions of the fluid k ( k = r , g or b ) are f ‡ i , k ( x , t ) = ρ k ρ f † i ( x , t ) + (cid:88) l , l (cid:44) k β kl w i ρ k ρ l ρ n kl · e i , (13)where f ‡ i , k is the recolored distribution functions of the fluid k , and β kl is a segregationparameter related to the thickness of the k - l interface. It should be noted that β kl = β lk in order to conserve mass and momentum during the recoloring process. (cid:2026) (cid:3045)(cid:3034) (cid:2026) (cid:3045)(cid:3029) (cid:2026) (cid:3034)(cid:3029) (cid:2030) (cid:3045)(cid:3034) (cid:2030) (cid:3034)(cid:3029) (cid:2030) (cid:3045)(cid:3029) Figure 1: Neumann’s triangle
For the ternary fluids and when three interfacial tensions satisfy a Neumann’s tri-angle (see Fig. 1), the equilibrium contact angle ϕ kl will be formed between the fluidsin three-phase junction, and it is related to the interfacial tensions bycos( ϕ kl ) = σ mk + σ ml − σ kl σ mk σ ml . (14)6pencer et al. [49] theoretically showed that in three-phase junction, there should be arelationship between ϕ kl and the (relative) interface thickness, which is controlled bythe segregation parameter β kl . Hence, it is of great importance to select a proper β kl inthree-phase simulations. Several di ff erent forms of β kl have been provided in literature.Spencer et al. [49] proposed the first expression for the segregation parameters, whichis given by β rg = β β rb = β (cid:34) + ρ r ρ g ρ b ρ (cid:16) sin ϕ gb − (cid:17)(cid:35) β gb = β (cid:34) + ρ r ρ g ρ b ρ (sin ϕ rb − (cid:35) , (15)where β is the reference segregation parameter. Clearly, the segregation parameters inEq. (15) will degenerate into β kl = β at an interface where only two fluids are present.So it is suggested to take β = . β kl = β for the largest ϕ kl in the Neumann’striangle, β kl = β kl with ϕ max β + β C t (cid:2) sin ( π − ϕ max − ϕ kl ) − (cid:3) otherwise , (16)where ϕ max = max( ϕ kl ) and C t = min (cid:16) ρ r ρ g ρ b ρ , (cid:17) . Leclaire et al.[50] also mentionedto use β kl = β when the Neumann’s triangle does not exist. Clearly, Eq. (16) willdegradate to Eq. (15) when ϕ max = ϕ rg . Althogh Eqs.(15) and (16) work to someextent especially when the Neumann’s triangle exists, they cannot accurately simulatethe critical state where the largest interfacial tension equals the sum of the other two,which will be shown later. Recently, Fu et al.[51] seemed to have also noticed thatEqs.(15) and (16) do not always produce convincing results in three-phase simulations,so they simply selected a constant β kl , i.e. β kl = β . (17)It is evident that the dependence of β kl on ϕ kl is not considered in Eq.(17), and thusincorrect results may be obtained, e.g. in the critical state.To overcome the aforementioned drawbacks associated with the existing β kl , anovel form of segregation parameters is proposed. First, we determine whether theNeumann’s triangle exists by calculating X kl = σ mk + σ ml − σ kl σ mk σ ml . (18)It is easily seen from Eq.(14) that the Neumann’s triangle will exist if | X kl | < kl . Then, the segregation parameter β kl is defined as a continuous function of X kl : β kl = β + β min (cid:32) ρ r ρ g ρ b ρ , (cid:33) g ( X kl ) , (19)7here g ( X kl ) = X kl < − − sin (arccos ( X kl )) − ≤ X kl < X kl )) − ≤ X kl ≤ − < X kl . (20)It should be noted in three-phase junction that Eqs.(19) and (20) are derived basedon the following relationship: β rg sin( ϕ rg ) = β rb sin( ϕ rb ) = β gb sin( ϕ gb ) , (21)which is consistent with the nature of di ff use interfaces, thus leading to more accurateresults than using other forms of β kl . Moreover, the proposed β kl works well no matterif the Neumann’s triangle exists or not.After the recoloring step, the red, green and blue distribution functions propagateto the neighboring lattice nodes, known as the propagation or streaming step: f i , k ( x + e i δ t , t + δ t ) = f ‡ i , k ( x , t ) , k = { r , g , b } (22)with the post-propagation distribution functions used to compute the densities of col-ored fluids by ρ k = (cid:80) i f i , k .
3. Numerical Validations
We first consider two separate static droplets immersed in another fluid (say bluefluid) to validate the present model for interface capturing. Initially, a red droplet anda green droplet, both having equal radius R =
20, are placed in a 200 ×
100 latticedomain, and their centers are located at ( x r , y r ) = (50 ,
50) and ( x g , y g ) = (150 , y =
50 canbe analytically expressed as[58] ρ r ρ ( x ) = . + . R − (cid:112) ( x − x r ) ξ , (23a) ρ g ρ ( x ) = . + . R − (cid:112) ( x − x g ) ξ , (23b) ρ b ρ ( x ) = − ρ r ρ ( x ) − ρ g ρ ( x ) , (23c)8or the red, green and blue fluids, respectively. Here, the parameter ξ is a measure ofthe interface thickness related to β by ξ = / (6 k β ) [59], and k is a geometric constantthat is determined by [58] k = (cid:88) i w i e i e i | e i | . (24)For the D2Q9 model, one can obtain k ≈ . ξ ≈ . β = .
7. The simulation is run with the interfacial tensions σ rg = σ rb = σ gb = . ν r = ν g = ν b = .
1. Periodic boundary conditions are applied inboth the x and y directions. Fig. 2 shows the simulated density distributions of the red,green and blue fluids along y =
50 in the steady state, and the corresponding analyticalsolutions, given by Eq.(23), are also shown for comparison. Clearly, the simulateddensity distributions are all in good agreement with the analytical solutions, indicatingthat the present color-gradient LBM can correctly model and capture phase interfaces. x (cid:85) k
50 100 15000.20.40.60.811.2 (cid:85) b simulated (cid:85) g simulated (cid:85) r simulated (cid:85) b analytical (cid:85) g analytical (cid:85) r analytical Figure 2: The equilibrium density distributions of three di ff erent fluids for two separate static dropletsimmersed in a third fluid. A compound droplet, which consists of an inner droplet encapsulated by anotherimmiscible fluid, suspended in a third fluid, is simulated to assess whether the interfa-cial tensions are correctly modelled. The computational domain is taken as 160 × ff erent fluids, which are initialized as ρ r = , ρ g = ρ b = x − + ( y − ≤ R r ρ g = , ρ r = ρ b = R r < ( x − + ( y − ≤ R g ρ b = , ρ r = ρ g = R g = R r . This gives the initial condition that a compound droplet is located in thecenter of the computational domain. The interfacial tensions and the fluid viscositiesare all kept the same as those used in Section 3.1, and the periodic boundary conditionsare used in both x and y directions. According to the Young-Laplace’s law, when thesystem reaches the equilibrium state, the pressure di ff erence ∆ p across an interface isrelated to the interfacial tension σ by ∆ p = σ R , (26)where R is the radius of the interface curvature. Eq.(26) allows us to quantify themodeling accuracy of interfacial tensions through the relative error (cid:15) = (cid:12)(cid:12)(cid:12)(cid:12) ∆ p gb R g + ∆ p rg R r − (cid:16) σ gb + σ rg (cid:17)(cid:12)(cid:12)(cid:12)(cid:12) σ gb + σ rg × . (27) Table 1: The relative errors of interfacial tensions for various values of R r . R r
15 20 25 30 (cid:15)
Table 2: The maximum spurious velocities ( | u | max ) obtained with two di ff erent forcing methods for various R r . R r
15 20 25 30 | u | max × Present forcing method 1.68 1.69 1.70 1.71Forcing method of Leclaire et al. 3.15 3.28 3.27 3.29
Table 1 shows the relative errors of interfacial tensions for di ff erent values of R r . Allthe relative errors (cid:15) are below 1 . ff ects can also be realized by the forcingmethod of Leclaire et al.[50], i.e. Eqs.(4) and (5). It is of interest to compare thee ff ect of these two di ff erent forcing methods on spurious velocities. Table 2 showsthe maximum spurious velocities ( | u | max ) for various R r , where the values of | u | max aremagnified by 10 times. It is seen that the maximum spurious velocities are almostindependent of R r for either forcing method, and that the present spurious velocitiesare always smaller than those obtained with the forcing method of Leclaire et al. [50]. To assess the overall performance of the proposed model, we simulate the spreadingof a droplet between two other immiscible fluids. The computational domain is set tobe 160 ×
160 lattices. Initially, a red circular droplet with the radius R =
20 is placed inthe center of the computational domain, and the green and blue fluids are allocated tothe lower and upper halves of the computational domain outside the droplet. Periodicboundary conditions are used in both the x and y directions. Depending on the values10 (cid:2869) (cid:2016) (cid:2870) Red fluidBlue fluidGreen fluid (cid:1860) (cid:2869) (cid:1860) (cid:2870) (cid:1830)
Figure 3: The shape of a liquid lens at equilibrium. of the interfacial tensions, two di ff erent spreading phenomena can be observed, i.e.partial spreading and complete spreading.We first consider the partial spreading, where three interfacial tensions yield a Neu-mann’s triangle. In a partial spreading, the droplet can eventually reach a steady lensshape, which is often characterized by the lens length D and the heights h and h (seeFig. 3). The lens length and heights can be analytically given as [60, 52] D = (cid:118)(cid:117)(cid:117)(cid:117)(cid:116) A (cid:80) i = θ i (cid:16) θ i sin θ i − cos θ i (cid:17) , (28a) h i = D (cid:32) − cos θ i sin θ i (cid:33) with i = , , (28b)where A is the area of the red droplet; θ = ϕ rg and θ = ϕ rb are the equilibrium contactangles that can be calculated from Eq.(14). Four groups of interfacial tensions aresimulated with a constant σ gb of 0 .
01 but varying σ rb and σ rg , i.e., (a) σ rb = .
01 and σ rg = .
01, (b) σ rb = . σ rg = . σ rb = . σ rg = .
02, (d) σ rb = . σ rg = . .
1, and the finalfluid distributions are shown in Fig. 4. As expected, the droplet exhibits a lens shapein each of the cases considered, and the geometrical sizes ( D , h and h ) of the lensare case dependent. Based on the fluid distributions, we also quantify the geometricalsizes of the lens, and compare the simulated results with the analytical predictions fromEq.(28). It is seen in Table 3 that the simulated results (denoted by D s , h s and h s ) agreewell with the analytical predictions (denoted by D a , h a and h a ) with the relative errors(defined by E ( χ ) = | χ a − χ s | χ a × χ = D , h or h ) all around 1% exceptin the cases of small contact angles. The increased errors at small contact angles areattributed to the low resolution in sharp corners, which were also found by Jiang andTsuji[52].We then consider the complete spreading, where three interfacial tensions cannotyield a Neumann’s triangle. Two di ff erent cases of complete spreading are simulatedfor σ rg = .
01 and σ rg = .
015 at σ gb = σ rb = . σ rg = σ gb + σ rb in thefirst case, which corresponds to the critical state; whereas σ rg > σ gb + σ rb in the secondcase, which corresponds to the supercritical state. Fig. 5 shows the time evolution of11 a) (b) (c) (d)Figure 4: Final fluid distributions in the cases of partial spreading for (a) σ rb = . σ rg = .
01; (b) σ rb = . σ rg = . σ rb = . σ rg = . σ rb = . σ rg = . σ gb = . Case D a h a h a D s h s h s E ( D ) E ( h ) E ( h )(a) 55.34 15.97 15.97 53.59 15.92 15.92 3.17% 0.37% 0.38%(b) 65.18 8.73 18.81 61.79 8.79 18.97 5.21% 0.59% 0.82%(c) 45.84 13.23 22.92 45.39 13.35 22.51 0.99% 0.85% 1.80%(d) 50.99 6.83 25.49 50.93 6.77 25.22 0.13% 0.94% 1.09% the interface in both cases for a constant fluid viscosity of 0 .
05. We can see that in thecritical state, the red droplet sits exactly on the gb interface in the end; whereas in thesupercritical state, it bounces o ff the gb interface and rises up to the blue fluid.
4. Structure and stability of double droplets
Double emulsions have received considerable attention because of their poten-tial applications in food science, cosmetics, pharmaceuticals and medical diagnostics.Since emulsion properties and functions are related to the droplet geometry and com-position, it is of great importance, from a numerical point of view, to accurately predictthe topological structure of double droplets when dispersed in another immiscible fluid.
Consider a pair of equal-sized droplets, consisting of red and green fluids and ini-tially sitting next to each other, immersed in the third fluid (blue fluid). Based on thetheoretical analysis, Guzowski et al.[15] presented a stability diagram that describesthe possible topologies of double droplets and their transitions in terms of two ratiosof the interfacial tensions (see the left panel of Fig. 6). In the stability diagram, seventypical cases (represented by the solid points) are simulated to examine if the presentmodel is able to reproduce the correct morphologies of double droplets. These typicalcases are (i) σ gb σ rg = . σ rb σ rg = .
5, (ii) σ gb σ rg = σ rb σ rg =
1, (iii) σ gb σ rg = . σ rb σ rg = .
4, (iv) σ gb σ rg = . σ rb σ rg = .
5, (v) σ gb σ rg = σ rb σ rg =
1, (vi) σ gb σ rg = a)(b)Figure 5: Time evolution of the interface in the cases of complete spreading for (a) σ rg = .
01 and (b) σ rg = . σ gb = σ rb = . t = σ rb σ rg =
2, and (vii) σ gb σ rg = . σ rb σ rg = .
7, which cover all the possible morphologiesidentified by Guzowski et al.[15].The computational domain is taken to be [1 , × [1 , ρ r ( x , y ) = . + . (cid:34) R − √ ( x − . + ( y − . − R ) ξ (cid:35) , (29) ρ g ( x , y ) = . + . (cid:34) R − √ ( x − . + ( y − . + R ) ξ (cid:35) , (30) ρ b ( x , y ) = − ρ r ( x , y ) − ρ g ( x , y ) , (31)where the droplet radius R =
20 lattices. The periodic boundary conditions are used inboth the x and y directions. All the fluids are assumed to have equal viscosity of 0 .
1, andthe interfacial tension σ rg is fixed at 0 .
01. The simulations are run until an equilibriumstate is reached, and the equilibrium morphologies of double droplets for the sevencases are depicted in the right panel of Fig. 6. It is seen that seven di ff erent equilibriummorphologies are exhibited and they can be described as complete engulfing of greenfluid by red fluid (i), critical engulfing of green fluid by red fluid (ii), separate dispersionor non-engulfing (iii), kissing (iv), partial engulfing (v), critical engulfing of red fluidby green fluid (vi), and complete engulfing of red fluid by green fluid (vii). Thesesimulation results are consistent with the theoretical predictions by Guzowski et al.[15]. Among the seven morphologies shown in Fig. 6, the double droplets with partialengulfing morphology are often referred to as the Janus droplet. When the interfacial13 igure 6: Stability diagram representing possible morphologies of double droplets (left panel) and equi-librium shapes of the droplets for the typical cases marked in the stability diagram (right panel). The redlines represent the critical morphologies or the transitions between the regions of complete engulfing, partialengulfing and non-engulfing. tension between the constituent fluids is negligibly small, the Janus droplet forms aperfect circle, which is known as perfect Janus droplet (PJD) [15]. Di ff erentiatingfrom the PJD, the Janus droplet that does not exhibit a perfect circle is called as thegeneral Janus droplet (GJD).A Janus droplet, consisting of red and green fluids, is immersed in a static bluefluid. Fig. 7 shows the equilibrium geometries of a GJD and a PJD, as well as thecorresponding force balances at three-phase junctions. In this figure, R r , R g and R b are the curvature radii of the rb , gb and rg interfaces respectively ( kl interface refersto the interface between fluid k and fluid l ); θ r , θ g and θ b are the half of the centralangles subtended by the chord AB ; and d rg ( d gb ) is the distance between the centers O g and O r ( O b ). For a GJD, provided that four independent geometric parameters,e.g. R r , R g , R b and d rg , are given, one can analytically obtain all the other geometricparameters, including d gb , θ r , θ g and θ b , and the relative magnitudes of σ rb , σ gb and σ rg . Specifically, the half of the central angle θ g can be first calculated by θ g = arccos R g + d rg − R r R g d rg , (32)which is used to calculate the other two angles, θ r and θ b , and R b according to R g sin θ g = R b sin θ b = R r sin θ r , (33)14 (cid:3029) (cid:1841) (cid:3045) (cid:1841) (cid:3034) (cid:2016) (cid:3029) (cid:2016) (cid:3045) (cid:2016) (cid:3034) (cid:1838) (cid:3045)(cid:3034) (cid:1838) (cid:3029)(cid:3045) (cid:1838) (cid:3029)(cid:3034) (cid:2026) (cid:3029)(cid:3045) (cid:2026) (cid:3029)(cid:3034) (cid:2026) (cid:3045)(cid:3034) (cid:1844) (cid:3034) (cid:1844) (cid:3045) (cid:1844) (cid:3029) (cid:1827) (cid:1828) (cid:2026) (cid:3029)(cid:3045) (cid:2026) (cid:3029)(cid:3034) (cid:2026) (cid:3045)(cid:3034) (cid:2030) (cid:3045)(cid:3034) (cid:2030) (cid:3029)(cid:3034) (cid:2030) (cid:3029)(cid:3045) (a) General Janus droplet (cid:2026) (cid:3029)(cid:3034) (cid:2026) (cid:3029)(cid:3045) (cid:1844) (cid:3034) (cid:1844) (cid:3029) (cid:1841) (cid:3034) (cid:1841) (cid:3029) (cid:1838) (cid:3029)(cid:3045) (cid:1838) (cid:3045)(cid:3034) (cid:1838) (cid:3029)(cid:3034) (cid:2016) (cid:3029) (cid:2016) (cid:3034) (cid:2026) (cid:3029)(cid:3034) (cid:2026) (cid:3029)(cid:3045) (cid:2026) (cid:3045)(cid:3034) (cid:1372) (cid:882) (cid:1827) (cid:1828) (b) Perfect Janus dropletFigure 7: Equilibrium geometry and force balance at a three-phase junction for (a) a general Janus droplet(GJD) and (b) a perfect Janus droplet (PJD). and the distance d gb is then obtained as d gb = R g cos θ g + R b cos θ b . (34)Next, we determine the angles ϕ rg , ϕ rb and ϕ gb through the geometric relationship andthe Neumann’s triangle shown in Fig. 7(a). For example, when R g cos θ g < d gb and R g cos θ g ≥ d rg , these angles can be calculated by ϕ rb = π − θ r + θ g ϕ gb = π − θ r − θ b ϕ rg = π − ϕ rb − ϕ gb ; (35)and on the other hand, when R g cos θ g < d gb and R g cos θ g < d rg , we have ϕ rb = θ b + θ g ϕ gb = θ r − θ b ϕ rg = π − ϕ rb − ϕ gb . (36)Finally, one can obtain the relative magnitudes of the interfacial tensions by the law ofSines: σ rg sin ϕ rg = σ rb sin ϕ rb = σ gb sin ϕ gb . (37)In other words, all the interfacial tensions can be determined from Eq.(37) if one ofthem is also given, as we shall do below. 15y contrast, the geometry of a PJD is only determined by two areas of the dispersedfluids, i.e. A r and A g , and its analytical solution is given by R g = R r = (cid:114) A g + A r π R b = R g tan θ g A r A r + A g = θ g − sin θ g cos θ g + tan θ g (cid:16) π − θ g − sin θ g cos θ g (cid:17) π d gb = R g cos θ g + R b cos θ b . (38)The above equation suggests that one can obtain all the other geometric parameters,such as R b , θ g and d gb , if the area ratio A r A r + A g and R g are given.To test the accuracy of the present model for Janus droplets, we conduct two groupsof simulations with one for GJD and the other for PJD. The size of the computationaldomain is set as [1 , × [1 , ν k = .
1. Inthe GJD simulations, we select σ rg = . R r = R g =
80 and R b = d rg from 40 to 120 with an increment of 20. Using these parameters, wecan compute the geometric parameters d rg and d gb as well as the interfacial tensions σ gb and σ rg through Eqs. (32) to (37), which are presented in Table 4. We initialize thefluid distribution such that it follows the given and analytically computed geometricparameters described above, and assume that the circles for gb , rb , and rg interfacesare initially centered at (150 . , R g + . , R g + + d rg ) and (150 . , R g + + d gb ),respectively. In the PJD simulations, we select σ rb = σ gb = . σ rg = × − and Table 4: The interfacial tensions σ gb and σ rb and the distance d gb calculated from Eqs. (32)-(37) for GJDswith σ rg = . R r = R g =
80 and R b =
160 at di ff erent values of d rg . d rg σ gb σ rb d gb
20 0.01 0.02 –40 0.01699 0.01497 204.0860 0.01182 0.01261 201.8180 0.00797 0.00973 207.52100 0.00583 0.00812 216.63120 0.00449 0.00712 227.67140 0.00357 0.00643 240.00 R r = R g =
80, and vary the area fraction A r A r + A g from 0.1 to 0.5 with an incrementof 0.1. These parameters allow us to analytically compute all the other geometricparameters of a PJD, e.g. R b and d gb , which are listed in Table 5. We follow theanalytical geometric parameters to initialize the fluid distribution, and assume that thecircles for gb , rb and rg interfaces are initially located at (150 . , . . , . . , . + d gb ), respectively. In particular, we note that A r A r + A g = . R b → ∞ and d gb → ∞ , suggesting that the interface rg is theoretically a straight linelocated at y = .
5. 16 able 5: The geometric parameters R b and d gb calculated from Eq.(38) for PJDs with R r = R g =
80 atdi ff erent area fractions. A r A r + A g R b ∞ d gb ∞ All of the simulations are run until a steady state is reached. Figs. 8 and 9 showthe comparison between the analytical and simulated results for the GJDs and PJDs. Ineach of the figures, the analytical interface profiles are represented by the white linesof di ff erent patterns, while the red, green and blue fluids are indicated in red, green andblue, respectively. It is seen that our simulation results agree well with the analyticalones for various geometry configurations of GJD and PJD. (a) d rg =
40 (b) d rg =
60 (c) d rg =
80 (d) d rg =
100 (e) d rg = σ rg = . R r = R g =
80 and R b =
160 at di ff erent values of d rg . The analytical interface profiles are represented by the whitelines of di ff erent patterns, while the simulated red, green and blue fluids are indicated in red, green and blue,respectively.(a) A r A r + A g = . A r A r + A g = . A r A r + A g = . A r A r + A g = . A r A r + A g = . R r = R g =
80 at di ff erentarea fractions. The analytical interface profiles are represented by the white lines of di ff erent patterns, whilethe simulated red, green and blue fluids are indicated in red, green and blue, respectively. For double droplets immersed in a static matrix, the critical state occurs when thelargest interfacial tension equals the sum of the other two. As previously shown inFig. 6, the critical state of double droplets can be subdivided into the kissing state ((iv)in Fig. 6) and the critical engulfing state ((ii) or (vi) in Fig. 6). For the convenienceof description, we define the near-critical state as the state where the largest interfacial17ension is close to the sum of the other two. It is challenging to accurately simulatethe critical and near-critical states, where a slight inaccuracy in modeling could lead tosignificant simulation errors.To highlight the strength of the present model for critical scenarios, we consider thekissing / near-kissing states and the critical / near-critical engulfing states in a rectangulardomain of [1 , × [1 , / near-kissing states, a pair of equal-sizeddroplets with the radii of R r = R g =
60 are initially placed with a distance of d rg ,and they are symmetric with respect to the centerline y = .
5. The simulations areperformed for a constant σ rg of 0 .
01 but varying σ gb ( = σ rb ), which is varied aroundthe critical value of 0 .
005 with an increment of 2 × − . Note that the initial distance d rg depends on the value of σ gb , and is given by its analytical value in equilibrium as d rg = R r σ rg /σ gb if σ gb > . R r + R g =
120 otherwise . (39)In the critical / near-critical engulfing states, we consider a green droplet with R g =
80 entirely or partially engulfing a red droplet with R r =
60 for σ gb = σ rg = . σ rb is varied around the critical value of 0.02 with an increment of 2 × − . With theseparameters, we are able to analytically compute other geometric parameters, which aregiven by d rg = (cid:113) R r + R g − R r R g cos α , R b = R g sin θ g sin θ b and d gb = R g cos θ g − R b cos θ b for σ rb < .
02, and by d rg = R g − R r =
20 for σ rb ≥ .
02. Herein, cos α = σ rb σ gb , θ g = arccos R g + d gr − R r R g d gr and θ b = θ g + α . Again, we initialize the fluid distribution suchthat it follows the analytical geometric parameters. (cid:86) gb L r g TheoryThe present modelThe model of Fu et al.The model of Leclaire et al.
The exact kissing state (a) (cid:86) rb L r b TheoryThe present modelThe model of Fu et al.The model of Leclaire et al.
The critical engulfing state (b)Figure 10: (a) The interface length L rg as a function of σ gb in the kissing / near-kissing states; (b) the interfacelength L rb as a function of σ rb in the critical / near-critical engulfing states. In addition to the present model, we also use the model of Fu et al.[51] and themodel of Leclaire et al.[50] for the simulations. When the simulations reach the steady18tate, we quantify the interface lengths L rg in the kissing / near-kissing states and L rb inthe critical / near-critical engulfing states. Fig. 10 compares the simulated results fromthe present model with those from the model of Fu et al.[51] and the model of Leclaireet al.[50], and the analytical solutions. It is seen that for either kissing / near-kissingstates or critical / near-critical engulfing states, the simulated results from the presentmodel are in good agreement with the analytical solutions, while the simulated resultsfrom the other two models significantly deviate from the analytical solutions. Fig. 11shows the final fluid distributions obtained by the present model and the model of Fuet al.[51], in both kissing and critical engulfing states. Note that the fluid distributionfrom the model of Leclaire et al.[50] is not shown in the figure, since it produces almostthe same results as the model of Fu et al.[51]. Clearly, both critical states are correctlyreproduced by the present model but not by the model of Fu et al.[51]. These resultsindicate that the present model is advantageous to simulate critical state in ternaryfluids. (a) (b) (c) (d)Figure 11: The final fluid distributions obtained by (a) the present model and (b) the model of Fu et al. [51]for the kissing state, and by (c) the present model and (d) the model of Fu et al. [51] for the critical engulfingstate.
5. conclusions
A LB color-gradient model is proposed to simulate immiscible ternary fluids with afull range of interfacial tensions. An interfacial force formulation for N -phase ( N ≥ ff ectively reduce spurious velocities. A recoloring algorithm pro-posed by Spencer et al.[49] is applied to produce the phase segregation and ensurethe immiscibility of three di ff erent fluids, where a novel form of segregation parame-ters is proposed by considering the existence of Neumann’s triangle and the e ff ect ofequilibrium contact angle in three-phase junction. The model’s capability in capturinginterfaces and modeling interfacial tensions is first validated by the simulation of thetwo separate static droplets and the Young-Laplace test for a compound droplet. Theoverall performance of the model is then assessed by simulating the spreading of adroplet between two stratified fluids, and both the partial and complete spreadings arepredicted with satisfactory accuracy.Finally, the present model is used to study the stability and structure of doubledroplets in a static matrix over a wide range of interfacial tensions. By changing two19atios of the interfacial tensions, seven possible equilibrium morphologies are success-fully reproduced, which are consistent with the theoretical stability diagram by Gu-zowski et al.[15]. For various geometry configurations of general and perfect Janusdroplets, good agreemento between simulated results and analytical solutions showsthe present model is accurate when three interfacial tensions yield a Neumann’s trian-gle. In addition, we also simulate the near-critical and critical states of double droplets,which is challenging since the outcomes are very sensitive to the model accuracy. Itis found that the simulated results from the present model agree well with the ana-lytical solutions, while the simulated results from the existing color-gradient modelssignificantly deviate from the analytical solutions, especially in critical states. In sum-mary, the present work provides the first LB multiphase model that allows for accuratesimulation of ternary fluid flows with a full range of interfacial tensions. Acknowledgements
This work was supported by the National Natural Science Foundation of China(Nos. 51506168, 51711530130), the National Key Research and Development Projectof China (No. 2016YFB0200902), the China Postdoctoral Science Foundation (No.2016M590943), Guangdong Provincial Key Laboratory of Fire Science and Technol-ogy (No. 2010A060801010) and Guangdong Provincial Scientific and TechnologicalProject (No. 2011B090400518). Y. Yu was supported by the China Scholarship Coun-cil for one year study at the University of Strathclyde, UK. H. Liu gratefully acknowl-edges the financial supports from Thousand Youth Talents Program for DistinguishedYoung Scholars, the Young Talent Support Plan of Xi’an Jiaotong University.
References [1] A. S. Utada, E. Lorenceau, D. R. Link, P. D. Kaplan, H. A. Stone, D. A. Weitz,Monodisperse double emulsions generated from a microcapillary device, Science308 (2005) 537–541.[2] N. Bhatia, S. Pandit, S. Agrawal, D. Gupta, A review on multiple emulsions,International Journal of Pharmaceutical Erudition 3 (2013) 22–30.[3] H. Lamba, K. Sathish, L. Sabikhi, Double emulsions: Emerging delivery systemfor plant bioactives, Food and Bioprocess Technology 8 (2015) 709–728.[4] D. Chong, X. Liu, H. Ma, G. Huang, Y. L. Han, X. Cui, J. Yan, F. Xu, Ad-vances in fabricating double-emulsion droplets and their biomedical applications,Microfluidics and Nanofluidics 19 (2015) 1071–1090.[5] M. A. Augustin, Y. Hemar, Nano- and micro-structured assemblies for encapsu-lation of food ingredients, Chemical Society Reviews 38 (2009) 902–912.[6] V. B. Patravale, S. D. Mandawgade, Novel cosmetic delivery systems: an appli-cation update., International Journal of Cosmetic Science 30 (2008) 19–33.207] E. E. Ekanem, S. A. Nabavi, G. T. Vladisavljevi, S. Gu, Structured biodegrad-able polymeric microparticles for drug delivery produced using flow focusingglass microfluidic devices, ACS Applied Materials & Interfaces 7 (2015) 23132–23143.[8] B. Ahmed, D. A. Barrow, T. Wirth, Enhancement of reaction rates by segmentedfluid flow in capillary scale reactors, Advanced Synthesis & Catalysis 348 (2006)1043–1048.[9] L. Chen, Y. Li, J. Fan, H. K. Bisoyi, D. A. Weitz, Q. Li, Microshells: Pho-toresponsive monodisperse cholesteric liquid crystalline microshells for tunableomnidirectional lasing enabled by a visible lightdriven chiral molecular switch(advanced optical materials 9 / / navier-stokes model for the simulation of three-phase flows, Transport in Porous Media82 (2010) 463–483.[23] J. Kim, J. S. Lowengrub, Phase field modeling and simulation of three-phaseflows, Interfaces and Free Boundaries 7 (2005) 435–466.[24] J. Kim, Phase field computations for ternary fluid flows, Computer Methods inApplied Mechanics and Engineering 196 (2007) 4779–4788.[25] J. Kim, A generalized continuous surface tension force formulation for phase-field models for multi-component immiscible fluid flows, Computer Methods inApplied Mechanics and Engineering 198 (2009) 3105–3112.[26] J. Kim, Phase-field models for multi-component fluid flows, Communications inComputational Physics 12 (2012) 613–661.[27] H. Liu, Y. Ba, L. Wu, Z. Li, G. Xi, Y. Zhang, A hybrid lattice boltzmann andfinite di ff erence method for droplet dynamics with insoluble surfactants, Journalof Fluid Mechanics 837 (2018) 381–412.[28] H. Liu, A. J. Valocchi, Q. Kang, Three-dimensional lattice boltzmann model forimmiscible two-phase flow simulations, Physical Review E 85 (2012) 046309.[29] X. He, L.-S. Luo, A priori derivation of the lattice boltzmann equation, PhysicalReview E 55 (1997) R6333.[30] S. Chen, G. D. Doolen, Lattice boltzmann method for fluid flows, Annual reviewof fluid mechanics 30 (1998) 329–364.[31] C. K. Aidun, J. R. Clausen, Lattice-boltzmann method for complex flows, Annualreview of fluid mechanics 42 (2010) 439–472.[32] A. K. Gunstensen, D. H. Rothman, S. Zaleski, G. Zanetti, Lattice boltzmannmodel of immiscible fluids., Physical Review A 43 (1991) 4320–4327.[33] T. Reis, T. N. Phillips, Lattice boltzmann model for simulating immiscible two-phase flows, Journal of Physics A 40 (2007) 4033–4053.[34] X. Shan, H. Chen, Lattice boltzmann model for simulating flows with multiplephases and components, Physical Review E Statistical Physics Plasmas Fluids &Related Interdisciplinary Topics 47 (1993) 1815.[35] X. Shan, H. Chen, Simulation of nonideal gases and liquid-gas phase transitionsby the lattice boltzmann equation, Physical Review E 49 (1994) 2941–2948.[36] M. Sbragaglia, R. Benzi, L. Biferale, S. Succi, K. Sugiyama, F. Toschi, General-ized lattice boltzmann method with multirange pseudopotential., Physical ReviewE 75 (2007) 026702. 2237] M. Sbragaglia, R. Benzi, M. Bernaschi, S. Succi, The emergence of supramolec-ular forces from lattice kinetic models of non-ideal fluids: applications to therheology of soft glassy materials, Soft Matter 8 (2012) 10773–10782.[38] B. Dollet, A. Scagliarini, M. Sbragaglia, Two-dimensional plastic flow of foamsand emulsions in a channel: experiments and lattice boltzmann simulations, Jour-nal of Fluid Mechanics 766 (2015) 556–589.[39] M. R. Swift, E. Orlandini, W. Osborn, J. Yeomans, Lattice boltzmann simulationsof liquid-gas and binary fluid systems, Physical Review E 54 (1996) 5041.[40] Y. Wang, C. Shu, H. Huang, C. Teo, Multiphase lattice boltzmann flux solver forincompressible multiphase flows with large density ratio, Journal of Computa-tional Physics 280 (2015) 404–423.[41] Y. Wang, C. Shu, L. M. Yang, An improved multiphase lattice boltzmann fluxsolver for three-dimensional flows with large density ratio and high reynolds num-ber, Journal of Computational Physics 302 (2015) 41–58.[42] X. He, S. Chen, R. Zhang, A lattice boltzmann scheme for incompressible multi-phase flow and its application in simulation of rayleigh-taylor instability, Journalof Computational Physics 152 (1999) 642–663.[43] H. Chen, B. M. Boghosian, P. V. Coveney, A ternary lattice boltzmann model foramphiphilic fluids, Proceedings of The Royal Society A: Mathematical, Physicaland Engineering Sciences 456 (2000) 2043–2057.[44] O. Shardt, J. J. Derksen, S. K. Mitra, Simulations of janus droplets at equilibriumand in shear, Physics of Fluids 26 (2014) 106–114.[45] H. Liang, B. C. Shi, Z. Chai, Lattice boltzmann modeling of three-phase incom-pressible flows., Physical Review E 93 (2016) 013308.[46] C. Semprebon, T. Kruger, H. Kusumaatmaja, Ternary free energy lattice boltz-mann model with tunable surface tensions and contact angles, Physical ReviewE 93 (2016) 033305.[47] M. W¨ohrwag, C. Semprebon, A. M. Moqaddam, I. Karlin, H. Kusumaatmaja,Ternary free-energy entropic lattice boltzmann model with high density ratio,arXiv preprint arXiv:1710.07486 (2017).[48] Y. Shi, G. Tang, Y. Wang, Simulation of three-component fluid flows using themultiphase lattice boltzmann flux solver, Journal of Computational Physics 314(2016) 228–243.[49] T. Spencer, I. Halliday, C. M. Care, Lattice boltzmann equation method for mul-tiple immiscible continuum fluids, Physical Review E 82 (2010) 066701.[50] S. Leclaire, M. Reggio, J. Trepanier, Progress and investigation on lattice boltz-mann modeling of multiple immiscible fluids or components with variable densityand viscosity ratios, Journal of Computational Physics 246 (2013) 318–342.2351] Y. Fu, S. Zhao, L. Bai, Y. Jin, Y. Cheng, Numerical study of double emulsionformation in microchannels by a ternary lattice boltzmann method, ChemicalEngineering Science 146 (2016) 126–134.[52] F. Jiang, T. Tsuji, Estimation of threephase relative permeability by simulatingfluid dynamics directly on rockmicrostructure images, Water Resources Research53 (2017) 11–32.[53] H. Liu, Y. Zhang, Modelling thermocapillary migration of a microfluidic dropleton a solid surface, Journal of Computational Physics 280 (2015) 37–53.[54] Z. Guo, C. Zheng, B. Shi, Discrete lattice e ff ects on the forcing term in the latticeboltzmann method., Physical Review E 65 (2002) 046308.[55] I. Halliday, C. M. Care, A. P. Hollis, Improved simulation of drop dynamics ina shear flow at low reynolds and capillary number, Physical Review E 73 (2006)056708.[56] U. d’Ortona, D. Salin, M. Cieplak, R. B. Rybka, J. R. Banavar, Two-color nonlin-ear boltzmann cellular automata: surface tension and wetting, Physical Review E51 (1995) 3718.[57] M. Latva-Kokko, D. H. Rothman, Di ffff