Multipseudopotential interaction : A consistent study of cubic equations of state in lattice Boltzmann models

A method is developed to analytically and consistently implement cubic equations of state into the recently proposed multipseudopotential interaction (MPI) scheme in the class of two-phase lattice Boltzmann (LB) models [S. Khajepor, J. Wen, and B. Chen, Phys. Rev. E 91, 023301 (2015)]. An MPI forcing term is applied to reduce the constraints on the mathematical shape of the thermodynamically consistent pseudopotentials; this allows the parameters of the MPI forces to be determined analytically without the need of curve fitting or trial and error methods. Attraction and repulsion parts of equations of state (EOSs), representing underlying molecular interactions, are modeled by individual pseudopotentials. Four EOSs, van der Waals, Carnahan-Starling, PengRobinson, and Soave-Redlich-Kwong, are investigated and the results show that the developed MPI-LB system can satisfactorily recover the thermodynamic states of interest. The phase interface is predicted analytically and controlled via EOS parameters independently and its effect on the vapor-liquid equilibrium system is studied. The scheme is highly stable to very high density ratios and the accuracy of the results can be enhanced by increasing the interface resolution. The MPI drop is evaluated with regard to surface tension, spurious velocities, isotropy, dynamic behavior, and the stability dependence on the relaxation time.


I. INTRODUCTION
The lattice Boltzmann method has emerged from the kinetic-theory approach and successfully found many applications in various aspects of computational fluid dynamics (CFD) such as turbulence [1], flows in porous media [2,3], blood rheology [4,5], and biopolymer translocation [6,7] while having the capability of being combined with the traditional CFD methods [8,9].Filling the gap between continuum and molecular levels, the LB method facilitates incorporating mesoscale interaction potentials.An outstanding outcome of that feature is the capability of modeling two-phase flows, which has always been a cutting edge area of study for scientists and engineers, without the necessity of paying the computational price of interface tracking or capturing methods which are applied in traditional CFD methods.
Various ideas of modeling vapor-liquid equilibrium (VLE) systems such as the color-gradient model [10], pseudopotential model [11,12], and free-energy model [13,14] had been introduced into the LB equation since its development as the successor of lattice-gas cellular automata (LGCA) [15].Among them, the pseudopotential model draws much attention due largely to its straightforward definition of particle interactions and ease of implementation.The original pseudopotential or Shan-Chen (SC) model has several shortcomings, such as being unstable at high density ratios (limited to density ratios of the order of 10) [12], being thermodynamically inconsistent with practical equations of state (EOSs) [16], VLE densities are dependent on viscosity [17], and interface thickness and interface tension cannot be tuned independently.Addressing such issues is the core motive for many researchers in the LB community.
Recent studies showed that surface tension can be controlled either by considering next-nearest interactions on the * b.chen@hw.ac.uk lattice [18,19] or incorporating an additional interaction force, or source term that affects only the interface region [20,21].Packing a practical EOS in a single interaction pseudopotential, which is a way of reaching high density ratios, was examined [16] for several cubic EOSs, and showed that the Carnahan-Starling (CS) EOS and Peng-Robinson (PR) EOS can enhance the stability of the model and reach the density ratio of the order of 1000.Nonetheless, all of the EOSs show thermodynamic inconsistency to different extents since only one specific pseudopotential is consistent [22,23].
As a workaround, it is found that the different forcing schemes (the ways in which the SC force is inserted in the LB equation) affect pressure tensor, thermodynamic consistency, and stability of the simulated system differently [17,24,25].These studies, adopting the CS EOS, showed that the forcing schemes proposed by Guo et al. [26] and He et al. [27] can similarly recover Navier-Stokes equations to the second order and they reproduce the predicted interface tensions well.Moreover, these forcing schemes [26,27] and exact difference method (EDM) [21] could effectively separate VLE densities from the change of viscosity.However, while EDM can reach to the density ratio of about 1000, the forcing schemes of Refs.[26,27] can reach density ratio of about 100.Wagner [28] compared the two methods of simulating nonideal fluids: the force method [29] and the pressure method [30].After successfully identifying spurious interfacial terms with the aid of Taylor expansion to fifth order and removing them from the pressure tensor in the force method, both pressure method and force method showed similar consistent results.Li et al. [31] showed that the forcing scheme proposed by Wagner [28] is identical to the one proposed by Guo et al. [26].Li et al. [24] successfully identified the effective part of EDM, which improves the instability, and implemented it into the forcing scheme proposed by Guo et al. [26] to reach higher density ratios (with the aid of the CS EOS), and then used a free parameter to approximately fit the mechanical stability on the solution of the Maxwell equal area.They further combined the scheme with multi-relaxation-time LB to improve the numerical stability [32].The study showed that the interface thickness can be controlled by the attraction parameter of the CS EOS.Lycett-Brown and Luo [33] analyzed the LB with a general forcing scheme with the aid of Taylor expansion to the third order and used the third-order errors to counteract the thermodynamic inconsistency.They concluded that the scheme can simulate arbitrary density ratios and increasing interface width can improve the thermodynamic consistency.However, it is inevitable for the methods to fit mechanical stability on the Maxwell construction for the entire coexistence curve or a particular thermodynamic state.
In a previous paper [34], we have shown that using single pseudopotential to recover practical EOSs can cause nonphysical interactions to happen in the simulation domain and, moreover, mechanical stability of the system is significantly affected by the lattice scaling and type and parameters of EOSs.Such a behavior can be found in Ref. [33], where they increase the attraction parameter of the EOS while favoring shortening the interface width; the vapor equilibrium density is considerably affected even at relatively low density ratios.Furthermore, we have introduced a multipseudopotential scheme which is thermodynamically consistent and can be initialized with the desired VLE state, independent of lattice spacing and the type of EOSs being recovered at least at low density ratios.
In this study, we propose a multipseudopotential scheme to incorporate most of the popular cubic EOSs into the LB equation consistently.While in the previous work multipseudopotential interaction (MPI) forces are determined to achieve a specific thermodynamic coexistence of interest, in this study they are analytically set to represent cubic equations of state including van der Waals (VW), CS, PR, and Soave-Redlich-Kwong (SRK) in a wide range of temperatures.It is shown that the scheme is numerically stable and high density ratios are achievable.Assigning each part of the EOSs (attraction and repulsion) to the consistent pseudopotentials is discussed.The EOS parameters that govern the interface width are identified and applied to improve the accuracy of VLE simulations at very high density ratios.The simulation results are supported by theoretical analysis which helps present them in reduced formats and link them to real world thermodynamic systems.

II. METHODOLOGY
The Bhatnagar-Gross-Krook lattice Boltzmann equation (LBE) for simulation of a flow field with a general forcing term is defined as where f i are particle distribution functions at position x and time t, i = 0 • • • q − 1 the indexes of neighboring nodes, e i the discrete microscopic velocity vectors to neighboring nodes, τ is the nondimensional relaxation time, f eq i the local equilibrium distribution function obtained from the Chapman-Enskog expansion of the Maxwellian to the second order at constant temperature [15,35], and F i the forcing term to be defined in Sec.II C. Most of the lattice Boltzmann models are developed based on Hermite quadratures [35] to reproduce the integrals of the equilibrium distribution function over the whole momentum space.However, the Gauss-Laguerre quadrature method can be employed to accurately calculate the integral of the equilibrium distribution function over the octant of the momentum space which helps to implement diffuse reflection boundary conditions at microscopic scales, i.e., high Knudsen numbers Kn 0.1 (for more details refer to [36][37][38]).For a flow field near its equilibrium state without internal or external forces, the momentum flux tensor can be obtained by employing the equilibrium distribution function and considering its isotropy on the lattice, where α and β are Cartesian coordinates, θ lattice temperature, and u velocity.The lattice temperature depends on the underlying lattice structure.For the case of the D2Q9 lattice, the discrete velocities are e 0 = 0, and e 1 = (ς,0), e 2 = (0,ς ), e 3 = (−ς,0), e 4 = (0,−ς ), e 5 = (ς,ς), e 6 = (−ς,ς), e 7 = (−ς,−ς ), e 8 = (ς,−ς ), where ς = x/ t is microvelocity and x is lattice spacing.This choice of lattice sets the lattice temperature as It can be found from Eq. ( 2) that such a system is governed by the ideal gas pressure.To model nonideal fluids and twophase systems, Shan and Chen [11] introduced a lattice force into the LBE, to describe the interaction potentials (pseudopotentials) between the fluid particles, where G is the amplitude of the force, N the number of neighboring nodes that interact with the node of interest, c i the vectors linking to those neighbors, and w the weight function that makes the interaction force dependent on the distance.The distance to the nearest neighbor is c, whose magnitude is equal to ς if the time step is set as t = 1.The SC force can be incorporated into the LBE through the forcing term and equilibrium distribution functions.

A. Multipseudopotential interaction scheme
The authors introduced the concept of the multipseudopotential interaction (MPI), in comparison with the singlepseudopotential interaction (SPI), as the pair interactions scheme to describe the hydrodynamic properties of real fluids, under the condition that each potential satisfies the thermodynamic requirement [34,39] F total = F (1) + F (2)  (5) where n is the number of pseudopotentials, G j the amplitude, and ψ j the consistent pseudopotential of the j th part of the force where C j is a constant and ε = (6e 4 − 2)/(6e 4 + 1) in which e 4 is the parameter can be obtained from fourth-order isotropy of the force on the lattice [22], For example, in the case of nearest-node interactions w(|c i | 2 ) = 0 for |c i | 2 > 2 and from fourth-order isotropy of the force on the lattice w(1) = 1/3, w(2) = 1/12, and e 4 = 1/3.Weights for |c i | 2 > 2 come into effect when a higher order of isotropy is needed (for more details refer to [18,40]).
Here, the interparticle potential is presented with the contributions from different interactions.Therefore, the interaction force between a considered node and its neighbors is composed of a set of subforces at various potentials, which are functions of particle densities, i.e., the inverse of the mean free path of particles.All interactions of the MPI are at the same range of the lattice, for example, the range to the nearest neighbors is selected in the simulations of this study.It is noteworthy that the concept and physics of multipseudopotential interactions are different to those of the multirange interaction (or pseudopotential with midrange interactions) model proposed by Sbragaglia et al. [18] and Falcucci et al. [19] in which different sets of forces are introduced to the scheme with different ranges (or cutoff radii) linking to the neighborhood nodes.Such a multirange interaction model provides capabilities for the pseudopotential models, such as flexibility of setting surface tension independent of EOS [18], achieving stable spraylike fluid [41,42], and increasing numerical stability [19,43].However, it has to be noted that the multirange interaction model [18] does not take the thermodynamic consistency into consideration, as all of its ranges utilize the same inconsistent pseudopotential ] where ρ 0 is a constant.

B. Analysis of MPI
It has been shown that the MPI scheme is intrinsically consistent with thermodynamics to the second order of spatial derivatives [34].The pressure tensor of the MPI system can be exactly expressed by the MPI force on a discrete lattice for the nearest and next-nearest interactions.In the case of D2Q9 lattice and continuum limit where all gradients are in the x direction, the contribution of all potentials of the MPI force to the normal component of the pressure tensor can be found in the x − y plane (to the second order) [34], where a 1 = 1 − 3e 4 , a 2 = 1 + 6e 4 , and ε = −2a 1 /a 2 .The system results in bulk pressure of In this equation of consistent pressure, G j , λ j , C j are free parameters; however, ε is constrained by the order of isotropy of the interaction force on the lattice [see Eqs. ( 6) and ( 7)].Therefore, as discussed in the previous study [34] the functions in the formats, such as ρ h and [ρ/(1 + Cρ)] h , where h is an arbitrary number, are not accepted in the MPI EOS while they are the necessary parts of well-known cubic equations of state including the VW, CS, SRK, and PR.In order to eliminate this constraint, in the following section the forcing scheme is modified to make ε flexible.

C. A more flexible MPI scheme
To incorporate the lattice force, shifting the particles' momentums in the equilibrium distribution functions was proposed by Shan and Chen [11].However, Guo et al. [26] demonstrated that the following forcing term is necessary to consider the discrete effects of the lattice: where the Einstein summation convention is adopted.This equation gives rise to the stress tensor in the continuum limit, as follows: where v is the actual fluid velocity, applied in the equilibrium distribution function, and defined as , Navier-stokes stress tensor is recovered [26] as, Li et al. [24] showed that considering the additional term F α F β /ψ 2 to the stress tensor can effectively modify ε which appears in the thermodynamic consistency condition (6).The leading term of Here, for the MPI, we use a similar term and set tensor C αβ as where s j are arbitrary constants.Substituting ( 14) into (10), we obtain As a result, the macroscopic approximation of the MPI-LB equation leads to where p is the ideal gas pressure.The excess terms on the right-hand side of ( 16) can be considered a part of the pressure tensor, where Eq. ( 13) is used.For the case of a one-dimensional planar interface, the component of the pressure tensor which is normal to the interface now takes the form [compare with Eq. ( 8)] Therefore, ε in Eq. ( 6) now becomes which is the relation derived to release the constraint of the MPI, as by which {ε j : j = 1,2, . . .,n} become flexible and independent of each other in the pseudopotentials of the MPI.The MPI EOS is, therefore, defined as As expected, in comparison with (9), the additional term in ( 14) modifies the thermodynamic condition only through ε j and the rest stays intact.In other words, ε j is a function of the free parameter s j .To obtain the interface shape, the equation (15) of our previous paper [34] becomes The surface tension of the two-phase MPI system can be presented by the mechanical route definition; in the case of a planar interface, it is calculated from σ = +∞ −∞ (p N − p T )dx, where p N is the local pressure normal to the interface and p T the local pressure tangential to the interface, where K j = −6e 4 + 12 tG j s j and ∂ x ρ is solved by (22).Considering Eq. ( 20), we have

D. Mapping cubic equations of state onto the MPI scheme
Here we show how the cubic equations of state such as the VW, CS, SRK, and PR can be implemented analytically and consistently into the MPI-LB system.
For the VW, the EOS, exists as the basis of many other EOSs of real fluids.The fluid's critical point (critical density ρ c = 1/3b, critical temperature T c = 8a/27b, and critical pressure p c = a/27b 2 ) is deemed to be the reference point for normalization.
According to (21), the functions such as ρ h or [ρ/(1 + Cρ)] h are thermodynamically consistent as pseudopotentials, and the parameters of pseudopotentials of the MPI can be adjusted; in a sense, the MPI EOS matches the VW EOS.We describe the procedures step by step as follows by means of three pseudopotentials (n = 3).Considering ( 21), firstly we remove the lattice ideal pressure (ρθ) from the MPI EOS; therefore, the first pseudopotential (j = 1) is assigned to it: From this equality we find G 1 = −2/3, C 1 = 0, ε 1 = 2, and λ 1 = 1/2.The second pseudopotential (j = 2) represents the repulsion term of the VW EOS, from which we obtain G 2 = 2T /c 2 , C 2 = −b, ε 2 = 2, and λ 2 = 1/2.The third pseudopotential (j = 3) demonstrates the attraction term of the VW EOS, then we find It should be noted that s j can be found from (20) and applied in Eq. ( 15).Consequently, it can be seen that each pseudopotential is assigned to the physics; one describes the dynamics of individual molecule, the second creates solely short range repulsion interactions, and the third mimics the long range attraction of particles.All parameters are listed in Table I for reference.Because the procedure of determining the MPI forces parameters for the VW EOS, conditions ( 26)- (28), is mostly the same as those to be used for the other EOSs, no details will be given in the following discussions.The Carnahan-Starling (CS) EOS, which modifies the hard sphere repulsion part of the VW equation, is an EOS widely used in the pseudopotential LB community.Here The CS EOS can be rewritten as Therefore, the CS EOS is in accordance with the MPI EOS.The relevant parameters of the MPI EOS are listed in Table I.It should be noted that the first pseudopotential has two parts: The first part is utilized to remove the lattice ideal pressure (ρθ) while the second part exerts the CS fluid ideal pressure (ρT ).
The Soave-Redlich-Kwong (SRK) EOS is the well-known two-parameter cubic equation of state and the first modern EOS which is widely applied to he design of hydrocarbontreatment plants [44].The SRK is the significant modification of RK EOS proposed by Soave [45]: where , and ω is Pitzer's acentric factor.α is determined from experimental vapor pressures of nonpolar substances.The properties of the critical point can be found from the fact that the first and second derivatives of pressure with respect to density at the critical point are zero.In contrast to RK EOS, the reduced format of the SRK depends on the acentric factor ω, which helps to treat nonpolar substances . By decomposing the second term of the SRK EOS, we obtain Now, the SRK EOS is well consistent with the MPI EOS.The relevant parameters of pseudopotentials of the MPI scheme are listed in Table I.
The Peng-Robinson (PR) EOS [46] is devised to overcome the weakness of the SRK in predicting liquid phase density [47]: where The reduced form of the PR EOS is The PR EOS can be rearranged into the following shape: Therefore, each term of the PR EOS now is in agreement with the MPI EOS.The corresponding pseudopotential parameters are listed in Table I.
It is worth mentioning that in contrast to many other studies [17,21,24,33], where thermodynamic consistency is sought by trial and error methods, curve fitting, or solving a nonlinear equation, the consistent MPI forces proposed here exactly represent the selected EOSs.Moreover, for the sake of brevity only the above-mentioned cubic EOSs are focused on, but the principle of the method is not limited to, the cubic EOSs; it is applicable to the EOSs in other function formats, such as polynomial and virial functions, and the versions of EOSs developed based on them (see Appendix A).
Using the MPI scheme, the interaction force between two nodes is calculated by summation of n (number of pseudopotentials) subforces.For the case of the two-dimensional nearest-node interactions, a node experiences 8n interactions.Therefore, increasing the number of pseudopotentials, we expect more computational costs.For example, the MPI code running the SRK or VW EOSs is faster than the one implementing the PR or CS EOS (see Table I); however, it is slower than a nonmodified SPI code [16] running the SRK, VW, PR, or CS.

III. RESULTS AND DISCUSSION
To verify the theoretical derivations and demonstrate how the proposed method works, simulations of planar interfaces are carried out.In Sec.III A the parameters and details of simulations are described.The thermodynamic consistency and interface shape are discussed in Sec.III B, and the mechanism of the scheme in controlling the interface width is studied in Sec.III C. The capability of the method in reproducing Maxwell construction curves is investigated in Sec.III D. The MPI droplet is thoroughly studied in Sec.III E.

A. Simulation setup
All simulations are run on a two-dimensional square lattice including nine velocities (D2Q9).The computational domain for one-dimensional planar interfaces contains 1 × 500 nodes and for two-dimensional cases contains 300 × 300 in most of the simulations, while, for the simulations of very wide interfaces a larger number of lattice nodes is utilized, for which the details will be given individually.The nondimensional relaxation time is set to unity τ = 1 unless otherwise stated.The nearest-node interactions are considered which set w(1) = 1/3, w(2) = 1/12, a 1 = 0, a 2 = 3.The MPI forces are calculated by use of Eqs. ( 5) and ( 19) whose parameters, for the selected EOSs, can be found in Table I.The internal forces are embedded in the LBE using Eq. ( 15), which reproduces Navier-Stokes equations to the second order.For the case of the planar interface, half of the domain is filled with liquid and the other half with the vapor.The periodic condition is applied at all boundaries.ρ l and ρ v (the densities of liquid and vapor) are initially set, in a sense, to ensure that the system is in the saturation state.To avoid an initial instability, a diffuse interface should be adopted at the beginning of simulations.We utilize the suggestion of Ref. [17], where W controls the interface width and, in all simulations, is set at W = 5; R is the distance between the interface and the reference point; and r is the variable distance from the reference point.For the case of the planar interface, it can be one of the boundaries parallel to the phase interface.Unless otherwise stated, the lattice spacing parameter is set to unity c = 1 and EOS parameters for all types are set at a = 0.01 and b = 0.2.The simulations are run for at least 10 5 steps and the pressure tensor profiles in VLE systems are monitored to verify that the equilibrium state is reached.

B. Thermodynamic consistency
To investigate the thermodynamic consistency, the first set of simulations is performed to measure the total normal  component of the pressure tensor, p 0 = P xx + ρθ.Based on the mechanical stability of the flat interfaces, p 0 must be constant across the phase interface (for more details refer to [28]).P xx is calculated from Eq. ( 18) where both first and second spatial derivatives of pseudopotentials are computed by use of the second-order central difference approximations.The results are shown in Fig. 1.For the selected EOSs and specified temperatures the errors are small and less than 0.15% which validates the mechanical stability of the proposed scheme.The errors increase with the decrease of temperature which will be discussed in Sec.III D.
The gradient of density along the phase interface can be calculated analytically from Eq. ( 22) as the right-side integral has an analytical solution.Then, by numerically solving x = ρ (ρ v +ρ l )/2 (∂ x ρ) −1 dρ, the spatial profile of density can be achieved.Figure 2 depicts the phase interfaces for the selected EOSs and the comparisons of the simulation results with the theoretical predictions.The main graphs of Fig. 2 focus on liquid phase densities and the inset graphs emphasize vapor phase densities; the simulations are in good agreement with the theory which means the interface shape can be ably predicted.
When a high density ratio system is desired an initial guess of Eq. ( 38) may not be of help and the simulation crashes as a result of instabilities generated at the interface.Therefore knowing the equilibrium interface shape in advance can effectively enhance the starting point stability of the system in practical cases.

C. Scalability and interface width control
Herein, for the sake of brevity we merely focus on the SRK EOS and investigate EOS parameters thoroughly.If we divide both sides of Eq. ( 22) by p c /ρ c and use the MPI SRK pseudopotentials of Table I, the equation could be normalized.Therefore, after some lengthy but simple algebra, we obtain It can be seen that the right-hand side of (39), which is the Maxwell construction integral, is in the reduced format which means it is independent of EOS parameters a and b, and lattice scaling parameter c.Therefore the reduced VLE densities are not affected by change of those parameters.It is understood that the terms, which are multiplied by the gradient of density, are impacting the interface shape.On the left-hand side of ( 39), all terms are in reduced format as well, except the first term in the brackets which is a function of ratio parameter FIG. 3. The phase interface is widened by decreasing χ = bc 2 /a for a particular VLE, the SRK EOS at T R = 0.80.Parameters are set at b = 0.2, c = 1, and a = 0.001, 0.01, 0.09, 0.25.The theoretical profiles are obtained from Eq. ( 22).The theoretical profile of density at χ = 0.8 is incomplete, which is explained in Sec.III C. χ = bc 2 /a.Therefore, it is expected that this parameter, χ , is the sole parameter responsible for the change in the width of the phase interface.As the VLE system is independent of χ in Eq. ( 39), the sound speed in liquid and vapor is independent of χ as well and is consistent with the EOS.It is worth mentioning that by changing χ from 200 to 2.2 the interface is perfectly predicted by (22); however, for the case of χ = 0.8 a very thin interface is partly predicted from vapor phase to liquid phase (see Fig. 3).This can be understood by checking the right-hand side of Eq. ( 22), which is the Maxwell construction integral and is always negative for ρ v < ρ < ρ l .In order to have a meaningful gradient of density, the left-hand side of Eq. ( 22) should be negative throughout the given density interval.Therefore, in the case of χ = 0.8, we cannot predict interface shape near the liquid phase as the left-hand side of Eq. ( 22) becomes positive there.Moreover, slight oscillations of density near the liquid phase of the simulated system are observed.We only experienced the partial prediction of the interface in very sharp and thin interfaces.
In χ another parameter, c, being the numerical resolution to capture the physical details, makes a contribution to the interface width.For the pseudopotential LB models, it is predicted that the interface width is varied directly as the lattice scaled by means of c.This behavior might not be of interest as the interface width of the fluids such as liquid water with its vapor is in the scope of molecular scale and LB simulations aim at microscale phenomena; therefore with scaling an LB FIG. 4. The phase interface preserves its shape through lattice scaling by keeping χ = 20.The theoretical profiles are obtained from Eq. ( 22).simulation domain in the interval of microscale, a change in the interface shape should not be seen.The results shown in Fig. 4 are those from the simulations with the domain being scaled at c = 0.1, 1, 10 but parameters a and b are utilized to keep χ = 20.It is demonstrated that the same interface shape is obtained.The simulation results are supported by the prediction of theory (39).
It should be noted the VW, CS, and PR EOSs have the same behavior as the SRK EOS and a normalized governing equation similar to (39) for the phase interface.Therefore, from the simulations and theoretical analysis, the parameter χ is identified as the only parameter that governs the interface width for the selected cubic EOSs.
In addition to the discussion on the effects of EOS parameters a and b, and lattice spacing parameter c on the numerical simulations of the phase interface, it is also important to investigate the effects of those parameters on the VLE state of interest.Such a discussion, in our previous study [34], through examining the single pseudopotential scheme recovering cubic EOSs showed that the mechanical stability condition is remarkably sensitive to the changes in these parameters.Figure 5 shows the errors for vapor density deviation of simulations, using the MPI scheme proposed in this study, Maxwell construction due to the change of a, b, and c.It can be seen from Fig. 5 that error in simulating vapor density is increasing with the increase of lattice scaling parameter, c, increase of a, and decrease of b.However, all the simulations errors are less than 0.6% for the SRK EOS at T R = 0.80 which means at relatively low density ratios (high temperatures) the thermodynamic state is satisfactorily independent of the change of those parameters.

D. Liquid-vapor saturation curves
To further test the proposed scheme, VLE systems are simulated from a temperature near the critical point to the lowest stable temperature for planar interfaces whose widths are controlled via parameter χ .As the interface width theoretically goes to infinity in the case of diffuse interfaces, the nondimensional interface width is defined as l 2% = |x ρ=0.98ρ l − x ρ=1.02ρ v |/ x.The one-dimensional domain is long enough to make sure the vapor and liquid densities reach plateau profiles.The VW, CS, SRK, and PR EOSs are examined in the sense that parameter χ is changed at each temperature to have constant interface widths of l 2% = 20, 30, 60.The χ values are listed in Appendix B.
Figure 6 compares the MPI simulations of VLEs with the predictions of the Maxwell construction.It can be found that all EOSs almost perfectly simulate the liquid phase.Regarding the vapor branch, in all cases we have satisfactory results at high temperatures and different interface widths.In the case of short interface width l 2% = 20 the increase of error is more obvious by decreasing temperature.Widening interface width l 2% = 30, 60, the error significantly decreases, which is more apparent at lower temperatures.Therefore, at lower temperatures a higher-resolution interface is necessary for the MPI to be in accordance with the thermodynamic requirement.
The errors due to deviation of the vapor and liquid branch from the Maxwell construction are plotted vs temperature for the selected EOSs in Figs.7(a) and 7(b).The errors of the liquid branch for all cases are less than 0.01% and generally decrease to a local minimum and, after a rise, decrease continuously again by decrease of temperature.In general, the errors at the liquid branch decrease as the interface widens at a given temperature; however, there are some exceptions such as SRK EOS at T R = 0.80.With regard to the vapor branch shown in Fig. 7(a), in all cases the error increases but faces a fall immediately and then again increases limitlessly when temperature recedes from the critical point.Widening the interface, the vapor branch error can be effectively reduced but it might not be very helpful near the local minimum errors, such as for VW and CS at T R = 0.7, and SRK and PR at T R = 0.8.
When checking the liquid branch error vs density ratio in Fig. 7(d) the local minimum errors can be identified, which are located around a density ratio of 20 for all EOSs (except PR EOS) at l 2% = 20, 30.The local minimum errors move to a density ratio of 500 when interface width is broadened to l 2% = 60.Concerning the vapor branch error shown in Fig. 7(d), for all cases the local minimum errors are at a density ratio around 25. Beyond the density ratio of 50, regardless of the type of EOS, the errors are only a function of the interface width whose increase reduces the error significantly.
Considering Eqs. ( 39) and ( 22), and knowing different EOSs having different VW-loop shape, the geometrical shape (or ∂ x ρ) of the interfaces should vary even at the same interface width.Therefore, we expect the different EOSs at the same density ratios and interface widths to show relatively different behavior.
The current MPI scheme is basically constructed on the previous idea [34].However, from the stability viewpoint, in comparison to that, the current scheme provides an independent parameter χ , to control the interface width which helps us straightforwardly enhance the stability of simulations.Therefore, the maximum achievable density ratio from order 10 in the previous paper could successfully reach to 10 6 .

E. Circular droplet
MPI drops are analyzed regarding measuring surface tension, spurious velocities, isotropy, dynamic behavior, and assessing the stability dependent to the relaxation time.
To understand the effect of χ on surface tension, the SRK EOS is considered.Therefore, by substituting the SRK MPI parts from Table I into Eq.( 23) and normalizing surface tension by σ R = σ/(cp c ), we obtain (40) where K 1 = K 2 = K 3 = −5 is used.All terms depend only on SRK constant parameters except c dρ R /dx and the first term inside the brackets which is a function of χ = bc 2 /a.In fact, while widening the interface, χ helps attraction interaction forces around the interface.
We make use of the Laplace test to validate Eq. ( 40).For measuring the surface tension in the Laplace test, it is necessary to find the accurate radius of drops in simulations.We consider a domain where density at the center (r = 0) is ρ l and it continuously declines to reach ρ v at r = R; we assume it is equivalent to a domain where a drop of radius R m with a homogeneous density of ρ l is immediately in contact with a vapor phase with density ρ v .By use of the mass equality, we obtain the mean radius . The simulations are run in a 300 × 300 domain and SRK EOS at T R = 0.80 is employed.Figure 8 shows the results of the Laplace test for drops with different radii.The surface tensions of the drops with interface widths χ = 2.5, 10, 25 are satisfactorily predicted by Eq. ( 40).As can be seen, surface tension increases with a slow rate proportional to χ ; thus we advise χ is better used for adjusting the interface width rather than tuning surface tension, which will be investigated in the future.
As had been discussed in the previous study [34], the MPI scheme reduces the spurious velocities as the MPI is thermodynamically consistent.In this study, the effects of interface width, which can be adjusted freely by χ , on the generation of spurious velocity are further discussed.Two two-dimensional (2D) drops with a density ratio of 1000 are set up using CS EOS at T R = 0.486 and VW EOS at T R = 0.37.Similar to Ref. [33], we define the distortion ratio Simulations are run for moving circular drops to study the dynamic behavior and assess the Galilean invariance of the proposed scheme using SRK EOS.For all three cases, the MPI To access the Galilean invariance, the distortion ratio is defined, DR 90 , as the ratio of longest diameter of a drop to the diameter perpendicular to the longest one [48].The MPI SRK drop experiences three stages: the rest, acceleration, and the long rest for each case run.In the beginning, the drop is kept to stay at the rest state for 5 × 10 4 steps to ensure it reaches the equilibrium shape.Then the drop is accelerated diagonally by a x = a y = 0.0001 for time period 500 t (case 1), 1000 t (case 2), and 2000 t (case 3), respectively.The accelerated drop is then released for 5 × 10 5 time steps, for all three cases, to remove the transit effects and enable the drop along with its surroundings to move diagonally at constant velocities.The simulation results are discussed as follows.
At equilibrium state, the MPI SRK drop keeps a perfect circle drop with DR 90 = 1.0 and has density of ρ R l = 3.195 in contact with vapor density of ρ R v = 5.505 × 10 −3 .The MPI drop and the pattern of its spurious velocities are shown in Fig. 10(a).
As the drop and its surroundings reach higher velocities, the shape eventually distorted from the circle to DR 90 = 1.004 (case 1), DR 90 = 1.018 (case 2), to DR 90 = 1.087 (case 3), which is demonstrated in Figs.Considering case 3, which is a critical case demonstrated in Fig. 10(d), the drop is disfigured noticeably and the system highly violates Galilean invariance as the fluid is unable to damp spatial velocity differences by time.This means that the proposed MPI scheme is well Galilean invariant for |u|/c < 0.15 but violates Galilean invariance beyond this ratio.
As had been discussed in [49][50][51], the mechanism of violation of Galilean invariance for LBM is due to the insufficiency of the equilibrium distribution function, which is Maxwellian expansion to the second order at constant temperature.It introduces the cubic velocity term, τ ∇ • (ρuuu), to the viscous momentum flux of the LB system.In this study, the equilibrium distribution function is employed without any modification; therefore the term, τ ∇ • (ρuuu), is the only term which violates the second-order Galilean invariance.If |u| 2 /θ or |u|/c is very small the cubic velocity term is negligible in comparison with the Navier-Stokes viscous term.That is why when the bulk velocity is increased, the spatial velocity differences are not dissipated by time and the drop shape is distorted.
The MPI SRK drop is tested to obtain its minimum stable relaxation time τ min .The drop, at different temperatures and interface widths, is located at the center of the computational domain.The relaxation time is set at τ = 1 initially; after the drop reaches its equilibrium state, τ is gradually, and in a quasiequilibrium process, decreased to achieve new lower viscose systems.The process is monitored to record the point before the system crashes at which the relaxation time is taken as the minimum possible relaxation time τ min .As shown in Table II, the decrease of the temperature and interface width weaken the τ -related stability.

IV. CONCLUSION
To summarize, the MPI scheme is extended, which is based on the original SC model, in order to eliminate the numerical errors as a result of the discretization effect of the lattice when cubic equations of state are implemented.The thermodynamic consistency condition has been made more flexible by modifying the forcing scheme by which in addition to cubic EOSs other types such as a virial EOS can be employed analytically and consistently in pseudopotential models.In such a way, different pseudopotentials can be identified by the  The most basic and popular cubic EOSs are adopted including the VW, CS, PR, and SRK to demonstrate the validity of the proposed MPI scheme by performing a set of simulations of planar interfaces.For all the cases the normal component of the pressure tensor is satisfactorily constant along the flat interfaces.
The LB system is analyzed in a reduced format which can be compared with real physical systems.The equilibrium interface shapes are predicted from mechanical stability condition and found in good agreement with the simulations.Interface width could be systematically adjusted with the aid of the EOS parameters, while they can positively affect the accuracy of the VLE densities.Such a feature helps to control the interface width in the case of grid refinements.
The proposed MPI scheme provides stable two-phase systems even for very high density ratios.The liquid branch of the Maxwell construction curve is almost perfectly achieved.
Regarding the vapor branch, the errors are negligible and small at low and midrange density ratios (less than 1000) but, at higher density ratios, exponentially grow.The errors can be suppressed by broadening the interface width, i.e., increasing the interface resolution.
We have shown that the isotropy of the circular drops are satisfactorily preserved in the process of shortening the interface width.An MPI SRK drop is moved diagonally in a periodic domain which showed its Galilean invariant at relatively low bulk velocities in addition to being stable.The minimum available relaxation time for the stationary MPI drop is obtained which increases with the decrease of temperature and interface width.
a = a T 2 c /P c , b = b T c /P c , bρ c = c , a = 0.457 236, b = 0.077 796, c = 0.253 077, and α is similar to that in the SRK equation but m is correlated with the aid of vapor-pressure data from normal boiling point to critical point: m = 0.374 64 + 1.542 26ω − 0.269 92ω 2 .(

FIG. 1 .
FIG. 1. Error due to the change of the normal component of pressure tensor along flat interfaces at different temperatures for the selected EOSs: (a) VW, (b) CS, (c) SRK, and (d) PR.Pressure is normalized with the aid of the saturation pressure predicted by the Maxwell equal-area rule.In graph (a) EXP-MPI symbols represent the MPI results of our previous paper[34] for reference only.The errors are smoothly changing from one phase to the other one.
Four simulations are run for the different attraction parameter values a = 0.001, 0.01, 0.09, 0.25 while the SRK EOS is set at T R = 0.80, b = 0.2, and c = 1.The results are shown in Fig. 3.It is found that the increase of parameter a (decrease of χ ) can effectively reduce the interface width.Another set of simulations have been run by keeping a = 0.01, c = 1 and changing b to have the same χ values as the previous the results are exactly the same as those of Fig. 3 where the interface width increases as parameter b rises.

FIG. 5 .
FIG. 5. Error in vapor density due to change of EOS parameters a and b, and lattice spacing parameter c.The SRK EOS at T R = 0.80 is applied.
10(b)-10(d).The changes in densities are neglectfully small, in comparison with those of the equilibrium case, for all cases smaller than 0.03% for liquid and 3.2% for vapor.The flow pattern, however, is more sensitive to the changes in bulk velocity whose increase from |u b |/c = 7.071 × 10 −2 (case 1) to |u b |/c = 1.414 × 10 −1 (case 2) and then to |u b |/c = 2.818 × 10 −1 (case 3) causes the flow to deviate more from the uniform distribution.

w 2 attraction
and repulsion parts of cubic EOSs without the need for curve fitting or trial and error methods.
interested readers, values of the χ parameter used in reproducing the Maxwell construction curves at different interface widths demonstrated in Sec.III D are shown in Tables IV and V.

TABLE I .
The MPI forces parameters which reproduce the selected EOSs in the LB system.
EOS j th pseudopotential

TABLE II .
The minimum possible relaxation time of the MPI SRK drop at different temperatures and interface widths.The interface width is measured from the flat interface simulations and can be obtained by setting χ at a given temperature found in Appendix B.

TABLE III .
The MPI forces parameters which reproduce virial and VW-like EOSs in the LB system.