The SC LBM has been implemented in two dimensions for the multiphase system. The distribution function is introduced for fluid components. And it satisfies the following lattice Boltzmann equation: f i ( r + c i Δ t , t + Δ t ) − f i ( r , t ) = − 1 τ [ f i ( r , t ) − f i e q ( r , t ) ] , (1) where f i ( r , t ) is the density distribution function which represents fractions of fluid molecules with a specific velocity c i at a specific position and time t when time and space are discretized by Δ t and Δ x ; τ is a relaxation time which is related to the kinematic viscosity as follows: μ = c s 2 ρ ( τ − 0.5 ) . (2)

The equilibrium distribution function f i e q ( r , t ) can be calculated as f i e q ( r , t ) = ω i   ρ [ 1 + 3 ( c i ⋅ u eq ) c s 2 + 9 ( c i ⋅ u eq ) 2 2 c s 4 − 3 u e q 2 2 c s 2 ] , (3) where c i is the discrete velocity. For the D2Q9 model, they are given by [ c 0 , c 1 , c 2 , c 3 , c 4 , c 5 , c 6 , c 7 , c 8 ] = c [   0   1   0 − 1             0           1         − 1       − 1           1 0   0   1       0         − 1           1                 1           − 1   − 1 ] . And where ω i is the weighting factor:   ω i = { 1 / 4 ,       i   =   0   1 / 18 , i = 1 ∼ 6     1 / 36 , i   =   7 ∼ 18 , (4) where c s = c / 3 , c = Δ x / Δ t is the ratio for lattice spacing Δ x , and the time step is Δ t . In Eq. 2, ρ is the density of fluid, which can be obtained from ρ = ∑ i f i . Macroscopic velocity u e q is given by u i e q = u + τ i F i ρ i , (5) where u is a velocity which is defined as u = ∑ i f i c i τ i ∑ i ρ i τ . (6)

In Eq. 5, F i = F i , c is the force acting on the fluid component, including fluid–fluid cohesion F i , c and excluding fluid–solid adhesion.

The cohesive force acting on the fluid component is defined as F i , c ( r , t ) = − G c ρ c ( r , t ) ∑ i ω i ρ c ' ( r + c i Δ t , t ) c i , (7) where ρ c and ρ c ' denote two different fluid components, and G c is a parameter that controls the strength of the cohesion force.

In this study, the initial density is set to be 1. The periodic boundary conditions are used for all boundaries in the 301 × 301 mesh. Initially, droplets are distributed uniformly in the domain. Small droplets grow up gradually by dispersion leading to a high degree of heterogeneity because larger droplets are able to minimize surface energy. Small droplets became less by coalesce due to dispersion and local flow behaviors. In the entire domain, the interface between two phases was decomposed into several high-density zones and low-density zones. Later, phase zones can further coalesce and grow up resulting minimized total lengths of the interface under surface tension. When the computations reached to the equilibrium state, two phases will separate completely. Small structures tend to have a spherical shape with a smooth surface before evaporating by diffusion (Reis and Phillips, 2007).

Various interactive forces have effects on speeds of phase separations. The density ratio and viscosity radio between two phases are represented as ρ ∗ and τ ∗ , respectively. As expected, higher interactive strengths lead to more clear and quick separations. A two-phase system with small interaction may fail to separate. In Figure 1, the system with the intermediate interactive strength is able to represent similar structures like those in the system with a low interactive strength earlier. At the final stage (Figure 1D), phase separation was shown as top and bottom bands. Initially, intermediate and final stages of phase separations with the highest interaction show the most clear interface boundary. Initial random structures have grown up into larger band-like structures and then coalesce. Under surface tension, the bands reshape into small droplets or bubbles for minimizing surface energy. At last, the domain was full of stably moving big droplets. If the simulation time is long enough, all droplets will merge into one. When we observe the separation behaviors at one specific time in a row, the interaction with a small interactive strength requires more time to evolve until the separation process is complete.

Phase separation happens in a system with an intermediate specified interactive force; see Figure 2. A clearer structure involving separated behaviors can be seen quickly. Distinctive band structures are formed with random swings but are hardly evolved into further structures within the simulation time.

When the interactive strength is risen to a high level with the magnitude of -100 shown in Figure 3, the exceptional separated structures are shown. Many smaller band structures are displayed in the earlier stage (Figure 3B). Subsequently, such “bands” are further combined (Figure 3C) and then reshaped into some round bubbles (Figure 3D).

Structure factor is used by the Fourier transform of spatial distributions of droplets or bubbles to present instabilities observed. Here, the first structure factor evolution will be showed as the function of simulation time. The structure factor shows heterogeneity changes due to phase separation (Zou et al., 1994). The structure factor is defined based on the fluctuation–dissipation theorem within integration of the imaginary part of the density response function. For a discrete system we are studying, it is given by S ( k , t ) = 1 N | ∑ r [ q ( r , t ) − q   ( t ) ¯ ] e i k ⋅ r | , (8) where N is the total number of grid points in the domain. q ( r , t ) = n 1 ( r , t ) − n 2 ( r , t ) ; k = 2 π L ( m i + n j ) , m . n = 1,2 , … , L , and L is the linear lattice size. q   ( t ) ¯ is the special average of q ( r , t ) at time t.

Structure factor is a quantitative measure of the structure along a prescribed direction. A value at a given wave number indicates the presence of heterogeneity in the periodic domain. Here, we present the first structure factors along the x and y directions. These structure factors are associated with heterogeneities developed to the size of the computational domain.

The onset of phase separation is defined as the time when the structure factor is above zero. In order to evaluate fluctuations from phase separation behaviors, graphic distributions are used corresponding to structure factors at selective times. Initially, random uniformed structures are specified in the domain, with a low value of structure factors; see Figure 4A. Then, small droplets start to merge with the neighboring ones, and therefore, relative bigger droplets or chunks were distributed with an increasing fluctuation presenting by higher structure factors, in Figure 4B. In Figure 4C, the structure factor of low interactive strengths is observed to be a peak. As separation continues, some droplets become larger than others, which results in uneven structures in the domain. These larger and more nonuniformed structures were verified by maximum of the mean structure factor value at the approximate 300th step. In Figure 4D, high- and low-density phases at the 1800th step were alternately distributed, and values of mean structure factors in Figure 4 remain.

Structure factors of the system with medium interactive strengths are shown in Figure 5. Similarly, the initial random uniformed distribution is shown in Figure 5A, which corresponds to low structure factors in Figure 5. And such low-level fluctuations are able to stay longer until approximately 200 steps as shown in Figure 5B due to increasing repulsive interaction compared to the low interactive strength. When large droplets further combine with neighboring same materials, some laminated structures appeared, and mean structure factors reach a peak at the approximate 800th step (Figure 5C). Although similar band structures are seen as those in the system with the low interactive strength, a clearer boundary between two phases is revealed.

When the interactive strength is increased to 100, quite clear phase separation is observed as shown in Figure 6. Highest values of mean structure factors among three systems are shown. From the beginning, distinct boundaries enclosing large droplets or chunks were taken shape as shown in Figures 6A and B. But, structure factors are low because degree of heterogeneity at these stages is relatively small. At the approximate 300th step (Figure 6C), structure factors are gradually fortified as a variety of sizes of droplets or chunks. Additionally, structure factors in this systems exhibit a smooth profile that is different from previous two systems with curved trends. This suggests that a two-phase system with strong repulsive interaction involves continuous fluctuations rather than abrupt changes.