*2.1. Lattice Boltzmann Method*

The lattice Boltzmann method is based on the well-known Boltzmann Equation, where particles' motion and interaction in space and time are represented by the particle distribution function *f* , as shown in Equation (1) [30].

$$\frac{\partial f}{\partial t} + \mathfrak{d} \cdot \nabla f = \mathcal{Q} \tag{1}$$

where *t*, *v*ˆ and *Q* are time, particle microscopic velocity and the collision operator. *Q* is originally found in a complex integral form which can be approximated using the Bhatnagar–Gross–Krook (BGK) relaxation, namely the single relaxation time (SRT) [31]. It represents the particles' relaxation from non-equilibrium to the local equilibrium state during collision, through the time *τ* as shown in Equation (2) [30], where the equilibrium distribution function *f equ* can be formulated in the scope of the kinetic theory of gases using the Maxwell–Boltzmann distribution.

$$Q = -\frac{1}{\pi} [f - f^{equ}] \tag{2}$$

Discretizing the particle distribution function into a subspace of discrete lattice velocities *i*, mimicking the degrees of freedom of where the particles can propagate in space, the distribution function can be shown as *fi*(*x*ˆ, *t*). This can be explained as the probability of finding a particle in position vector *x*ˆ at time *t* with a microscopic speed in the direction *i*, using the explicit time integration, which when based on the trapezoidal rule [32] maintains the second order accuracy of the LBM in time [14]. The lattice Boltzmann equation (LBE) finally reads as:

$$f\_i(\mathbf{\hat{x}} + \mathbf{\hat{e}}\_i \Delta t, t + \Delta t) = f\_i(\mathbf{\hat{x}}, t) - \frac{1}{\pi} (f\_i(\mathbf{\hat{x}}, t) - f\_i^{c\eta u}(\mathbf{\hat{x}}, t)) \tag{3}$$

The above equation carries the core procedure of the LBM, showing the propagation of particles from one lattice position to the other and the collision operation, i.e., relaxation [33]. The equilibrium distribution function is based on the Hermite series expansion of the Maxwell–Boltzmann distribution with the truncation up to the second order macroscopic velocity, which is valid for low lattice Mach numbers and can be written as:

$$f\_i^{\text{equ}} = w\_i \rho \left[ 1 + 3 \frac{\text{\hat{e}} \cdot \text{\hat{n}}}{c^2} + \frac{9}{2} \frac{\left( \text{\hat{e}} \cdot \text{\hat{n}} \right)^2}{c^4} - \frac{3}{2} \frac{\text{\hat{n}} \cdot \text{\hat{n}}}{c^2} \right] \tag{4}$$

where *e*ˆ*i*, *wi*, *c*, *ρ* and *u*ˆ are the lattice unit vectors, weighting parameter, lattice speed (lattice spacing/lattice time step size Δ*x*/Δ*t* = 1), macroscopic density and velocity, respectively.

The macroscopic density and velocity can be evaluated from the 0th and 1st moment of the distribution function, respectively, as follows:

$$\rho = \sum\_{i=0}^{8} f\_i \; ; \quad \mathfrak{h} = \frac{1}{\rho} \sum\_{i=0}^{8} f\_i \mathfrak{e}\_i \tag{5}$$

The common D2Q9 is used here representing a 2D case and 9 lattice directions. The weighting parameters for the D2Q9 according to the short-range interaction from Figure 1 (presented in red) are listed in Table 1. Although LBM in accordance with the SC model (Section 2.2) is easily represented in 3D, the interest in this work is limited to the 2D case. Taking into consideration that the 2D mid-range interaction model (Section 2.5) requires the information from 24 neighboring nodes, the 3D case would require 92 in order to maintain the same 8th order of isotropy [21], which means much higher computational cost.

**Table 1.** Weighting parameter *wi* values for the D2Q9 lattice.


The LBM relaxation time can be related to the fluid kinematic viscosity through the Chapman–Enskog expansion, and comparing it with the incompressible Navier–Stokes equation leads to the following equation:

$$
\pi = \mathfrak{w}\_{LII} + 0.5 \tag{6}
$$

where *vLU* stands for the kinematic viscosity in lattice unis (LU). The Mach number *M*, lattice speed of sound *cs* and static pressure *p* based on the isothermal ideal equation of state can be presented, respectively, as follows:

$$M = \frac{||\hat{u}||}{c\_s}, \quad c\_s = \frac{c}{\sqrt{3}}, \quad p = \rho c\_s^{-2} \tag{7}$$

**Figure 1.** D2Q9 Lattice for short range nodes (red only). D2Q25 Lattice for short- and mid-range nodes (red and blue).

#### *2.2. Shan-Chen Multiphase Multicomponent Model*

The main benefit of the Shan-Chen [10,11] multiphase multicomponent model is that it does not dictate tracking any interface, since the separation automatically occurs due to the repulsive and attraction forces between components and phases, respectively. The SC is a pseudopotential method based on the potential function *ψ* as shown in Equation (8), which leads to the non-ideal equation of state as shown in Equation (9). The potential function can be written in two cases, and case *A* is reported to be more stable than case *B* [19] although case *B* will be used in a specific case as explained in Section 2.5.

$$\psi^{(\sigma)} = \begin{cases} \rho\_o \left[ 1 - e^{\frac{-\rho^{(\sigma)}}{f^{\rho}}} \right] & \text{case } A \\\rho^{(\sigma)} & \text{case } B \end{cases} \tag{8}$$

$$p = c\_s^2 \sum\_{\sigma} \left( \rho^{(\sigma)} + \frac{1}{2} \mathbf{G}\_{\sigma \sigma} \cdot \boldsymbol{\Psi}^{(\sigma)^2} \right) + \frac{c\_s^2}{2} \sum\_{\sigma, \{\overline{\sigma}\}} \mathbf{G}\_{\sigma \overline{\sigma}} \cdot \boldsymbol{\Psi}^{(\sigma)} \cdot \boldsymbol{\Psi}^{(\overline{\sigma})} \tag{9}$$

where *<sup>ρ</sup>o*, *<sup>σ</sup>*, ,*<sup>σ</sup>* and *<sup>G</sup>* represent the reference density which is taken as unity, the respected component, opposite components indices and the interaction parameter between phases and components, respectively. Here, both the interaction between phases and components are included. Only two components are considered in this work, 1 and 2. Setting the value of *G*<sup>11</sup> = 0, *G*<sup>12</sup> = *G*<sup>21</sup> > 0, and *G*<sup>22</sup> < 0 will lead to the presence of the two components with the coexistence of the component 2 in both phases. Tuning the two non-zero parameters will lead to controlling both the density ratio and surface tension separately.

The local SC interaction force for each component with respect to the neighboring nodes will then be presented as:

$$\mathcal{F}\_{\mathbf{SC}}^{(\sigma)}(\hat{\mathbf{x}}) = -\boldsymbol{\psi}^{(\sigma)}(\hat{\mathbf{x}}) \left( \mathbf{G}\_{\sigma\overline{\sigma}} \sum\_{i=0}^{8} w\_{i} \,\boldsymbol{\psi}^{(\overline{\sigma})}(\hat{\mathbf{x}} + c\_{i}\Delta t) \,\hat{\mathbf{z}}\_{i} \, \right. \\ \left. + \,\left. \mathbf{G}\_{\sigma\sigma} \sum\_{i=0}^{8} w\_{i} \,\boldsymbol{\psi}^{(\sigma)}(\hat{\mathbf{x}} + c\_{i}\Delta t) \,\hat{\mathbf{z}}\_{i} \right) \right. \\ \left. \left. \left( \mathbf{1} \right) \left( \mathbf{G}\_{\sigma\sigma} \sum\_{i=0}^{8} w\_{i} \,\boldsymbol{\psi}^{(\sigma)}(\hat{\mathbf{x}} + c\_{i}\Delta t) \,\hat{\mathbf{z}}\_{i} \right) \right) \right. \\ \left. + \left. \mathbf{G}\_{\sigma\sigma} \sum\_{i=0}^{8} w\_{i} \,\boldsymbol{\psi}^{(\sigma)}(\hat{\mathbf{x}} + c\_{i}\Delta t) \,\hat{\mathbf{z}}\_{i} \right) \right) \right)$$

For the short range interactions, the potential function *ψ*(*x* + *ci*Δ*t*) requires the density information from the first level of neighbors, which are 8 for the case of D2Q9.

Two LBE are then required for both components with two sets of distribution functions, which can be written as:

$$f\_i^{\ \ \ \ \left(\sigma\right)}\left(\hat{\mathfrak{x}} + \hat{\mathfrak{e}}\_i \Delta t, \ t + \Delta t\right) = f^{\ \left(\sigma\right)}\left(\hat{\mathfrak{x}}, t\right) - \frac{1}{\tau^{\left(\sigma\right)}} \left(f\_i^{\ \left(\sigma\right)}\left(\hat{\mathfrak{x}}, t\right) - f\_i^{\alpha \mu^{\left(\sigma\right)}}\left(\hat{\mathfrak{x}}, t\right)\right) + S\_i^{\left(\sigma\right)}\left(\hat{\mathfrak{x}}, t\right) \tag{11}$$

where *Si* (*σ*)(*x*ˆ, *t*) is a source term that incorporates the interaction and any other external forces based on the selected forcing scheme shown in Section 2.3, while *τ*(*σ*) is responsible for the kinematic viscosity of each component. The equilibrium distribution function for each component can be written as:

$$f\_i^{c\eta\mu^{(\sigma)}} = w\_i \rho^{(\sigma)} \left[ 1 + 3\frac{\mathfrak{e}\_i \cdot \hat{\mathfrak{u}}'}{c^2} + \frac{9}{2} \frac{\left(\mathfrak{e}\_i \cdot \hat{\mathfrak{u}}'\right)^2}{c^4} - \frac{3}{2} \frac{\hat{\mathfrak{u}}' \cdot \hat{\mathfrak{u}}'}{c^2} \right] \tag{12}$$

where *<sup>ρ</sup>*(*σ*) and *<sup>u</sup>*ˆ**'** are the density of each component evaluated as *<sup>ρ</sup>*(*σ*) <sup>=</sup> <sup>8</sup> ∑ *i*=0 *fi* (*σ*) and the corrected macroscopic velocity of the mixture, which is shown in Equation in 14 in Section 2.3.
