Ternary free-energy entropic lattice Boltzmann model with high density ratio

A rigorous free energy model for ternary fluid flows with density ratio up to of order $O(10^3)$ is presented and implemented using the entropic lattice Boltzmann scheme. The model is thermodynamically consistent and allows a broad range of surface tension ratios, covering both partial wetting states where Neumann triangles are formed, and full wetting states where complete encapsulation of one of fluid components is observed. We further demonstrate that we can capture the bouncing, adhesive and insertive regimes for the binary collisions between immiscible droplets suspended in air. Our approach opens up a vast range of multiphase flow applications involving one gas and several liquid components.

Multiphase flows comprising of one gas and several liquid components are of considerable scientific interest due to their broad range of applications. The collision between oil and water droplets is a key ignition step in combustion engines, where the collision parameters can be varied to control the effective burning rate [1]. The presence of an immiscible crude oil layer on the sea surface alters the processes occurring during raindrop impact, with consequences for marine aerosol creation and oil spill dispersal [2]. In advanced oil recovery, considerable gain can be achieved by alternately displacing the oil by air and water in the so-called immiscible wateralternating-gas displacement process [3]. Infusing porous materials with lubricants results in composite surfaces, known as lubricant impregnated surfaces [4][5][6][7], with superior non-wetting and drag-reduction properties.
Despite the wide-ranging applications, suitable quantitative models for studying these phenomena are surprisingly still lacking. Most simulations to date have focussed on either single-component liquid-gas systems with high density ratio [8][9][10][11][12] or multicomponent flows with equal (or similar) density ratio [13][14][15][16][17][18][19]. In contrast, our aim here is to demonstrate an accurate and flexible model that can predict complex interfacial dynamics of ternary systems with significant ratio between the liquid and gas densities, up to of order O(10 3 ). This enables a new class of multiphase problems to be simulated, which are not previously possible. While we focus on one gas and two liquid components, the model can be extended to include more liquid components.
Our approach is based on the lattice Boltzmann method (LBM) [20,21], which has been shown to deliver reliable results, with quantitative agreement against experiments and other simulation methods, including on droplet dynamics [22][23][24], liquid phase separation [25,26] and flow through porous media [27]. In LBM interfacial forces can be implemented without explicit tracking of the interfaces, making it an elegant choice for studying mesoscopic interface dynamics in complex geometries.
Our key contribution over existing LBM models is a novel free energy functional that combines optimal equation of state for liquid-gas systems with double-well potentials to introduce multiple liquid components. The former, combined with the use of entropic lattice Boltzmann scheme [8], allows us to introduce significant density ratios, compared to other ternary free-energy LBM models [17,18]. The free energy formulation also ensures our model is thermodynamically consistent, unlike alternative approaches [14,28,29]. The capabilities of our new model are demonstrated using several static and dynamic problems. Firstly, we find excellent agreement between the numerical and analytical liquid-gas coexistence curves as a function of temperature, proving the thermodynamic consistency of the model. Secondly, we illustrate how the liquid-liquid and liquid-gas surface tensions can be flexibly tuned by simulating liquid lenses with varying Neumann angles. Finally, we simulate binary collisions between two immiscible droplets and show we capture many relevant features reported in experiments [30][31][32][33][34].
We introduce a free energy functional that consists of two parts, the bulk and interfacial contributions: The bulk free energy density f B is designed to allow three distinct minima, corresponding to one gas and two liquid components, as illustrated in Fig. 1. Ψ eos (ρ) can be derived from integrating the liquid-gas equation of state (EOS), p eos = ρ(dΨ eos /dρ)−Ψ eos , with coexisting liquidgas densities at ρ l and ρ g . For concreteness, here we use Carnahan-Starling EOS, but our approach is flexible, and in SI [35] we describe the implementation of Peng-Robinson and van der Walls EOS. For Carnahan-Starling arXiv:1710.07486v2 [physics.flu-dyn] 11 May 2018 Contour plot of the bulk free energy density fB as a function of two order parameters, ρ and φ. Three distinct minima exist, corresponding to a gas component at (ρg, 0), and two liquid components at (ρ l , +χ) and (ρ l , −χ). [36]: The constants C and Ψ 0 are chosen such that Ψ eos (ρ g ) = Ψ eos (ρ l ) = Ψ 0 , ensuring common tangent construction is met between all coexisting fluid phases. We use a = 0.037, b = 0.2 and R = 1. The critical temperature is T c = 0.3373 a bR , and the temperature T governs the liquid-gas density ratio.
The second and third terms in Eq. (2) have the form of double well potentials with C l1 and C l2 the relative concentrations of the two liquid components. Established works on critical phenomena show that such form is universal to describe the physics of continuous phase transitions close to the critical point [37], including for fluid mixtures. Away from the critical point, additional terms may be needed. However, so long as the details of the equation of state of the fluid mixtures is not important for the problem at hand, a large body of work in diffuse interface models for binary fluids has shown a double well potential is sufficient to capture interfacial dynamics with excellent agreement against both theory and experimental results, such as for droplet dynamics [22,23] and coarsening in phase separation [25,26]. This is the case for the examples studied here, and the double well potentials are therefore chosen as the simplest model possible.
Each double well term has two minima at C l# = 0 (component absent) and C l# = 1 (present). We also define the relative concentration of the gas phase as C g = (ρ − ρ l )/(ρ g − ρ l ), which is 0 for ρ = ρ l and 1 for ρ = ρ g . Given the constraint C g + C l1 + C l2 = 1, there are two independent order parameters: the density ρ and the phase field φ. The relative concentrations are related to the density and phase field via C l1 = 1 , with χ a constant scaling parameter for φ. Our free energy functional has three minima at (ρ g , 0), (ρ l , +χ) and (ρ l , −χ).
For the interfacial free energy density, f I , all three terms in Eq. (3) are necessary because there are three independent surface tensions in ternary systems. Upon expanding C l# in terms of ρ and φ, we can rewrite f I as We vary the λ parameters in Eq. (2) and κ parameters in Eq. (3) to tune the surface tensions and interfacial widths of the three fluid interfaces. The continuum equations of motion for the fluid are the continuity, Navier-Stokes, and Cahn-Hilliard equations: v is the fluid velocity and η is the dynamic viscosity that depends on the local density and phase field. For simplicity, we employ constant mobility parameter M , though in general it can depend on the local density and phase field [38]. The thermodynamics of ternary fluids, described by the free energy functional F in Eq. (1), enter the equations of motion through the chemical potentials, µ ρ = δF/δρ| T,φ and µ φ = δF/δφ| T,ρ , and the pressure tensor, ∇ · P = ρ∇µ ρ + φ∇µ φ . To solve the equations of motion, we introduce two sets of distribution functions in our LBM scheme, evolving the density and phase field. For the former, we employ the entropic lattice Boltzmann method, augmented with an exact-difference forcing term [8]. For the latter we use a standard BGK scheme [21]. We provide the details of our LBM implementation in SI [35], including the expressions for the chemical potentials and pressure tensor.
To demonstrate the accuracy and broad range of surface tension ratios allowed in our model, we simulate a liquid lens, where a droplet of liquid 1 is suspended at the interface between liquid 2 and the gas phase, as shown in Fig. 2(a). We show the profiles of the density ρ and phase field φ across the liquid lens configuration in Fig.  2(b). At the interface between liquids 1 and 2, ρ remains constant at ρ l , while φ transforms smoothly between −χ and χ. Both ρ and φ vary at the interface between any of the liquids and the gas.
At equilibrium, force balance between the surface tensions at the three phase contact line leads to a distinct set of angles known as the Neumann angles. Mathematically, γ 12 / sin(θ 3 ) = γ 2g / sin(θ 2 ) = γ 1g / sin(θ 1 ). To test this relation we vary the value of κ 3 , while keeping λ 3 = 3.125 × κ 3 and other simulation parameters (see caption of Fig. 2) constant. For all simulations shown here, we also set the kinematic viscosity ν = η/ρ = 0.167, and mobility parameter M = 0.5. Fig. 2(c) shows the Neumann angles calculated in two different ways. Firstly, we measure the Neumann angles geometrically (diamond symbol) from our liquid lens simulations. Secondly, we use Laplace pressure tests to independently measure surface tensions for all permutations of the interfaces (see SI [35]), and subsequently compute the expected Neumann angles (square symbol). The agreement is excellent, with typical deviations of < 3 • . Similar agreement is observed upon varying other parameters. In addition to partial wetting states, where the Neumann triangle is formed, our model allows simulations of full wetting states. To demonstrate this, in SI [35], we present simulation results of two droplets where γ 1g + γ 12 < γ 2g . The droplets are initialised such that they are just touching each other. As dictated by thermodynamics, the simulation shows that the liquid 2 droplet becomes fully encapsulated by the liquid 1 droplet.
A wide range of density ratios can be simulated by tuning the temperature T in the equation of state, Eq. (4). Fig. 3 shows the coexistence curve for the Carnahan-Starling EOS. The left (right) branches correspond to the gas (liquid) densities. Good agreement is obtained between the analytical solution from Maxwell construction (line) and the numerical results (dots). The lowest temperature we can robustly simulate is T = 0.61T c , corresponding to a numerical density ratio of O(10 3 ). In SI [35] we also show that high density ratios can be achieved with Peng-Robinson and van der Walls EOS.
We now present simulation results of collisions between two immiscible droplets. In comparison to the more  commonly studied problem of collisions between miscible droplets of the same materials (e.g. [8,9]), the collision outcomes for immiscible droplets are much richer. Here we show three regimes observed in experiments: bouncing, adhesive and insertive collisions, and their transitions. To our best knowledge, this is the first time they have been simulated using LBM. We will focus on generic features of the drop collisions. Systematic studies, including parameter matching against experiments, will be presented elsewhere.
We first consider bouncing collision. Fig. 4(a) shows an experimental example where the two droplets are water and diesel oil [31]. As the droplets collide (columns ii and iii), we observe compression in the drop shapes parallel to the collision direction and radial expansions perpendicular to the collision direction. This is followed by retraction in the radial direction (column iv), and if there is sufficient kinetic energy, the two droplets bounce off and become separated (column v). Our simulations show this sequence is ubiquitous for head-on bouncing collisions. Fig. 4(b) shows one such case at We 1 = We 2 = 20.8 and Re 1 = Re 2 = 72.0, where We i = ρ i V 2 r D i /γ ig , Re i = ρ i V r D i /η i , and V r is the relative droplet velocity.
Here the two droplets have symmetric properties and we use γ 12 /γ 2g = 1.33. We set T = 0.65T c for the rest of the paper, corresponding to a density ratio of ρ l /ρ g 150. For the cases shown here, the results do not sensitively depend on the density ratio beyond ρ l /ρ g ∼ 100. This is illustrated explicitly in SI [35] by comparing the results in Fig. 4 to those obtained using T r = 0.61 (ρ l /ρ g 1000).
By reducing the droplets' velocities, we observe a transition from bouncing to adhesive collision, shown in Fig.  4(c) for We 1 = We 2 = 5.6 and Re 1 = Re 2 = 36.0. Qualitatively the initial collision dynamics is similar between rows (b) and (c). However, at column (iv) there is not enough kinetic energy for the droplets to detach. Subsequently the compound droplet oscillates until it relaxes to its equilibrium configuration, determined by the Neumann triangle. Animations of the drop collisions in Fig. 4(b) and (c) are provided as SI [35]. Adhesive collision between two immiscible droplets with similar liquidgas surface tension has been observed experimentally for diesel and ethanol droplets [32].
A powerful advantage of our model is that it covers a wide range of surface tension ratios. We can now consider the asymmetric case where the liquid-gas surface tension of droplet 2 is considerably larger than that for droplet 1, yet it does not correspond to the full wetting state. Fig. 5(b-e) shows the case where γ 12 /γ 2g = 0.54 and γ 1g /γ 2g = 0.49, with normalised spreading parameter S = 1 − (γ 1g + γ 12 )/γ 2g = −0.029. In agreement with experimental observations [30], we observe a transition between adhesive and insertive collisions by varying the impact velocities.
An experimental example of insertive collision is shown in Fig. 5(a) for water and n-hexadecane [30]. For comparison, Fig. 5(b) and (c) show the typical dynamical sequence observed in our simulations, with We 1 = 16.4, We 2 = 6.1 and Re 1 = Re 2 = 45.0. Upon collision, the composite droplet expands radially (column iii), followed by contraction in the radial direction (column iv) and elongation in the collision axis (column v). The oscillation between the prolate and oblate shapes can sustain several periods (videos in SI [35]), accompanied by the propagation of the three-phase contact line until the high surface tension droplet is fully encapsulated (column vi).
The transition from insertive to adhesive collision can be induced by decreasing the droplets' velocities. In Fig.  5(d) and (e), we present the case where We 1 = 4.4, We 2 = 2.2 and Re 1 = Re 2 = 27.0. Initially the contact line propagates to cover the high surface tension droplet as the composite droplet oscillates between the prolate and oblate shapes (videos in SI [35]). Since the kinetic energy is insufficient to drive full encapsulation, the contact line eventually recedes and the droplet relaxes to its equilibrium shape (column vi). In SI [35] we further show the critical velocity for the transition between insertive and adhesive collisions increases as the normalised spreading parameter becomes more negative.
To conclude, we presented a strategy for modelling ternary multiphase multicomponent flows by combining a novel free energy formulation and the use of entropic LBM scheme. Our approach allows significant density ratios, up to of order O(10 3 ), and a broad range of surface tension ratios, covering both partial and full wetting states, to be simulated. These flexibilities open up a number of applications, which are not previously possible.
As an example, we demonstrated the bouncing, adhesive and insertive regimes for binary collisions between immiscible droplets. Our method can meet the gap in systematic computational work for such collision dynamics, to complement the rich body of existing experimental studies [30][31][32][33][34]. Other applications are numerous, including drop impact on immiscible liquid film [2], advanced oil recovery [3], and liquid impregnated surfaces [4][5][6][7].
Here we have assumed the liquids to have the same density. This is justifiable in most water-oil-gas systems where the liquid-liquid density ratio is several orders of magnitude smaller than the liquid-gas density ratio. A useful future extension is to allow all density ratios to be varied independently. Our model can also be generalised to include more liquid components, by introducing additional double well potential and gradient terms in the bulk and interfacial free energy densities respectively. Another key avenue for future work is the interactions between ternary flows and complex solid surfaces. Our model is compatible with various approaches to introduce wetting boundary conditions [39][40][41].