*2.3. Forcing Scheme*

The adopted forcing scheme is the one by Guo et al. [34], as it gave the most stable simulations, lowest interfacial spurious currents and, hence, the highest stable interaction parameter when compared to the simple method by Buick and Greated [35] or the velocity shifting method, which was adopted in the original SC model [10,11]. The scheme by Guo et al. [34] is originally evaluated while considering the discrete lattice effects in order to recover the macroscopic continuity and momentum with external force. It was also reported by Yu and Fan [36] that the model by Guo can produce a *τ*-independent (viscosity independent) surface tension.

The evaluation of the source term can be written in the terms of the SC interaction force as follows:

$$S\_i^{(\sigma)} = \left(1 - \frac{1}{2\ \tau^{(\sigma)}}\right) R\_i^{(\sigma)}; \quad R\_i^{(\sigma)} = \frac{w\_i}{c\_s^2} \left(F\_{\text{SC}\_a}^{\text{C}\_{(\sigma)}} c\_{i\mathfrak{a}} + \frac{e\_{i\mathfrak{a}} c\_{i\mathfrak{f}} - c\_s^2 \delta\_{\mathfrak{a}\mathfrak{f}}}{c\_s^2} F\_{\text{SC}\_a}^{\text{C}\_{(\sigma)}} u' \beta\right) \tag{13}$$

where *α*, *β* = 1, 2 for the 2D case, while the corrected macroscopic velocity of the mixture *u*ˆ**'** shall written as a function of the interaction (or any external forces) as follows:

$$\hat{\mathfrak{u}}' = \frac{1}{\rho} \sum\_{\sigma} \left( \sum\_{i=0}^{8} f\_i^{(\sigma)} \mathfrak{k}\_i + \frac{\mathfrak{k}\_{\text{SC}}^{(\sigma)}}{2} \right); \quad \rho = \sum\_{\sigma} \rho^{(\sigma)} \tag{14}$$

Any external gravity or wall forces can be also added directly to the SC interaction forces.

#### *2.4. Two Relaxation Time*

Two phase flows always suffer from instabilities, especially when reaching the limiting cases of interaction parameters which are essential in order to tune density ratio and the surface tension together. Additionally, high viscosity ratio is a very high source of instability. This will always lead to the two relaxation time (TRT), which gives an extra degree of freedom for the simulations, yet in a simple way. This is available by splitting the distribution functions to symmetric *fi* (*σ*) <sup>+</sup> = *fi* (*σ*) <sup>+</sup> *<sup>f</sup>*−*<sup>i</sup>* (*σ*) /2 and asymmetric parts *fi* (*σ*) − = *fi* (*σ*) <sup>−</sup> *<sup>f</sup>*−*<sup>i</sup>* (*σ*) /2, where −*i* stands for the opposite lattice velocity direction. Each of the new distribution functions is relaxed with two separated relaxation times *τ*(*σ*)<sup>+</sup> and *τ*(*σ*)<sup>−</sup> , and controlled through the magic parameter Λ, as shown in Equation (15). It is reported that the value of Λ = 1/4 provides the highest stability [14].

$$
\Lambda = \left(\tau^{(\sigma)^+} - 0.5\right) \left(\tau^{(\sigma)^-} - 0.5\right); \quad \tau^{(\sigma)^+} = 3\mathfrak{v}\_{LII}{}^{(\sigma)} + 0.5\tag{15}
$$

while the source term according to Guo et al. [34] shall also be modified as follows [37,38]:

$$\mathcal{S}\_{i}\,^{(\sigma)}\_{TRT} = \left(1 - \frac{1}{2\tau^{(\sigma)}}\right) \mathcal{R}\_{i}\,^{(\sigma)} + \left(1 - \frac{1}{2\tau^{(\sigma)}}\right) \mathcal{R}\_{i}\,^{(\sigma)} \tag{16}$$

where *Ri* (*σ*) <sup>+</sup> = *Ri* (*σ*) <sup>+</sup> *<sup>R</sup>*−*<sup>i</sup>* (*σ*) /2 and *Ri* (*σ*) − = *Ri* (*σ*) <sup>−</sup> *<sup>R</sup>*−*<sup>i</sup>* (*σ*) /2 following same idea of splitting the distribution functions. Apart from the collision and forcing steps, the algorithm is exactly the same as the SRT.

#### *2.5. Mid-Range Interaction Model*

A positive value of disjoining pressure Π is a major macroscopic property of stable thin films, due to the existence of repulsive forces between protein or surfactant molecules [3]. The disjoining pressure can be related to the film tension *γ<sup>f</sup>* as follows:

$$\gamma^f = 2\gamma - \int\_{\infty}^{l} \Pi dl + \Pi l \tag{17}$$

where *γ* and *l* are the surface tension and the film thickness, respectively.

The D2Q25 model is used in order to reach positive values of disjoining pressure, hence stable films and foam [20,21]. The neighboring nodes, far till <sup>√</sup><sup>8</sup> from the central one, are considered as represented in Figure 1 (both red and blue colors). This will yield to a new interaction parameter responsible for the repulsive forces between the same component, which can be called *<sup>G</sup>*,<sup>11</sup> <sup>=</sup> *<sup>G</sup>*,<sup>22</sup> <sup>&</sup>gt; 0. The extra SC force shall be added to Equation (10), and can be written as:

$$\hat{\vec{F}}\_{\text{SC}}^{(\sigma)}(\hat{\mathfrak{x}}) = -\,\,\Psi^{(\sigma)}(\hat{\mathfrak{x}}) \tilde{G}\_{\sigma\sigma} \sum\_{\substack{\bar{i}=0\\\bar{i}=0}}^{24} \bar{w}\_{\frac{\bar{i}}{\bar{i}}} \Psi^{(\sigma)}(\hat{\mathfrak{x}} + c\_{\bar{i}} \Delta t) \,\,\mathfrak{k}\_{\bar{i}} \tag{18}$$

where *<sup>e</sup>*ˆ,*<sup>i</sup>* and *<sup>w</sup>*, ,*<sup>i</sup>* represent the mid-range lattice unit vectors and the weighting parameter as shown in Table 2.

**Table 2.** Weighting parameter *<sup>w</sup>*, ,*<sup>i</sup>* values for the D2Q25 lattice.


The full non-ideal static pressure (EOS) shall then be modified using the scope of Equation (9) as follows:

$$p = c\_s^2 \sum\_{\sigma} \left( \rho^{(\sigma)} + \frac{1}{2} \left( G\_{\sigma \sigma} + \tilde{G}\_{\sigma \sigma} \right) \Psi^{(\sigma)^2} \right) + \frac{c\_s^2}{2} \sum\_{\sigma, \bar{\sigma}} G\_{\sigma \bar{\sigma}} \, \Psi^{(\sigma)} \, \Psi^{(\bar{\sigma})} \tag{19}$$
