Taylor instability, Droplet coalescence, Droplet deformation.

Multiphase and multicomponent flows are universal phenomena of scientific investigation and technological practice, from oil field exploitation to spray droplet formation, droplet transport, fuel cells, boiling processes, and microfluidics. Multiphase systems include immiscible and miscible fluids. In turn, such fluids involve moving and deformable interfaces whose dynamics are governed by the combined action of inertia, viscosity, gravity, and surface tension. Thus, numerically predicting the dynamics of such flows is very difficult, since, in addition to the strong non-linearity of the system dynamics, the interface itself may be strongly deformed, stretched, and split. At the same time, the phases involved may differ greatly in density and viscosity. Therefore, an accurate description of interfacial forces and dynamics is crucial for physically accurate modeling of multiphase flows [1–4].

There are several macroscopic numerical models, including the Volume of Fluid (VOF), Level-Set, and front-tracking methods, that are widely used to simulate immiscible multiphase flows. Though those models have proved to be quite successful, they usually include explicit interface reconstruction, tracking, or reinitialization procedures. These procedures may be computationally expensive and numerically challenging when interfaces deform strongly and/or change their topology. The Lattice Boltzmann Method (LBM) is a mesoscopic model of flow, in which macroscopic dynamics emerge from the evolution of particle distribution functions. Due to its local computational nature, ease of parallel implementation, and ability to model complex interfacial configurations without explicit reconstruction, LBM has become a popular tool for simulating multiphase and multicomponent flows [5].

There are several multiphase LBM approaches available, namely the color-gradient, free-energy, phase-field, and pseudopotential models [2], [3]. Among those models, the pseudopotential approach introduced by Shan and Chen [6], [7] has attracted considerable attention for its simplicity, local character, and ability to induce phase segregation via an effective interaction force. In the Shan-Chen approach, phase segregation and interface dynamics emerge from interactions among distribution functions without any explicit interface-tracking procedure. Thus, it becomes possible to model complex processes such as droplet deformation, splitting, and multiphase multicomponent processes; a pseudopotential approach was developed [8].

However, there are several drawbacks of the original Shan-Chen multiphase LBM worth noting. Firstly, thermodynamic inconsistency, independent control deficiency of the surface tension and the interfacial width, substantial spurious currents near curved interfaces, and reduced numerical stability for large density and viscosity ratios have been observed [9–12]. All the abovementioned drawbacks are associated with the discrete approximation of intermolecular forces and the forcing term. In Czelusniak et al. [13], the authors investigated the force-based formulation of the pseudopotential LBM. They showed how incorporating force into the LBE influences simulation accuracy and stability. This means that the force formulation should be regarded not only as a numerical technique but also as an integral component of the pseudopotential multiphase models that define their physical correctness.

Spurious current generation is one of the most critical numerical problems of pseudopotential multiphase LBM. Nonphysical velocities of fluid particles are especially pronounced near curved interfaces and adversely affect predictions of equilibrium and dynamic multiphase phenomena. As shown by Peng et al. [9], isotropy of the discrete interaction force is related to the generation of the mentioned numerical phenomenon. Moreover, the authors have shown that better discretization of the interaction force leads to a noticeable reduction in spurious currents. Qin et al. [10] have shown that even more accurate discrete-difference schemes for gradient evaluation can suppress spurious currents and improve predictions of coexistence densities. The abovementioned examples clearly show the significance of the proper interaction force discretization. There are also problems concerning the thermodynamic consistency of pseudopotential models. The coexistence densities determined by the lattice model do not match the corresponding equation of state at high density ratios. Peng et al. [11] proposed a customized equation-of-state approach, which guarantees strict thermodynamic consistency and independent control of surface tension. The authors' results show that a proper modification of the equation of state can significantly improve the agreement between the numerical and thermodynamic coexistence densities. Recently, Li et al. [12] have systematically analyzed the thermodynamic properties of the pseudopotential LBM for planar and curved multiphase interfaces, and developed an improved formulation.

The fluid-fluid interactions, apart from force discretization and thermodynamic consistency, also depend heavily on the choice of interaction scheme and fluid properties. The studies by Li and Liu [14] on fluid-fluid interactions in pseudopotential lattice Boltzmann models have revealed that changes in the scheme and in fluid properties can significantly affect phase separation, interfacial behavior, and the characteristics of the resulting multiphase system. Such information is highly relevant to the current work, as the original Shan-Chen model has only a small number of tunable parameters, and interaction strength is a key parameter that affects phase separation, interface thickness, and interfacial tension. In this regard, the systematic characterization of the influence of parameters is required before applying the formulation to complex multiphase flows.

The use of multicomponent flow systems introduces additional complexity, as it requires an accurate description of interactions among various fluid components while maintaining numerical stability over a wide range of density and viscosity ratios. Zhao et al. [15] have performed a Chapman-Enskog analysis of the modified pseudopotential model of multicomponent flow and found that the original formulation may be limited for large viscosity ratios. In their study, the modified formulation has expanded the range of applicability of the viscosity ratio, especially when combined with the Multi-Relaxation-Time (MRT) collision operator. Yang et al. [16] have proposed a pseudopotential LBM formulation for multicomponent and multiphase flow at high density ratios, using the improved forcing scheme, the MRT collision operator, and the equation-of-state scaling coefficient.

The use of an MRT collision operator in multiphase simulations is increasingly desirable due to its flexibility in independently controlling kinetic moment relaxations and its greater numerical stability than the SRT model. These advantages are increasingly important for multiphase flow simulations with large density differences, interfacial deformations, and high Reynolds or Weber numbers. Hosseini et al. [17] verified the robustness of an LBM multiphase solver for multiphase flow conditions with relatively high Weber and Reynolds numbers. This work further brought to the fore ongoing developments toward the optimization of collision and forcing formulations for complex interfacial flow applications. However, the increased complexity of such advanced formulations points to the need for systematic evaluation of the fundamental Shan-Chen model.

The use of pseudopotential LBM has shown promising results in addressing various multicomponent transport phenomena encountered in real applications. For example, Moslemi et al. [18] used the multicomponent, multiphase pseudopotential LBM with an MRT collision scheme to investigate the transport of liquid water through compressed porous transport layers in PEM fuel cells. A forcing-treatment mechanism was incorporated into the computational model and verified using permeability measurements. Droplet-pressure analyses were performed, thereby demonstrating the framework's ability to simulate multicomponent, multiphase transport under real engineering conditions with large density and viscosity ratios. Javaherdeh et al. [19] also conducted a study on the transport of liquid water through compressed gas diffusion and microporous layers using the Shan–Chen multicomponent multiphase model coupled with MRT collisions. They studied the effects of compression, liquid water injection, capillary number, and contact angle. Additionally, LBM was employed by Ashorynejad et al. [20] to simulate liquid water transport in a wavy gas channel, and its capability to capture droplet-gas interaction processes was demonstrated.

Pseudopotential LBM has achieved encouraging results in solving multicomponent transport phenomena encountered in real applications. For example, the multicomponent, multiphase pseudopotential LBM with an MRT collision scheme was used by Moslemi et al. [18] to study the transport of liquid water through compressed porous transport layers in PEM fuel cells. Because of the rotational model, a significantly different mechanism was introduced and validated through permeability and droplet-pressure analyses, demonstrating the framework's ability to simulate multicomponent, multiphase transport under real engineering conditions with large density and viscosity ratios. Javaherdeh et al. [19] studied the transport of liquid water through compressed gas diffusion and microporous layers using the Shan–Chen multicomponent multiphase model combined with MRT collisions. They investigated the effects of compression, liquid water injection, capillary number, and contact angle. Furthermore, Ashorynejad et al. [20] employed LBM to simulate liquid water transport in a wavy gas channel, demonstrating its ability to capture droplet-gas interaction processes.

Furthermore, the validity of pseudopotential models has been demonstrated in simulations of confined multicomponent flows. Li and Liu [14] developed a pseudopotential LBM method for simulating immiscible multicomponent flows in microchannels. The method was demonstrated to predict interfacial behavior under microscale confinement. These simulations are highly relevant, as confinement can strongly influence the droplet shape, interfacial transport phenomena, and phase distribution within the domain. For example, Phelan et al. [21] investigated the droplet dynamics at supercritical capillary numbers using the LBM.

Alternative multiphase LBM formulations have also been developed to describe the dynamics of compound droplets. Chen et al. [22] developed a ternary phase-field reduced multiphase LBM and investigated the dynamics of compound droplets on a solid surface under shear flow. They showed that phase-field-based lattice Boltzmann formulations can effectively capture complex interfacial configurations and the deformation of compound droplets. These approaches are conceptually distinct from the pseudopotential formulation, but they provide useful comparisons regarding the difficulties of accurately resolving multiple interfaces and complex droplet dynamics. Therefore, benchmark problems of droplet deformation and breakup are used to stringently test the interfacial accuracy and stability properties of the Shan-Chen model.

Another widely adopted problem for validating multiphase numerical methods concerns the investigation of Rayleigh–Taylor instability, given its relevance to gravity-induced instabilities, interface deformations, and the complex dynamics of interfacial formation. Tavares et al. [23] proposed a multicomponent lattice Boltzmann model using the Shan–Chen pseudopotential theory to simulate immiscible Rayleigh–Taylor turbulent flows. They demonstrated the ability of this numerical approach to predict complex, evolving interfacial dynamics. Successful predictions of these instability issues are significant because they simultaneously account for the effects of gravity forces, interfacial tension, density differences, and nonlinear deformation behaviors. Although significant improvements have been achieved owing to innovations in forcing treatments, equations of state, MRT collision models, and multiphase approaches, the proper determination of the physical parameters in the base Shan–Chen model remains crucial. In particular, the intensity of intermolecular interactions significantly influences phase separation and interfacial characteristics. In contrast, the surface tension, interfacial thickness, and spurious velocities are usually strongly interrelated through model parameters. Hence, numerical investigations into the correlations among interaction intensity, surface tension, interface width, and spurious velocity are essential prior to employing the Shan–Chen model in complex flow regimes.

Evaluation of surface tension and validation of interfacial force formulations based on pressure jumps across curved interfaces have become de facto standard procedures [9–13], [23]. This is understandable, since determining the interface shape and stability requires an accurate prediction of the magnitude of the interfacial force. Several recent studies have highlighted the significance of this issue for force formulation, thermodynamic consistency, and interfacial behavior. Although there are many sophisticated multiphase LBMs for particular problems such as spurious currents, thermodynamically consistent treatments, high-density-ratio multiphase flows, or droplet dynamics, a systematic assessment of the basic Shan-Chen model with different benchmark problems still provides important insights into its fundamental properties. The adequacy of the simple Shan-Chen multicomponent model for practical purposes can be assessed by studying a wide variety of properties, e.g., equilibrium interfacial behavior, mechanisms of instability, droplet interaction phenomena, oscillatory flow in confinement, shear-driven deformation.

In view of the above background, the present study aims to conduct a comprehensive validation analysis to predict immiscible multiphase flow phenomena using the Shan-Chen multiphase model. First, we study how the interface thickness, spurious velocities, and interfacial properties of a stationary droplet depend on the strength of intermolecular interactions. Then, the surface tension will be measured by the Laplace test, and the interfacial force will be checked. Then, the accuracy and reliability of the Shan-Chen model will be tested using benchmark cases such as Rayleigh-Taylor instability, droplet collision and coalescence, oscillating droplets in a U-tube, and droplet deformation under shear flow. In particular, the performance of the proposed approach will be evaluated for predicting different regimes of droplet deformation across combinations of the Reynolds and capillary numbers in shear flows. All results obtained will be critically compared with existing theoretical, experimental, and numerical data in the literature to validate and assess the multiphase modeling scheme proposed herein quantitatively.

In this work, we use a multicomponent Shan-Chen pseudopotential Lattice Boltzmann Method (SC-LBM) with the SRT collision model and D2Q9 lattice configuration to simulate immiscible multiphase flows. In this setup, interactions between components are handled using the Shan-Chen pseudopotential approach, and the body forces follow the methodology of Guo et al. [24]. The formulation allows the inclusion of interfacial interaction forces in the LB equation without the need to identify the interface explicitly. The evolution equation for the distribution function of each fluid component is written as

k k 1 k k,eq k

fi (x + eiΔt, t + Δt)fi (x, t) − τ [fi (x, t) − fi (x, t)] + Δt, Fi , (1)

where fik(x, t) is the distribution function of component (k) in the (i)-th lattice direction, (ei) is the corresponding discrete lattice velocity, (τk) is the relaxation time, and (Fik) is the discrete forcing term associated with the total force acting on component (k). The relaxation time is related to the kinematic viscosity of each component according to

νk = cs 2 (τk − 1) Δt, (2)

where νk is the kinematic viscosity and (cs) is the lattice speed of sound. The equilibrium distribution function is obtained from the second-order expansion of the Maxwell–Boltzmann distribution:

k,eq ei ⋅ u (ei ⋅ u)2 u ⋅ u

fi = wiρk [1 + 2 + 4 2 ], (3)

cs 2cs 2cs where (wi) is the weighting coefficient associated with the (i)-th lattice direction, ρkis the density of component (k), and you denotes the common macroscopic velocity of the multicomponent system. The density of each component is calculated from the zeroth-order moment of its distribution function:

ρk = ∑ fik. (4)

Accordingly, the total density of the mixture is obtained from

ρ = ∑ ρk (5)

The macroscopic velocity is calculated by considering the half-time-step correction associated with the forcing term, as proposed by Guo et al. [24]:

k Δt ∑k ∑i fi ei + F

u = 2 , (6)

where (F) represents the total force acting on the multicomponent system and is obtained from

F = ∑ Fk. (7) For the immiscible multicomponent system considered in the present study, the force acting on each component is generated exclusively by the intermolecular interaction between different fluid components. Thus, no gravitational or other external body force is considered. The Shan–Chen interaction force acting on component (k) is expressed as

Fk SC(x) = Gkk′ψk(x) ∑ wi ψk′(x + eiΔt)ei k ≠ k′, (8)

where (Gkk') is the interaction strength between components (k) and (k'), and ψk and (ψk′) denote their corresponding pseudopotentialpseudopotentialpseudopotentialpseudopotential functions. The interaction parameter (Gkk′) controls the strength of the intermolecular interaction and consequently affects phase separation and the resulting interfacial properties. The pseudopotential function is defined as

ρk

ψk(ρk) = ρ0 [1 − exp (− ρ )], (9)

where ρ0 is a reference density. This nonlinear pseudopotential allows phase separation to emerge naturally from the intermolecular interactions without explicitly tracking the fluid interface. To incorporate the Shan– Chen interaction force into the lattice Boltzmann equation, the forcing term is formulated according to the Guo forcing scheme:

k ei − u (ei ⋅ u)ei SC

Fi = wi [ 2 + 4 ] ⋅ Fk . (10)

cs cs The above formulation provides a second-order accurate representation of the body-force contribution and improves the consistency between the forcing term and the macroscopic Navier–Stokes equations. For the D2Q9 lattice employed in this study, the discrete velocity set is defined as

e0 = (0,0),

e1−4 = c[(1,0), (0,1), (−1,0), (0, −1)], (11)

e5−8 = c[(1,1), (−1,1), (−1, −1), (1, −1)], where (c = Δx/Δt) is the lattice velocity. The corresponding lattice weighting factors are given by

4 1 1

w0 = , w1−4 = , w5−8 = . (12)

9 9 36 For the D2Q9 lattice, the lattice speed of sound is

cs 2 = c2. (13)

Therefore, the SRT relaxation time can alternatively be expressed as

νk 1

τk = 2 + . (14)

cs Δt 2 The mathematical framework of the multicomponent Shan–Chen lattice Boltzmann model used in the present study is given by Eqs. (1)–(14). The model is based on the SRT collision operator for computational simplicity, the D2Q9 lattice for two-dimensional simulations, the Shan-Chen intermolecular interaction force to trigger phase separation, and the Guo forcing scheme to consistently account for the interaction force in the lattice Boltzmann evolution equation.

A set of equilibrium and dynamic benchmarks validates the multicomponent Shan-Chen pseudopotential LBM. First, we consider the system's sensitivity to the range of the intermolecular interaction parameter (G). This parameter is important in the design of the Shan-Chen formulation and determines the phase separation behavior and interfacial properties. To this end, the influence of G on the interfacial width and density profile is investigated by simulating a stationary droplet. The next step is to determine the surface tension using the Laplace test. Subsequently, the prediction of more complex dynamic multiphase behaviors, such as Rayleigh-Taylor instability, droplet coalescence, oscillating flow in a U-bend geometry, and droplet deformation in shear flow, will be considered.

Fig. 2 shows the density distribution for the stationary droplet case. The density distribution remains nearly constant in the bulk phase region and changes abruptly across the interface. Increasing G from 1.6 to 2.5 results in a sharper density gradient and a reduced interface thickness. In contrast, further increasing G to 3.5 yields only marginal improvements. Thus, G=2.5 represents a reasonable trade-off and serves as the benchmark for future simulations.

This study has numerically evaluated the multicomponent Shan-Chen pseudopotential lattice Boltzmann model for simulating immiscible multiphase flows. Specifically, the study aimed to evaluate whether the original formulation can accurately predict equilibrium interfacial properties and complex transient behaviors in multiphase flows. An SRT collision operator based on the D2Q9 lattice structure and forced by the Guo formulation was adopted, and various benchmarks were selected to reflect increasing degrees of interfacial complexity .

The effect of the intermolecular interaction strength (G) was initially investigated via a stationary droplet problem. It was found that (G) played a crucial role in affecting the separation process and interfacial properties. Higher values of (G) yielded a sharper profile and thinner interfaces. Nonetheless, it should be noted that further enhancement of interfacial sharpness would be relatively minor once certain limits are exceeded, suggesting that excessively high interaction strength does not ensure proportional improvements in interface resolution. Considering the tested cases, the value of (G=2.5) offered a relatively appropriate balance of separation capability and sharpness.

Surface tension was calculated based on the Laplace criterion. The variation in the pressure difference between stationary circular droplets with respect to the inverse of the droplet radii showed a nearly linear trend, validating the Laplace relation over the studied range. It was observed that the surface tension depends on the strength of intermolecular interaction, wherein larger values of (G) yield higher surface tensions. Such observations emphasize the importance of calibrating the interaction parameter before using the Shan–Chen model to simulate dynamic multiphase problems in which capillary forces dominate.

The proposed model was assessed using multiple benchmark tests of dynamic behavior. In the case of Rayleigh–Taylor instability, the proposed approach was shown to predict evolving interface deformation, finger growth, nonlinear instability, and significant topological changes. Droplet coalescence accurately reproduced the sequence: droplet approach, formation of a liquid bridge between them, complete merger, and subsequent relaxation driven by capillary forces. Oscillations in a U-bend were used to verify the current model's ability to simulate the combined effects of inertia, restoration, and viscous dissipation, which attenuate oscillatory motion.

Droplet deformation under shear flow has been studied as a function of the capillary and Reynolds numbers. As expected, an increasing trend in droplet deformation with increasing capillary number has been observed, indicating greater relevance of viscous stress relative to surface-tension effects. Results from the current simulation have shown reasonable correlation with experimental data available in the literature.

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

The data that support the findings of this study are available within the article. No additional data are associated with this publication.

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.