Phase-field-based lattice Boltzmann model for immiscible incompressible N-phase flows

TL;DR

This paper introduces an efficient lattice Boltzmann model for simulating complex immiscible incompressible N-phase flows, validated against benchmarks and capable of handling multi-phase interactions.

Contribution

The paper develops a novel LB model based on phase field theory that reduces gradient calculations and ensures thermodynamic consistency for N-phase flows.

Findings

Accurate simulation of static droplets and liquid lens spreading.

Successful modeling of droplet collision dynamics for four phases.

Effective handling of complex interactions among multiple immiscible fluids.

Abstract

In this paper, we develop an efficient lattice Boltzmann (LB) model for simulating immiscible incompressible $N$-phase flows $(N \geq 2)$ based on the Cahn-Hilliard phase field theory. In order to facilitate the design of LB model and reduce the calculation of the gradient term, the governing equations of the $N$-phase system are reformulated, and they satisfy the conservation of mass, momentum and the second law of thermodynamics. In the present model, $(N-1)$ LB equations are employed to capture the interface, and another LB equation is used to solve the Navier-Stokes (N-S) equations, where a new distribution function for the total force is delicately designed to reduce the calculation of the gradient term. The developed model is first validated by two classical benchmark problems, including the tests of static droplets and the spreading of a liquid lens, the simulation results show that the current LB model is accurate and efficient for simulating incompressible $N$-phase fluid flows. To further demonstrate the capability of the LB model, two numerical simulations, including dynamics of droplet collision for four fluid phases and dynamics of droplets and interfaces for five fluid phases, are performed to test the developed model. The results show that the present model can successfully handle complex interactions among $N$ ($N \geq 2$) immiscible incompressible flows.

Tables

Table 1. Table 1: The equilibrium contact angles θ 1 subscript 𝜃 1 \theta_{1} and θ 2 subscript 𝜃 2 \theta_{2} with different surface tension ratios.

Surface tension	Present LB model		Analytical solutions		Relative errors
$\left({\sigma }_{12}:{\sigma }_{13}:{\sigma }_{23}\right)$	${\theta }_1$	${\theta }_2$	${\theta }_1$	${\theta }_2$	${\theta }_1$	${\theta }_2$
$1:\frac{4}{3}:1$	${84.6}^{\circ }$	${50.0}^{\circ }$	${83.5}^{\circ }$	${52.5}^{\circ }$	$1.2\%$	$0.2\%$
$1:1:1$	${60.7}^{\circ }$	${60.7}^{\circ }$	${60.0}^{\circ }$	${60.0}^{\circ }$	$1.2\%$	$1.2\%$
$0.6:0.6:1$	${34.1}^{\circ }$	${34.1}^{\circ }$	${34.0}^{\circ }$	${34.0}^{\circ }$	$1.5\%$	$1.5\%$

Table 2. Table 2: The equilibrium length d 𝑑 d and height ( h 1 , h 2 ) subscript ℎ 1 subscript ℎ 2 (h_{1},h_{2}) with different surface tension ratios.

Surface tension	Present LB model			Analytical solutions			Relative errors
$\left({\sigma }_{12}:{\sigma }_{13}:{\sigma }_{23}\right)$	$d$	$h_1$	$h_2$	$d$	$h_1$	$h_2$	$d$	$h_1$	$h_2$
$1:\frac{4}{3}:1$	$73.0$	$33.2$	$16.3$	$73.8$	$33.0$	$16.5$	$1.1\%$	$0.6\%$	$1.2\%$
$1:1:1$	$82.0$	$23.6$	$23.6$	$80.7$	$23.3$	$23.3$	$1.7\%$	$1.3\%$	$1.3\%$
$0.6:0.6:1$	$114.0$	$17.5$	$17.5$	$113.9$	$17.2$	$17.2$	$0.1\%$	$1.8\%$	$1.8\%$

Equations

$$V=V_1+V_2+\cdots +V_N.$$

$${\tilde{\gamma }}_i=\frac{1}{{\tilde{\rho }}_i},1\le i\le N;\Gamma =\sum _{i=1}^N{\tilde{\gamma }}_i;$$

$$c_i=\frac{V_i}{V}=\frac{M_i/{\tilde{\rho }}_i}{M_i/{\rho }_i}=\frac{{\rho }_i}{{\tilde{\rho }}_i}={\tilde{\gamma }}_i{\rho }_i,1\le i\le N;\sum _{i=1}^N{c}_i=\frac{{\sum }_{i=1}^N{V}_i}{V}=1;\rho =\sum _{i=1}^N{\rho }_i.$$

$$\frac{\partial {\phi }_i(\varphi )}{\partial t}+\mathbf{u}\cdot \nabla {\phi }_i(\varphi )=-\nabla \cdot {\mathbf{J}}_{ai},1\le i\le N-1,$$

$$\nabla \cdot \mathbf{u}=0.$$

$${\mathbf{J}}_{ai}=-{\tilde{m}}_i(\varphi )\nabla C_i,1\le i\le N-1,$$

$$\sum _{j=1}^{N-1}\frac{\partial {\phi }_j}{\partial {\varphi }_i}{C}_j=\frac{\partial W(\varphi ,\nabla \varphi )}{\partial {\varphi }_i}-\nabla \cdot \frac{\partial W(\varphi ,\nabla \varphi )}{\partial (\nabla {\varphi }_i)},1\le i\le N-1,$$

$$\rho (\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u})+\tilde{\mathbf{J}}\cdot \nabla \mathbf{u}=-\nabla p+\nabla \cdot \mathbf{S},$$

$$\tilde{\mathbf{J}}=\sum _{i=1}^{N-1}\left(1-\frac{N}{\Gamma }{\tilde{\gamma }}_i\right){\mathbf{J}}_{ai}.$$

$$\mathbf{S}=\mu (\varphi )\mathbf{D}(\mathbf{u})-\sum _{i=1}^{N-1}\nabla {\varphi }_i\otimes \frac{\partial W(\varphi ,\nabla \varphi )}{\partial \nabla {\varphi }_i},$$

$$\nabla \cdot \mathbf{u}=0,$$

$$\rho \left(\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}\right)+\tilde{\mathbf{J}}\cdot \nabla \mathbf{u}=-\nabla p+\nabla \cdot [\mu (\varphi )\mathbf{D}(\mathbf{u})]-\sum _{i=1}^{N-1}\nabla \cdot \left(\nabla {\varphi }_i\otimes \frac{\partial W}{\partial \nabla {\varphi }_i}\right),$$

$$\sum _{j=1}^{N-1}\frac{\partial {\phi }_i}{\partial {\varphi }_j}\left(\frac{\partial {\varphi }_j}{\partial t}+\mathbf{u}\cdot \nabla {\varphi }_j\right)=\nabla \cdot [{\tilde{m}}_i(\varphi )\nabla C_i],1\le i\le N-1,$$

$$\begin{array}{rl}W(\varphi ,\nabla \varphi )& =\sum _{i,j}^{N-1}\frac{{\lambda }_{ij}}{2}\nabla {\varphi }_i\cdot \nabla {\varphi }_j+\frac{{\beta }^2}{{\eta }^2}H(\varphi )\\ & =\sum _{i,j}^{N-1}\frac{{\lambda }_{ij}}{2}\nabla {\varphi }_i\cdot \nabla {\varphi }_j+\frac{{\beta }^2}{{\eta }^2}\sum _{k=1}^N{c}_k^2(1-c_k{)}^2,\end{array}$$

$${\phi }_i(\varphi )={\rho }_i-{\rho }_N=\frac{1}{2}({\tilde{\rho }}_i-{\tilde{\rho }}_N)+\frac{1}{2}({\tilde{\rho }}_i+{\tilde{\rho }}_N){\varphi }_i,1\le i\le N-1,$$

$$\sum _{i=1}^{N-1}(L_i^{kl}{)}^2{\lambda }_{ii}+\sum _{j=1}^{N-1}\sum _{i=1}^{j-1}2L_i^{kl}{L}_j^{kl}{\lambda }_{ij}=\frac{9}{2}\frac{{\eta }^2}{{\beta }^2}{\sigma }_{kl}^2,1\le k<l\le N,$$

$$L_{i}^{kN}=\left\{\begin{array}[]{l}1,\quad\quad\,if\ i=k,\\ \frac{\tilde{\rho}_{N}}{\tilde{\rho}_{i}+\tilde{\rho}_{N}},if\ i\neq k.\end{array}\right.1\leq k<l=N,$$

$$L_{i}^{kl}=\left\{\begin{array}[]{l}\frac{\tilde{\rho}_{k}}{\tilde{\rho}_{k}+\tilde{\rho}_{N}},\ \ \quad\,if\ i=k,\\ -\frac{\tilde{\rho}_{l}}{\tilde{\rho}_{l}+\tilde{\rho}_{N}},\,\quad if\ i=l,\\ 0,\qquad\qquad otherwise,\end{array}\right.1\leq k<l<N.$$

$$\nabla \cdot \mathbf{u}=0,$$

$$\rho \left(\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}\right)+\tilde{\mathbf{J}}\cdot \nabla \mathbf{u}=-\nabla p+\nabla \cdot [\mu (\varphi )\mathbf{D}(\mathbf{u})]-\sum _{i,j=1}^{N-1}\nabla \cdot \left({\lambda }_{ij}\nabla {\varphi }_i\nabla {\varphi }_j\right)+\mathbf{G}(\mathbf{x},t),$$

$$\frac{\partial {\varphi }_i}{\partial t}+\mathbf{u}\cdot \nabla {\varphi }_i=m_i{\nabla }^2{C}_i,1\le i\le N-1,$$

$$\left\{\begin{array}[]{l}\rho_{i}=\frac{1}{\Gamma}+\sum_{j=1}^{N-1}\left(\delta_{ij}-\frac{\tilde{\gamma}_{j}}{\Gamma}\right)\left[\frac{1}{2}(\tilde{\rho}_{j}-\tilde{\rho}_{N})+\frac{1}{2}(\tilde{\rho}_{j}+\tilde{\rho}_{N})\phi_{j}\right],\quad 1\leq i\leq N,\\ c_{i}(\phi)=\tilde{\gamma}_{i}\rho_{i}=\frac{\tilde{\gamma}_{i}}{\Gamma}+\sum_{j=1}^{N-1}\left(\tilde{\gamma}_{i}\delta_{ij}-\frac{\tilde{\gamma}_{i}\tilde{\gamma}_{j}}{\Gamma}\right)\left[\frac{1}{2}(\tilde{\rho}_{j}-\tilde{\rho}_{N})+\frac{1}{2}(\tilde{\rho}_{j}+\tilde{\rho}_{N})\phi_{j}\right],\quad 1\leq i\leq N,\\ C_{i}(\phi)=\frac{\delta W}{\delta\phi_{i}}=-\sum_{j=1}^{N-1}\lambda_{ij}\nabla^{2}\phi_{j}+\frac{\beta^{2}}{\eta^{2}}h_{i}(\phi),\quad 1\leq i\leq N-1,\\ \tilde{\mathbf{J}}=-\sum_{i=1}^{N-1}\left(1-\frac{N}{\Gamma}\tilde{\gamma}_{i}\right)\frac{\tilde{\rho}_{i}+\tilde{\rho}_{N}}{2}m_{i}\nabla\left[-\sum_{j=1}^{N-1}\lambda_{ij}\nabla^{2}\phi_{j}+\frac{\beta^{2}}{\eta^{2}}h_{i}(\phi)\right],\\ \rho=\sum_{k=1}^{N}\rho_{k}=\frac{N}{\Gamma}+\sum_{i=1}^{N-1}\left(1-\frac{N}{\Gamma}\tilde{\gamma}_{i}\right)\left[\frac{1}{2}(\tilde{\rho}_{i}-\tilde{\rho}_{N})+\frac{1}{2}(\tilde{\rho}_{i}+\tilde{\rho}_{N})\phi_{j}\right],\end{array}\right.$$

$$h_i(\varphi )=\frac{\partial H}{\partial {\varphi }_i}=\frac{{\tilde{\rho }}_i+{\tilde{\rho }}_N}{2}\sum _{k=1}^N\left({\tilde{\gamma }}_k{\delta }_{ki}-\frac{{\tilde{\gamma }}_k{\tilde{\gamma }}_i}{\Gamma }\right)c_k(1-c_k)(1-2c_k),1\le i\le N-1.$$

$$\nabla \cdot \mathbf{u}=0,$$

$$\frac{\partial (\rho \mathbf{u})}{\partial t}+\nabla \cdot (\rho \mathbf{uu})=-\nabla p+\nabla \cdot [\mu (\varphi )\mathbf{D}(\mathbf{u})]-\nabla \cdot (\tilde{\mathbf{J}}\mathbf{u})+{\mathbf{F}}_s+\mathbf{G},$$

$$\frac{\partial {\varphi }_i}{\partial t}+\nabla \cdot ({\varphi }_i\mathbf{u})=m_i{\nabla }^2{C}_i,1\le i\le N-1,$$

$$\begin{array}{rl}-{\lambda }_{ii}\nabla \cdot (\nabla {\varphi }_i\nabla {\varphi }_i)& =-{\lambda }_{ii}{\nabla }_{\beta }({\nabla }_{\alpha }{\varphi }_i{\nabla }_{\beta }{\varphi }_i)\\ & =-{\lambda }_{ii}[{\nabla }_{\beta }({\nabla }_{\alpha }{\varphi }_i){\nabla }_{\beta }{\varphi }_i+{\nabla }_{\alpha }{\varphi }_i{\nabla }_{\beta }^2{\varphi }_i]\\ & =-{\lambda }_{ii}[{\nabla }_{\alpha }(\frac{1}{2}{\nabla }_{\beta }{\varphi }_i{\nabla }_{\beta }{\varphi }_i)+{\nabla }_{\alpha }{\varphi }_i{\nabla }_{\beta }^2{\varphi }_i]\\ & =-{\lambda }_{ii}[\nabla (\frac{1}{2}{\left|\nabla {\varphi }_i\right|}^2)+\nabla {\varphi }_i{\nabla }^2{\varphi }_i],\end{array}$$

$$\begin{array}{rl}& -[{\lambda }_{ij}\nabla \cdot (\nabla {\varphi }_i\nabla {\varphi }_j)+{\lambda }_{ji}\nabla \cdot (\nabla {\varphi }_j\nabla {\varphi }_i)]\\ =& -[{\lambda }_{ij}{\nabla }_{\beta }({\nabla }_{\alpha }{\varphi }_i{\nabla }_{\beta }{\varphi }_j)+{\lambda }_{ji}{\nabla }_{\beta }({\nabla }_{\alpha }{\varphi }_j{\nabla }_{\beta }{\varphi }_i)]\\ =& -{\lambda }_{ij}[{\nabla }_{\beta }({\nabla }_{\alpha }{\varphi }_i){\nabla }_{\beta }{\varphi }_j+{\nabla }_{\alpha }{\varphi }_i{\nabla }_{\beta }^2{\varphi }_j+{\nabla }_{\beta }({\nabla }_{\alpha }{\varphi }_j){\nabla }_{\beta }{\varphi }_i+{\nabla }_{\alpha }{\varphi }_j{\nabla }_{\beta }^2{\varphi }_i]\\ =& -{\lambda }_{ij}\{{\nabla }_{\alpha }[({\nabla }_{\beta }{\varphi }_i)({\nabla }_{\beta }{\varphi }_j)]+{\nabla }_{\alpha }{\varphi }_i{\nabla }_{\beta }^2{\varphi }_j+{\nabla }_{\alpha }{\varphi }_j{\nabla }_{\beta }^2{\varphi }_i\}\\ =& -{\lambda }_{ij}\nabla [(\nabla {\varphi }_i)\cdot (\nabla {\varphi }_j)]-{\lambda }_{ij}[\nabla {\varphi }_i{\nabla }^2{\varphi }_j+\nabla {\varphi }_j{\nabla }^2{\varphi }_i],\end{array}$$

$$\begin{array}{rl}{\mathbf{F}}_s& =\sum _{i=1}^{N-1}{C}_i\nabla {\varphi }_i\\ & =\sum _{i=1}^{N-1}\left[-\sum _{j=1}^{N-1}{\lambda }_{ij}{\nabla }^2{\varphi }_j+\frac{{\beta }^2}{{\eta }^2}{h}_i(\varphi )\right]\nabla {\varphi }_i\\ & =-\sum _{i,j=1}^{N-1}{\lambda }_{ij}{\nabla }^2{\varphi }_j\nabla {\varphi }_i-\sum _{i=1}^{N-1}\frac{{\beta }^2}{{\eta }^2}{h}_i(\varphi )\nabla {\varphi }_i\\ & =-\sum _{i,j=1}^{N-1}{\lambda }_{ij}{\nabla }^2{\varphi }_j\nabla {\varphi }_i-\sum _{i=1}^{N-1}\frac{{\beta }^2}{{\eta }^2}\frac{\partial H}{\partial {\varphi }_i}\nabla {\varphi }_i\\ & =-\sum _{i,j=1}^{N-1}{\lambda }_{ij}{\nabla }^2{\varphi }_j\nabla {\varphi }_i-\sum _{i=1}^{N-1}\frac{{\beta }^2}{{\eta }^2}\nabla H,\end{array}$$

$$\frac{\partial f_k}{\partial t}+{\mathbf{c}}_k\cdot \nabla f_k={\Omega }_k+G_k$$

$$f_k(\text{x}+{\text{c}}_k\Delta t,t+\Delta t)-f_k(\text{x},t)={\int }_0^{\Delta t}{\Omega }_k(\text{x}+{\text{c}}_k{t}^{\prime },t+t^{\prime })d{t}^{\prime }+{\int }_0^{\Delta t}{G}_k(\text{x}+{\text{c}}_k{t}^{\prime },t+t^{\prime })d{t}^{\prime }.$$

$${\Omega }_k=-\frac{1}{{\tau }^{\prime }}(f_k-f_k^{eq}),$$

$$f_k(\text{x}+{\text{c}}_k\Delta t,t+\Delta t)-f_k(\text{x},t)=\frac{\Delta t}{2}\left[{\Omega }_k(\text{x}+{\text{c}}_k\Delta t,t+\Delta t)+{\Omega }_k(\text{x},t)\right]+{\int }_0^{\Delta t}{G}_k(\text{x}+{\text{c}}_k{t}^{\prime },t+t^{\prime })d{t}^{\prime }.$$

Full text

1 Introduction

$N$-phase flow problems ($N\geq 2$) are ubiquitous in nature and engineering, such as the emulsion and foam formation in the microfluidic devices, where one or more fluid phases are dispersed in another continuous phase [1, 2, 3, 4], the enhanced oil recovery, the geological CO2 sequestration in depleted oil/gas reservoirs [4, 5, 6], as well as many daily phenomena such as rain drops, spraying of pesticides, bubbles in water, etc. In many cases, the number of fluid phases is greater than or equal to $3$, and $N$-phase flow problems are usually accompanied by the droplet/bubble generation, coalescence, and breakup. There are roughly three types of methods for studying these complicated $N$-phase flow problems, namely theoretical approach, experimental approach and numerical simulation. Due to limitations of theoretical and experimental methods, numerical simulation plays an important role in the study of $N$-phase flows. However, numerical simulation of such $N$-phase problems ($N\geq 3$) is still very challenging because of the inherent nonlinearities, complex topological changes and the complexity of moving interfaces. Due to the importance and complexity of the $N$-phase problem ($N\geq 3$), our main focus in the present paper will be on situations involving more than two fluid phases.

Compared to the numerous researches of two-phase fluid flows [7, 8, 9, 10, 11, 12, 13], there are fewer theoretical and numerical studies on fluid flows that involves three or more phases. Nonetheless, some efforts have been made to develop efficient numerical methods for dealing with $N$-phase flows, including the volume of fluid (VOF) method [14, 15], level set method [16], front tracking method [17], smoothed particle hydrodynamics method [18, 19], diffuse interface method (phase field method) [20, 21, 22, 23, 5, 24], and lattice Boltzmann (LB) method [27, 28, 29, 30, 31, 32, 33, 34, 6, 35]. In this work, we will focus on the phase-field-based LB method for immiscible incompressible $N$-phase flows ($N\geq 3$).

In the phase field method, the thickness of the interface between the two immiscible flows is supposed to be very small but nonzero. Then different phases can be characterized by phase variables (order parameters) which vary continuously across thin interfacial layers. One can derive a set of governing equations for the whole computational domain from the free energy, where the order parameters satisfy an advection-diffusion equation (usually the Cahn-Hilliard equation), and this equation is coupled with the Navier-Stokes equations. With the Cahn-Hilliard phase field approach, Boyer et al. [23] proposed a diffuse interface model for the study of three immiscible component incompressible viscous flows. The effects of different forms for the bulk free energy have been investigated. Kim et al. [24] proposed a phase field model for the immiscible incompressible ternary fluid flows with interfaces. In these models, the pairwise surface tensions $\sigma_{ij}$ are decomposed into three positive phase-specific surface-tension coefficients. However, this decomposition will encounter some difficulties when $N\geq 4$, as the number of pairwise surface tensions $\frac{N(N-1)}{2}$ would be greater than the number of mixing energy density coefficients $N$, which leads to an overdetermined system [25]. To overcome this shortcomings, Kim et al. [25] proposed the generalized continuous surface tension force model for multi-component fluid flows. Boyer and his collaborator [26] have proposed a generalized model for the $N$-phase mixtures using the consistency principle which can be described as that the $N$-phase model exactly coincide with the classical two-phase model when only two phases are present in the system. Recently, Dong [4] have derived a physical formulation and a numerical algorithm for the mixture of N ($N\geq 2$) immiscible incompressible fluids in the thermodynamics framework. In his model, the mixture velocity is the volume-averaged velocity, and the mixing energy density coefficients involved in the $N$-phase model can be determined based on the pairwise surface tensions among the N fluids.

As a mesoscopic level method, the LB method has made rapid progress since its appearance in the late 1980s due to its simplicity, scalability on parallel computers, and ease to handle complex geometries [44, 45]. Based on different physical perspectives, the existing LB models for multiphase flows can be generally classified into four categories: color-gradient model [36], pseudo-potential model [37], free-energy model [38], and phase-field-based model [39, 40, 41, 42, 43]. These methods have gained great success in the study of two-phase flow problems, and they have been extended to simulate immiscible ternary fluids. For instance, Lamura et al. [27] proposed a LB model to simulate oil-water-surfactant mixtures based on a Ginzburg-Landau free energy. However, this model can only be used to simulate ternary flows where an amphiphile phase is located at oil-water interface, and cannot be extended to arbitrary ternary flows. Chen et al. [29] presented a LB model for amphiphilic ternary fluid dynamics, which can be regarded as an extension of the pseudo-potential model. His model suffers from the inherent defects in the pseudo-potential model, such as high spurious velocities and lack of flexibility in adjusting surface tension [34]. Leclaire et al. [32] developed a LB model for the simulation of multiple immiscible fluids based on the improved color-gradient model, where three subcollision operators are also applied. We note that thermodynamic consistency, which is distinctly important in the physical formulation, are lack in the above three LB models. Recently, an alternative LB ternary model was proposed by Liang et al. [6] based on the Cahn-Hilliard phase field theory, which provides a firm physical foundation on the dynamics of the interfaces among three fluids. Shi et al. [35] extended the LB flux solver which was originally proposed for simulating incompressible flows of binary fluids based on two-component Cahn-Hilliard model to three-component fluid flows, while their models are limited to fluids with a very small density difference (no more than $1.25$). Although Shi et al. [35] pointed out that the model can be extended to $N$-phase $(N\geq 4)$ fluid flows by introducing the free energy $W(\phi,\nabla\phi)=\int_{\Omega}\left[0.25\sum_{i=1}^{N}\phi_{i}^{2}(1-\phi_{i}^{2})+\frac{\epsilon^{2}}{2}\sum_{i=1}^{N}|\nabla\phi_{i}|^{2}\right]\,d\mathbf{x}$, where $\phi_{i}$ is the order parameter and $\epsilon$ is a small parameter denoting the interface thickness, we note that this free energy expression is unreasonable as it does not involve the interaction between different phases. In fact, this interaction is related to the multiple pairwise surface tensions. Besides, the author does not give a numerical example when $N\geq 4$. Therefore, all the mentioned LB models are limited to the case of $N\leq 3$.

To our best knowledge, there is still no available work on the construction of LB model for fluid flows with $N\geq 4$. In this work, we develop an efficient LB model for simulating immiscible incompressible $N$-phase flows $(N\geq 2)$ based on the Cahn-Hilliard phase field theory. The proposed model has several distinct features. Firstly, the governing equations of the $N$-phase system proposed by Dong [4] are reformulated, and they also conserve the mass, momentum and second law of thermodynamics. In such a sense, the reformulated system is thermodynamically consistent. Secondly, the velocity in this model is the volume-averaged mixture velocity and is divergence free. Thirdly, the mixing energy density coefficients involve the interaction between different phases and can be determined by the pairwise surface tensions among the N fluids. Fourthly, in the present model, $(N-1)$ LB equations are employed to capture the interface, and another LB equation is used to solve the N-S equations. Fifthly, we introduce a new distribution function for the total force to reduce the calculation of the gradient term, and the governing equations can be recovered correctly from the present LB model through Chapmann-Enskog analysis.

The rest of present paper is organized as follows. In Sec. II, we describe the physical formulation of $N$-phase flows, and a phase-field-based LB model for immiscible incompressible $N$-phase flows is given in Sec. III. In Sec. IV, we use two classic numerical examples to validate our model. It is found that the numerical results agree well with analytical solutions. In Sec. V, some numerical applications for $N$-phase fluid flows ($N\geq 4$) are conducted. Finally, some conclusions are given in Sec. VI.

2 Physical formulation

In this section, we will describe the physical formulation of immiscible incompressible $N$-phase flows ($N\geq 2$) based on the mass conservation, momentum conservation, the second law of thermodynamics, and Galilean invariance. We refer to [4] for detailed derivations of this system.

Consider an incompressible system consisting of N immiscible Newtonian fluids. Let $\Omega$ denote the flow domain, $\tilde{\rho}_{i}$ ($1\leq i\leq N$) represent constant densities of these $N$ pure fluids, $V$ is an arbitrary control volume with the mass $M$, $V_{i}$ is the volume of fluid $i$ before mixing with the mass $M_{i}$, and $\tilde{\mu}_{i}$ expresses their constant dynamic viscosities. Here we assume that there is no volume loss or increase during the $N$-phase mixing process, which leads to

[TABLE]

For convenience, we introduce the following auxiliary parameters

[TABLE]

$\phi_{i}$ ($1\leq i\leq N-1$) are employed to denote ($N-1$) independent order parameters, and $-1\leq\phi_{i}\leq 1$. $\phi=(\phi_{1},\phi_{2},\cdots,\phi_{N-1})$ represents the vector of $(N-1)$ phase field functions, $c_{i}$ and $\rho_{i}$ ($1\leq i\leq N$) are used to represent the volume fraction and density of fluid $i$ within the mixture, and $\rho$ is introduced to denote the average density of the mixture. Then they satisfy the following relations [4, 46, 47]

[TABLE]

The mass conservation of $N$-phase flows can be given by the following ($N-1$) mass balance equations [4, 46, 47]

[TABLE]

where $\varphi_{i}(\phi)$ is defined as $\varphi_{i}\equiv\rho_{i}-\rho_{N}$ ($1\leq i\leq N-1$), and $\mathbf{u}$ is the volume averaged mixture velocity and can be rigorously proved that it is divergence free [4]

[TABLE]

$\mathbf{J}_{ai}$ is the relative differential flux between fluids $i$ and $N$, and it can be written as

[TABLE]

based on the second law of thermodynamics, where $\tilde{m}_{i}(\phi)\geq 0$ is the mobility. $C_{i}$ ($\leq i\leq N-1$) are ($N-1$) effective chemical potentials determined by the linear system

[TABLE]

where $W(\phi,\nabla\phi)$ represents the free energy density function.

The momentum conservation can be described by [4]

[TABLE]

where $\tilde{\mathbf{J}}$ is the flux and satisfies

[TABLE]

$p$ is the pressure, and $\mathbf{S}$ represents a trace-free stress tensor which is determined as

[TABLE]

based on the second law of thermodynamics, where $\mathbf{D(u)}=\nabla\mathbf{u}+\nabla\mathbf{u}^{T}$, $\mu(\phi)=\sum_{k=1}^{N}\tilde{\mu}_{k}c_{k}(\phi)$ is the viscosity, $\otimes$ denotes the tensor product.

With these constitutive relations, we obtain the following phase field system for $N$-phase flows [4]

[TABLE]

based on the mass conservation [Eq. (4)] and momentum conservation [Eq. (8)]. Once the free energy density $W(\phi,\nabla\phi)$ and the functions $\varphi_{i}(\phi)$ ($1\leq i\leq N-1$) are specified, we can compute all the other quantities in this model. Following [4], we chose $W(\phi,\nabla\phi)$ and $\varphi_{i}(\phi)$ as

[TABLE]

[TABLE]

where the constant $\beta^{2}$ represents the characteristic energy scale, and the constant $\eta>0$ denotes the characteristic scale of the interfacial thickness. The constants $\lambda_{ij}$ ($1\leq i,j\leq N-1$) denote the mixing energy density coefficients, and they can be determined by the pairwise surface tensions $\sigma_{kl}$ ($1\leq k<l\leq N$),

[TABLE]

where $\lambda_{ij}$ should be symmetry, and $L_{i}^{kl}$ ($1\leq i,j\leq N-1$) are determined as

[TABLE]

or

[TABLE]

With the expressions (12) and (13), Eqs. (11a)-(11c) will reduce to

[TABLE]

where we have taken into account an external body force $\mathbf{G}(\mathbf{x},t)$ in the momentum equation, and the constant $m_{i}$ ($m_{i}>0$) is the mobility.

For the convenience of calculation, the parameters $\rho_{i}$, $c_{i}$, $C_{i}$, $\tilde{\mathbf{J}}$, $\rho$ can be rewritten as

[TABLE]

where $\delta_{ij}$ is the Kronecker delta function, and $h_{i}(\phi)$ is given as

[TABLE]

Note that the $N$-phase system consisting of Eqs. (17a)-(17c) is complicated due to excessive gradient calculation in momentum equation, and it is not easy to solve this system directly using LB method. In order to facilitate the design of an efficient LB model, we introduce the following system of equations

[TABLE]

where $\mathbf{F}_{s}$ is the surface tension which can be given as $\mathbf{F}_{s}=\sum_{i=1}^{N-1}C_{i}\nabla\phi_{i}$. Next, we will show that Eq. (17b) and Eq.(20b) are equivalent.

Consider the term $-\sum_{i,j=1}^{N-1}\nabla\cdot\left(\lambda_{ij}\nabla\phi_{i}\nabla\phi_{j}\right)$. For $1\leq i=j\leq N-1$,

[TABLE]

where the first term $-\lambda_{ii}\nabla(\frac{1}{2}\left|\nabla\phi_{i}\right|^{2})$ can be absorbed into the pressure gradient term in momentum equation. For $1\leq i<j\leq N-1$, we have

[TABLE]

where the relation $\lambda_{ij}=\lambda_{ji}$ has been used, and the first term $-\lambda_{ij}\nabla[(\nabla\phi_{i})\cdot(\nabla\phi_{j})]$ can also be absorbed into the pressure gradient term in momentum equation.

Note that the surface tension $\mathbf{F}_{s}$ can be transformed as follows

[TABLE]

where the last term $-\sum_{i=1}^{N-1}\frac{\beta^{2}}{\eta^{2}}\nabla H$ can also be absorbed into the pressure gradient term. In contrast to Eqs. (21)-(23), we conclude that $-\sum_{i,j=1}^{N-1}\nabla\cdot\left(\lambda_{ij}\nabla\phi_{i}\nabla\phi_{j}\right)$ and $\mathbf{F}_{s}$ are equivalent except for the terms which are absorbed into $-\nabla p$. In addition, with the relation of $\tilde{\mathbf{J}}\cdot\nabla\mathbf{u}=\nabla\cdot(\tilde{\mathbf{J}}\mathbf{u})-(\nabla\cdot\tilde{\mathbf{J}})\mathbf{u}$ and the mass conservation equation $\frac{\partial\rho}{\partial t}+\mathbf{u}\cdot\nabla\rho=-\nabla\cdot\tilde{\mathbf{J}}$ [4], Eq. (17b) can be converted to Eq. (20b).

Thus, we have proved that the two immiscible $N$-phase systems [Eqs. (17a)-(17c) and Eqs. (20a)-(20c)] are equivalent. But Eqs. (20a)-(20c) are much simpler from the perspective of LB model construction. Since the chemical potential $C_{i}$ is a quantity that must be calculated, we only need to calculate gradient terms for $(N-1)$ times when calculating the surface tension $\mathbf{F}_{s}$, which is much simpler than calculating the term $-\sum_{i,j=1}^{N-1}\nabla\cdot(\lambda_{ij}\nabla\phi_{i}\nabla\phi_{j})$, especially when $N$ is relatively large. Simultaneously, the reformulated $N$-phase system Eqs. (20a)-(20c) can also satisfy the conservations of mass, momentum and the second law of thermodynamics. Therefore, next we will construct the LB model based on the system of Eqs. (20a)-(20c).

3 LB model for immiscible incompressible N-phase flows

In this section, we will present a LB model for simulating the incompressible $N$-phase flows. The governing equations Eqs. (20a)-(20c) can be regarded as the coupling between the incompressible N-S equations and the $(N-1)$ C-H equations. Next, we will take two sets of LB models to solve the N-S equations and the C-H equations separately.

3.1 LB model for the Navier-Stokes equations

To obtain the evolution equation, we integrate the following discrete velocity Boltzmann equation

[TABLE]

along a characteristic line $\mathbf{c}_{k}$ over a time interval $\Delta t$ [50, 51], then we get

[TABLE]

where $f_{k}(\mathbf{x},t)$ denotes particle distribution function for fluids $i$ with velocity $\mathbf{c}_{k}$ at position $\mathbf{x}$ and time $t$, $G_{k}(\mathbf{x},t)$ represents the force term, $\Omega_{k}(\mathbf{x},t)$ is the collision operator which can be approximated by

[TABLE]

where $\tau^{\prime}$ denotes the relaxation time and $f_{k}^{eq}(\mathbf{x},t)$ represents the equilibrium distribution function.

The trapezoidal rule is used to integrate the collision term, then Eq. (25) becomes

[TABLE]

Let $\bar{f}_{k}=f_{k}-\frac{\Delta t}{2}\Omega_{k}$, we can transform Eq. (27) into

[TABLE]

where $\tau_{g}=\frac{2\tau^{\prime}+\Delta t}{2\Delta t}$ denotes the dimensionless relaxation time, and $\bar{f}_{k}$ satisfies $\sum_{k}\bar{f}_{k}=\sum_{k}f_{k}$, and $\sum_{k}\mathbf{c}_{k}\bar{f}_{k}=\sum_{k}\mathbf{c}_{k}f_{k}$.

Through Taylor expansion of $G_{k}(\textbf{x}+\textbf{c}_{k}t^{\prime},t+t^{\prime})$ and neglecting the terms of order $O(\Delta t^{2})$, the last term in the right hand of Eq. (28) can be transformed into [51]

[TABLE]

where $D_{k}=\partial t+\mathbf{c}_{k}\cdot\nabla$.

The LB evolution equation with the BGK collision operator for the N-S equations can be written as

[TABLE]

For the term $D_{k}G_{k}(\mathbf{x},t)$ in Eq. (30), one can take different discretization schemes to deal with it. When the up-wind scheme is used, the evolution equation becomes

[TABLE]

To remove the implicitness, we introduce a modified particle distribution function,

[TABLE]

With some simple operations, the explicit evolution equation can be derived,

[TABLE]

where $g_{k}^{eq}$ denotes the new equilibrium distribution function which satisfies $g_{k}^{eq}=f_{k}^{eq}$.

To recover the N-S equations (20a)-(20b) through Chapman-Enskog analysis, the equilibrium distribution function $g_{k}^{eq}$ is delicately designed as

[TABLE]

with

[TABLE]

where $\rho_{0}$ is a constant, $\omega_{k}$ and $\mathbf{c}_{k}$ denote the weighting coefficient and the discrete velocity, and $c_{s}$ represents the speed of sound. We would like to point that the present model is based on the D$d$Q$q$ lattice with $q$ velocity directions in $d$-dimensional space. The values of $\mathbf{c}_{k}$ and $\omega_{k}$ depend on the lattice model used. In D$1$Q$3$ model, $\{\mathbf{c}_{k}\}=c(0,1,-1)$, $\omega_{0}=2/3$, $\omega_{1,2}=1/6$, $c_{s}=c/\sqrt{3}$, where $c=\Delta x/\Delta t$ , with $\Delta x$ and $\Delta t$ representing the spacing and time step, respectively; In the D$2$Q$9$ model, $\omega_{k}$ is chosen as $\omega_{0}=4/9$, $\omega_{1-4}=1/9$, $\omega_{5-8}=1/36$, $c_{s}=c/\sqrt{3}$, and $\mathbf{c}_{k}$ is defined as

[TABLE]

In the D$3$Q$15$ model, $\omega_{0}=2/9$, $\omega_{1-6}=1/9$, $\omega_{7-14}=1/72$, $c_{s}=c/\sqrt{3}$, and $\mathbf{c}_{i}$ is defined as

[TABLE]

In order to accurately recover the macroscopic equations, the design of the force distribution function is very important. In the present model, the distribution function for the total force is given by

[TABLE]

where $\mathbf{F}$ is the total force and is expressed as

[TABLE]

It can be rigorously proved that the present LB model can correctly recover Eqs. (20a)-(20b) through the Chapman-Enskog analysis (see A for the details) with the following fluid kinematic viscosity

[TABLE]

In the LB method, the macroscopic quantities, e.g. fluid velocity $\mathbf{u}$ and pressure $p$, can be calculated from the moment of the distribution function [41, 6],

[TABLE]

[TABLE]

Remark 1. The distribution function for the total force can also be designed as

[TABLE]

where $\mathbf{F}$ is redefined as

[TABLE]

We would like to point out that the first force distribution function [Eq. (36)] is better because we don’t have to calculate the gradient term $-\nabla\cdot(\tilde{\mathbf{J}}\mathbf{u})$. Therefore, in the following simulations, we will adopt Eq. (36) as the force distribution function.

3.2 LB model for the Cahn-Hilliard equations

To complete the modeling of incompressible $N$-phase flows, the N-S equations should be coupled with the system of $(N-1)$ C-H equations. Here the evolution equations for the C-H equations can be written as

[TABLE]

where $h_{k}^{i}(\textbf{x},t)$ denotes the distribution function of order parameter $\phi_{i}$, $\tau_{i}$ represents the non-dimensional relaxation time for $i$-phase which is related to the mobility, $h_{k}^{i,eq}(\textbf{x},t)$ is the local equilibrium distribution function, which is introduced as [58, 57, 59]

[TABLE]

where $\eta_{i}$ denotes an adjustable parameter that controls the mobility for $i$-fluid. $R_{k}(\textbf{x},t)$ is the source term and is defined as [41]

[TABLE]

In our simulations, the temporal derivative in Eq. (45) is calculated by the first-order Eulerian scheme, i.e.,

[TABLE]

In the present model, the order parameter $\phi_{i}$ is calculated by taking the zeroth moment of the order distribution function,

[TABLE]

Then, according to mass conservation, the density $\rho$ can be projected on the basis of $\phi_{i}$,

[TABLE]

It is shown using the Chapman-Enskog analysis [41] that the C-H equations [Eqs. (20c)] can be recovered with second-order accuracy and the mobility can be determined by

[TABLE]

In addition, to compute the first-order spatial derivatives and the Laplacian operators in the force term $G_{k}$, surface tension $\mathbf{F}_{s}$ and chemical potentials $C_{i}$ $1\leq i\leq N-1$, the following isotropic schemes are adopted if not specified [60]:

[TABLE]

where $\zeta$ is any macroscopic quantity. This scheme is referred to as isotropic central scheme, and it not only has a second-order accuracy in space, but also conserve the global mass of a $N$-phase system [41].

4 Model validation

In this section, we employ two classical benchmark problems, including tests of static droplets and the spreading of a liquid lens between the other two components, to demonstrate the capability and accuracy of the LB model. Some detailed comparisons of numerical results and the analytical solutions are conducted. Here we indicate that in all the simulations, the D$2$Q$9$ lattice structure is adopted, and the grid resolution test has been conducted. It is demonstrated that the grids used in the present work are accurate enough to give grid-independent results.

4.1 Static droplets

A basic test of static droplets is first performed to verify the present LB model. Initially, two droplets with radius $R=20$ are placed in a lattice domain of $NX\times NY=300\times 100$ and the periodic boundary condition is applied at all boundaries. The initial volume fractions are given by

[TABLE]

where $(x_{c_{1}},y_{c_{1}})$ and $(x_{c_{2}},y_{c_{2}})$ represent the coordinate of circular droplet and are fixed as $(x_{c_{1}},y_{c_{1}})=(50,50)$, $(x_{c_{2}},y_{c_{2}})=(150,50)$. In all the simulations, $c=\delta x=\delta t=1.0$. Some other parameters are given as $\rho_{1}:\rho_{2}:\rho_{3}=20:1:5$, $\tau_{g}=\tau_{1}=\tau_{2}=0.8$, $\sigma_{12}=\sigma_{13}=\sigma_{23}=0.01$, $m_{1}=m_{2}=0.001$. When choosing $\beta$ and $\eta$ in the free energy density expression, we follow such a relation

[TABLE]

where $\sigma_{min}=\min{\sigma_{ij}}$ is the minimum value among the $\frac{1}{2}N(N-1)$ surface tensions [4]. Here $\eta$ is set to $\eta=\sqrt{2}$, then $\beta=\sqrt{0.06}$. Fig. 1 shows the steady distributions of volume fractions obtained by the present LB model. It can be observed that the present LB model can accurately preserve the profiles of volume fractions with initial configurations for this $3$-phase problem. In order to give quantitative comparisons, we also plot the volume fractions along the centerline ($x=50$, $0\leq y\leq 300$) in Fig. 2. From this figure, one can observe that the numerical results of the volume fractions agree well with the theoretical results.

Besides, we also consider a $4$-phase problem of static droplets. In this test, three droplets with radius $R=20$ are placed in the same domain of $NX\times NY=300\times 100$. The initial volume fractions are expressed as

[TABLE]

with $(x_{c_{1}},y_{c_{1}})=(50,50)$, $(x_{c_{2}},y_{c_{2}})=(150,50)$ and $(x_{c_{3}},y_{c_{3}})=(250,50)$. Some other parameters are given as $\rho_{1}:\rho_{2}:\rho_{3}:\rho 4=20:1:10:5$, $\tau_{g}=\tau_{1}=\tau_{2}=\tau_{3}=0.8$, $\sigma_{12}=\sigma_{13}=\sigma_{14}=\sigma_{23}=\sigma_{24}=\sigma{34}=0.01$, $m_{1}=m_{2}=m_{3}=0.001$. The rest of the parameters are set as previous case. As shown in Fig. 3, the present LB model can also accurately preserve the profiles of volume fractions with initial configurations for $4$-phase problems. To further give a detailed comparison, we also plot the volume fractions along the center line and compare the results between analytical results and numerical results. The results show that the developed LB model is accurate enough to simulate multiphase problems.

Next, we will show that our LB model for $N$-phase flows can satisfy the reduction-consistent property. Here, the reduction-consistent property means that when one phase of the $N$-phase system disappears, the designed model can still be used to simulate the remaining ($N-1$)-phase system. Consider the case where the volume fraction of the third phase is [math] in the previous $4$-phase example. The remaining parts can be viewed as a $3$-phase system. For convenience, the $4$th-phase in the former case is marked as the $3$th-phase. Fig. 5 depicts the comparisons of volume fractions obtained from our $3$-phase model, $4$-phase model with one phase disappeared, and theoretical values along the centerline. It can be seen that the results from our $3$-phase model and $4$-phase model with one phase disappeared are almost identical, and they are also in good agreement with the theoretical values. Thus, we conclude that our LB model can satisfy the reduction-consistent property.

4.2 Spreading of a liquid lens

Here we consider a droplet located at the interface between the other two immiscible fluids. Under the influence of surface tension, the droplet will form a lens when the system reaches equilibrium. According to Neumann’s law [48], the following relationships are satisfied between the contact angles and the surface tensions at the equilibrium state,

[TABLE]

where $\theta_{i}$ ($i=1,2$) denote the contact angles as shown in Fig. 6. The relationship between the lens area $A$, its length $d$ (the distance between two triple junctions) and the contact angles is

[TABLE]

After some simple geometric transformations, we can derive the following expression

[TABLE]

Initially, a circular droplet with the radius $R=30$ is placed in the middle of the computational domain with $NX\times NY=150\times 150$. The initial volume fractions are given by

[TABLE]

where $(x_{c},y_{c})$ is the coordinate of circular droplet.

In this test, the periodic boundary conditions are employed in the $x$-direction, and the halfway bounce-back boundary conditions are enforced on the top and bottom walls. The density ratio of the three fluids is set to $\rho_{1}:\rho_{2}:\rho_{3}=10:1:5$. Some other physical parameters are given as $\tau_{g}=\tau_{1}=\tau_{2}=0.8$, $\sigma_{12}:\sigma_{13}:\sigma_{23}=1:\frac{4}{3}:1$, $1:1:1$, and $0.6:0.6:1$, $m_{1}=m_{2}=0.1$, $\eta=\sqrt{2}$, $\beta=\sqrt{0.06}$. It is shown in Fig. 7 that the droplets will undergo some deformations due to the effect of surface tension, and eventually form distinct equilibrium configurations. The shapes of the liquid interface in this work are in good agreement with the previous results [49, 23]. In order to give quantitative results, we also measured the contact angles $\theta_{1}$ and $\theta_{2}$ and compared them with the analytical solutions (see Table 1), where the numerical solutions $\theta_{1}$ and $\theta_{2}$ are calculated by

[TABLE]

From Table 1, one can observe that the numerical contact angles $\theta_{1}$ and $\theta_{2}$ are in good agreement with analytical solutions, where the maximum relative errors are no more than $1.5\%$. In addition, we also give a comparison between the numerical solutions and the analytical solutions of the length $d$ and the height $(h_{1},h_{2})$, see Table 2. It is shown that the maximum relative errors of $d$, $h_{1}$ and $h_{2}$ shall not exceed $1.8\%$, which means that the present LB model can accurately capture the interface of the incompressible three-phase flow without considering gravity.

Furthermore, now we consider the effect of gravity on the equilibrium configurations. To generate the gravitational effects, a body force, $\mathbf{G}=(0,\rho g)$, is added to the momentum equation, where $g$ is the gravitational acceleration. In this test, we vary the value of the gravitational acceleration while fixing all the other physical parameters. Due to the large deformation of the interface, the computational domain is set to $NX\times NY=400\times 200$. The physical parameters are given as $\rho_{1}:\rho_{2}:\rho_{3}=3:1:6$, $\tau_{g}=\tau_{1}=\tau_{2}=0.8$, $\sigma_{12}:\sigma_{13}:\sigma_{23}=1:1:1$, $m_{1}=m_{2}=0.1$. As shown in Fig. 8, when the gravity is relatively small and surface tension is dominant, the droplet forms a lens on the water surface. As the gravity increases, the droplet becomes flatter [see Fig. 8]. When the gravity is large enough and dominates over the surface tension, it will forms a puddle [see Fig. 8]. The current results are qualitatively consistent with the de Gennes theory [52]. To give a quantitative comparison, we measured asymptotic thickness $H$ as a function of gravity and compared it with the de Gennes theory, where the asymptotic puddle thickness can be given as [52, 53]

[TABLE]

Fig. 9 depicts the asymptotic thickness as a function of gravity between the present results and the de Gennes theory. One can observe that the present results agree well with the de Gennes theory when the gravity is large ($g\geq 0.5$ for this case). However, when the gravity is small, a large discrepancy between the numerical results and theoretical results, because the expression (59) is not valid for the case that $g$ is too small [4].

In addition, we also consider the effect of surface tension on the equilibrium configurations. In this test, we fix the gravitational acceleration $g$ and the surface tensions $\sigma_{12}$, $\sigma_{23}$ to $g=1\times 10^{-5}$, $\sigma_{12}=\sigma_{23}=0.01$ and vary the surface tension $\sigma_{13}$. The computational domain is set to $NX\times NY=300\times 150$. The other parameters are the same as the previous test. As shown in Fig. 10, the droplet thickness strongly depend on the surface tension. As the surface tension increases, the droplet thickness gradually increases. To quantitatively describe the relationship between asymptotic thickness and surface tension, we give a comparison of the results between the current simulation and the de Gennes theory [52] (see Fig. 11). From Fig. 11, one can conclude that $H$ and $\sigma_{13}$ increase simultaneously, and the simulation results are in good agreement with the de Gennes theory.

Therefore, all the results of this section, especially qualitative and quantitative comparisons with de Gennes theory, show that the current LB model is accurate and efficient for simulating the motion of incompressible multiphase systems.

5 Numerical applications

To further demonstrate the capabilities and performance of the LB model we proposed for $N$-phase flows ($N\geq 2$), two numerical applications involving four to five fluid phases are carried out. These applications include the dynamics of droplet collision for four fluid phases and dynamics of droplets and interfaces for five fluid phases. The simulation results show that the present LB model can be used to study the interactions among multiple types of fluid interfaces.

5.1 Dynamics of droplet collision for four fluid phases

In this subsection, we study a dynamic problem involving four fluid phases as another test for validating the developed LB method. This phenomenon is very common in engineering and industrial problems, as well as in nature, such as in spraying processes, in micro-reactors, and in kitchen [55, 56, 4]. The problem configuration is shown in Fig. 12, where the upper part of the domain is filled with phase $3$, and the lower part is filled with phase $4$. One droplet, referred to as phase $1$, is located inside phase $3$. Simultaneously, a bubble of the same size, denoted as phase $2$, is located in phase $4$. Then, the droplet and the bubble start to move from rest under gravity (buoyancy). In this simulation, the computational domain is chosen as $NY\times NX=320\times 200$. Periodic boundary conditions are applied in the horizonal direction, and the halfway bounce-back boundary conditions are enforced on the top and bottom walls. The initial volume fractions are given by

[TABLE]

where $(x_{c_{1}},y_{c_{1}})$ and $(x_{c_{2}},y_{c_{2}})$ are the center of two droplets, and $y_{0}=\frac{1}{2}NY$ is the initial surface position between phase $3$ and phase $4$. Both of phase $1$ and phase $2$ have a diameter of $D=60$. Different values of surface tension $\sigma_{13}$ are investigated. The other parameters are fixed as $\rho_{1}:\rho_{2}:\rho_{3}:\rho_{4}=6:1:3:2$, $\sigma_{12}=\sigma_{14}=\sigma_{23}=\sigma_{24}=\sigma_{34}=0.01$, $\tau_{g}=\tau_{1}=\tau_{2}=\tau_{3}=0.8$, $m_{1}=m_{2}=m_{3}=0.1$, $\eta=\sqrt{2}$, $\beta=\sqrt{0.06}$, $g=10^{-5}$, $(x_{c1},y_{c1})=(160,280)$, $(x_{c2},y_{c2})=(160,40)$. Fig. 13 shows the time sequence of evolution of the fluid interfaces with $\sigma_{13}=0.01$, where the solid line represents the contour levels of the volume fraction $c_{i}=0.5$. As time evolves, the significant deformations of phase $1$ and phase $2$ are generated, and the interface between the fluids $3$ and $4$ is curved due to the effect of the droplets on both sides (Fig. 13). Next, the films are formed between phase $1$ and phase $3$ and between phase $2$ and phase $3$, respectively (Fig. 13). And due to the disturbance caused by the movement of the droplets and bubbles, the interface of the two continuous phases produces a phenomenon similar to Rayleigh-Taylor (RT) instability. That is to say, the heavy and light fluids (phase $3$ and phase $4$) penetrate into each other and the penetration length increases with time, which then leads to the formation of the spike. As depicted in Fig. 13, at $t=12.25$, the interface between different fluids becomes more complicated. The bubble is squeezed into a strip, and two small drops of phase $3$ are trapped between the droplet (phase $1$) and the bubble (phase $3$).

We next look into the effect of the surface tension on the fluid interfaces. Here we set the surface tension between phase $1$ and phase $3$ to $\sigma_{13}=0.006$. Fig. 14 depicts the evolution of the fluid interfaces. One can observe that at $t=8.98$, the droplet (phase $1$) in Fig. 14 produces a larger deformation than in Fig. 13, and in Fig. 14, phase $3$ doesn’t break to form small droplets at $t=12.25$. To explain this, we introduce a dimensionless parameter, $Bo=\frac{\Delta\rho gD^{2}}{\sigma_{13}}$, which can be used to describe the ratio of droplet gravity to surface tension. In the former case, $Bo=10.8$, while in this case, $Bo=18.0$. As the $Bo$ number increases, the gravity of the droplet is more prominent, and the droplet is more susceptible to deformation (see Fig. 14). This deformation causes the droplet to be subjected to greater resistance during the falling process, resulting in a slower drop of the droplet. As shown in Fig. 15, when $\sigma_{13}=0.006$, the falling speed of the droplet $U$ is significantly smaller than the one when $\sigma_{13}=0.01$, where $U$ has been normalized by $U=\frac{U_{l}}{\sqrt{gD}}$, and $U_{l}$ is the velocity in the lattice unit. Eventually, a smaller drop velocity will cause the interface dynamics of this four-phase flow problem to be less complicated.

5.2 Dynamics of droplets and interfaces for five fluid phases

In this section, to demonstrate the ability to calculate the dynamic problem involving five fluid phases using the present model, we will consider a complex problem of droplet and interface motion involving five fluid phases. The initial setup of the physical problem is shown in Fig. 16. Similar to the previous example, the upper part of the domain is filled with phase $4$, and the lower part is filled with phase $5$. While in this simulation, two droplets (i.e., phase $1$ and phase $2$) of the same diameter are located inside phase $4$, and a bubble of the same size, referred to as phase $3$, is located in phase $4$. Then they start moving simultaneously under gravity (buoyancy). The simulation was carried out in a uniform mesh of $NY\times NX=320\times 200$ with periodic boundary conditions in the horizontal direction and halfway bounce-back boundary conditions in the vertical direction. The profiles of the volume fractions can be initialized by

[TABLE]

where $R$ is the radius of the circle drop for phase $1$, $2$ and $3$. $(x_{c_{1}},y_{c_{1}})$, $(x_{c_{2}},y_{c_{2}})$ and $(x_{c_{3}},y_{c_{3}})$ are the center of two droplets and the bubble, and $y_{0}=\frac{1}{2}NY$ is the initial surface position between phase $4$ and phase $5$. In this test, $R=30$, $(x_{c_{1}},y_{c_{1}})=(60,280)$, $(x_{c_{2}},y_{c_{2}})=(140,280)$, $(x_{c_{3}},y_{c_{3}})=(100,40)$. This complex five-phase problem has higher requirements for the selection of difference schemes. To further reduce the parasitic currents and guarantee the stability of present LB model, here we adopt the following isotropic mixed difference scheme [60]

[TABLE]

Although this scheme may not conserve the mass and momentum precisely due to the discretization errors [54], due to the better stability of this scheme, we will adopt this difference scheme for five-phase flows. Actually, the mass change can be controlled within $0.1\%$ after numerical tests, which can be ignored in the simulation. The other physical parameters in our simulations are set to be $\rho_{1}:\rho_{2}:\rho_{3}:\rho_{4}:\rho_{5}=6:4:1:2:3$, $\sigma_{ij}=0.01$, $(1\leq i\neq j\leq 5)$, $\tau_{g}=\tau_{1}=\tau_{2}=\tau_{3}=\tau_{4}=0.8$, $m_{1}=m_{2}=m_{3}=m_{4}=0.1$, $\eta=\sqrt{2}$, $\beta=\sqrt{0.06}$, $g=10^{-5}$. Fig. 17 depicts the evolution of the fluid interfaces for five fluid phases. It is found that after releasing at $t=0$, the two droplets of phase $1$ and phase $2$ fall freely in phase $4$, while the bubble of phase $3$ rise in phase $5$. Due to the higher density, phase $1$ falls faster than phase $2$, and the interface between phase $4$ and phase $5$ begins to deform due to the movement of the droplets [see Fig. 17]. Fig. 17 shows that the droplets and the bubble approach the surface between phase $4$ and $5$. As time goes by, phase $1$ continues to move downward, causing phase $4$ to penetrate into phase $5$ to form a spike [see Figs. 17 and 17]. At the same time, another droplet (phase $2$) and the bubble (phase $3$) are about to impact the interface between phase $4$ and phase $5$. At $t=10.21$, the droplet penetrates into phase $5$ from phase $4$. Meanwhile, phase $2$ and phase $3$ gradually become unstable, and the interface between them begins to incline to the right. Subsequently, phase $2$ continues to move upward due to the buoyancy, phase $3$ moves downward around phase $2$, and causes a small drop of phase $4$ to be trapped between phase $2$ and phase $5$.

Fig. 18 depicts the temporal sequence of snapshots of velocity fields at a few special moments. At $t=4.08$, the falling droplet of phase $1$ and the bubble (phase $3$) have induced a distinct velocity field. Subsequently, phase $1$ moves downward, penetrating into phase $5$ with a pair of vortices [see Fig. 18]. After that, the velocity field of phase $1$ gradually subsides [Figs. 18 and 18]. On the other hand, phase $2$ and phase $3$ produce a clockwise velocity field after contact at $t=12.25$, eventually causing phase $2$ and Phase $3$ to bypass each other and then move in the opposite direction.

We next investigate the effect of density ratio. In the present case, the density ratio is set to be $\rho_{1}:\rho_{2}:\rho_{3}:\rho_{4}:\rho_{5}=6:4:1:3:2$. As illustrated in Fig. 19, the current droplets and bubble (phases $1$, $2$ and $3$) move more slowly in contrast to the previous case, because the droplets and bubble are subject to smaller gravity (buoyancy) [Figs. 19-19]. Another significant difference is that phase $2$ and phase $3$ bypass each other in the opposite direction to the previous example [Figs. 19-19]. This phenomenon can be observed and explained more intuitively through the velocity field at time $t=12.25$ [Fig. 20]. At this moment, we find that the vortex between phase $2$ and phase $3$ is counterclockwise in the current situation, which is the opposite of the vortex direction in the previous example. Because the density of the upper fluid (phase $4$) is larger than that of the lower one (phase $5$), the interface between phase $4$ and phase $5$ undergoes a more distinct deformation [Figs. 19-19]. Finally, phase $1$ and phase $2$ fall to the bottom and phase $3$ rise to the top of the calculation domain.

6 Conclusions

In this study, a LB model based on the multicomponent phase field theory is proposed for $N$-phase flow systems ($N\geq 2$). In the proposed model, we reformulate the governing equations of the $N$-phase system proposed by Dong [4], and the new governing equations are also thermodynamically consistent. We point out that the velocity in this model is the volume-averaged mixture velocity and is divergence free. Besides, the mixing energy density coefficients in the free energy involve the interaction between different phases and can be determined by the pairwise surface tensions among the N fluids. In the present LB model, $(N-1)$ LB equations are employed to capture the interface, and another LB equation is used to solve the N-S equations. To reduce the calculation of the gradient term, a distribution function for the force term is carefully designed in the LB equation for flow field , and the governing equations can be recovered correctly through Chapmann-Enskog analysis.

The proposed model is first validated by two classical benchmark problems, including tests of static droplets and the spreading of a liquid lens, the simulation results show that the current LB model is accurate and efficient for simulating incompressible $N$-phase fluid flows. To further demonstrate the capability and performance of the LB model, we use two numerical simulations, including dynamics of droplet collision for four fluid phases and dynamics of droplets and interfaces for five fluid phases, to test the proposed model. The results show that the present model can be used to handle complex interactions among multiple types of fluid interfaces.

Acknowledgements

This work is financially supported by the National Natural Science Foundation of China (Grant Nos. 51576079, 51836003), and the National Key Research and Development Program of China (Grant No. 2017YFE0100100).

Appendix A Chapman-Enskog analysis of the present model

In this appendix, the Chapman-Enskog (C-E) expansion is performed to demonstrate the consistency of the present LB model with hydrodynamic equations [Eqs. (20)].

Before performing C-E expansion, the moment conditions are first given based on the expressions of the equilibrium and force distribution functions,

[TABLE]

where $\Delta_{\alpha\beta\gamma\theta}$ is given by $\Delta_{\alpha\beta\gamma\theta}=\delta_{\alpha\beta}\delta_{\gamma\theta}+\delta_{\alpha\gamma}\delta_{\beta\theta}+\delta_{\beta\gamma}\delta_{\alpha\theta}$, and Greek indices denote Cartesian spatial components.

In the C-E analysis, the time and space derivatives, as well as the force term can be expanded as

[TABLE]

where $\epsilon$ denotes a small expansion parameter. Using the Taylor expansion to Eq. (33), one have

[TABLE]

where $D_{k}=\partial_{t}+c_{k\alpha}\partial_{\alpha}$, and substituting Eq. (64d) into Eq. (65), the following multi-scale equations can be obtained,

[TABLE]

[TABLE]

where $D_{1k}=\partial_{t_{1}}+c_{k\alpha}\partial_{1\alpha}$.

Then, the substitution of Eq. (66b) into Eq. (66c) yields

[TABLE]

By summing Eq. (66b) and Eq. (66b)$\times c_{k\beta}$ over $k$, we can obtained the recovered equations at $\epsilon$ scale,

[TABLE]

Similarly, the recovered equations at $\epsilon^{2}$ scale can also be obtained from Eq. (67)

[TABLE]

where $\Lambda^{(1)}=\sum_{k}c_{k\alpha}c_{k\beta}g_{k}^{(1)}$ is the first-order momentum flux tensor.

Summing Eq. (32) and Eq. (32)$\times c_{k\alpha}$ over $k$, we can obtain the following relations

[TABLE]

which can be further recast as

[TABLE]

Substituting Eq. (71b) into Eqs. (68b) and (69b), one have

[TABLE]

[TABLE]

[TABLE]

Combining Eq. (72a) and Eq. (73a) at $\epsilon$ and $\epsilon^{2}$ scales with the help of Eqs. (63b) yields

[TABLE]

To derive the equation at $\epsilon^{2}$ scale, $\Lambda^{(1)}$ should be expressed as

[TABLE]

where Eq. (66b) has been used, and the term $\partial_{t_{1}}(\rho u_{\alpha}u_{\beta})$ can be written as

[TABLE]

where the Eqs. (72a) and Eq. (72b) are used, and the term of $O(Ma^{3})$ has been neglected.

Combining Eq. (72b) with Eq. (73b) at $\epsilon$ and $\epsilon^{2}$ scales, together with Eq. (75), Eq. (76) and Eq. (63b), we have

[TABLE]

with the relation of

[TABLE]

In the limit of a low Mach number, the dynamic pressure is assumed to be $\delta p\sim O(Ma^{2})$. As a result, Eq. (77) will reduce to

[TABLE]

Thus, Eqs. (74) and (79) clearly shows that the incompressible N-S equations can be exactly recovered from the present LB model.