Fluid–fluid interactions in pseudopotential lattice Boltzmann models: Effects of model schemes and fluid properties

This article demonstrates the characteristics of interparticle forces and forcing schemes in pseudopotential lattice Boltzmann simulations. The Yuan–Schaefer (YS), multipseudopotential interaction (MPI), and piecewise linear methods are examined as techniques of equation of state (EOS) inclusion in pseudopotential models. It is suggested here that it is important to have an understanding of the interparticle forces generated by the models in order to obtain good quality results. Poor choice of parameters can lead to generation of unphysical interactions. The piecewise linear method is found to perform well and to decouple parameters. It decouples the density ratio from the surface tension and from the collision operator relaxation rates. It is proposed that the decoupling occurs due to generation of lower values of high‐order error terms in the interfacial region by the piecewise linear EOS. In general, the multiple‐relaxation‐time (MRT) collision operator should be combined with the Huang–Wu forcing scheme for simulating high values of surface tension and with the Li–Luo method for simulating low values of surface tension. It is found that reducing kinematic viscosity is more detrimental to the stability of the simulations than increasing the density ratio. Introducing a kinematic viscosity ratio between the phases practically eliminates the influence of density ratio on spurious velocities. The factors affecting stability of dynamic simulations are examined. It is found that they have the following hierarchy from the greatest impact to the least: kinematic viscosity ratio between the phases, bulk viscosity, method of EOS inclusion, and reduced temperature/density ratio.


INTRODUCTION
The lattice Boltzmann method (LBM) 1 is a mesoscopic numerical method for simulating fluid dynamics. It has been developed into an attractive alternative to traditional continuum-based computational fluid dynamics (CFD) 2 and it has the potential to be useful to the engineering community for the simulation of physical and chemical processes. 3 LBM fills an important niche in the investigation of complex phenomena in fluid dynamics which are challenging for traditional CFD methods, especially multiphase and/or multicomponent flows such as the simulation of interfacial dynamics induced by fluid interactions. 4,5 Another area of application of the LBM is in the treatment of complex (e.g., discontinuous) interfaces and tortuous geometries. 4 The strength of LBM derives from its mesoscopic nature which places it between macroscopic and microscopic fluid simulation methods. 6 Kinetic origins and statistical nature allow it to retain physical insight without sacrificing efficiency. The seminal work of Gunstensen et al. 7 initiated the evolution of multiphase LBM models with the development of the color-gradient model. The method was based on simulating two different-colored (red and blue) fluids and its algorithm contained two collision operators and a recoloring step. Pressure inclusion in the model violated Galilean invariance requiring introduction of correction terms to fix the simulated macroscopic equations. 6 Significant pieces of research work dealing with this model include simulations of miscible flows 8 and tackling of Galilean invariance for simulation of multiphase flow in porous media. 9 The computational cost of the model can be described as average. 4 Despite being the first multiphase LBM model, the color-gradient approach is still in use and is still being developed. Notably, Montessori et al. 10,11 recently extended the model to include near contact interactions enabling studying of highly internal phase emulsions and foams.
The pseudopotential model, 12,13 which is the subject of this study, was the second multiphase LBM model to be published in the literature. It is a simple model with a "bottom-up" structure in which interparticle forces are introduced to bring about multiphase behavior by giving rise to a non-ideal pressure tensor. 14 An interesting feature of bottom-up models is that they work on a kinetic basis, just like the LBM itself, which may explain their algorithmic simplicity. The pseudopotential model has garnered widespread acceptance within the diverse research communities using the LBM. 15,16 It is particularly suited to dynamic problems, which happens to be the field where the LBM is most powerful. 6,17 The pseudopotential model has a sound physical basis with connections to the second virial coefficient and the radial distribution function in the Enskog kinetic theory. 18 The two main problems associated with the model are spurious velocities and thermodynamic inconsistency. The original pseudopotential model 12 has been extensively expanded and improved since its inception to expand the physical envelope which can be simulated accurately and efficiently.
The next multiphase LBM model to be developed was the free-energy model. 19,20 It represents a different type of multiphase method than the pseudopotential model. The starting point of the model is the thermodynamics of the system in the form of the free energy rather than the intermolecular interactions. Thus, macroscopic properties are the starting point of the model leading to the "top-down" structure. 14 A non-ideal pressure tensor is introduced by modifying the equilibrium distribution function. 6 As in the color-gradient model, this necessitates introduction of correction terms to restore Galilean invariance. 6 The phase-field models 21 are another group of multiphase LBM models. The interface dynamics simulated by these models are governed by the Cahn-Hilliard equation. 6 A shortcoming of these models is their computational cost which is greater than the computational costs of the color-gradient, pseudopotential or free-energy models. 4 A factor working in favor of the phase-field method is its suitability to dynamic cases which stems from the fact that, unlike the color-gradient or free-energy models, no correction terms to restore Galilean invariance are required. 6 The entropic model 22 is a recent addition to the list of multiphase methods for the LBM. Entropic models find their application mainly in the simulation of higher Mach number single-phase flows including supersonic flow. 23 The principle of the method is based on increasing the kinematic viscosity of the unstable regions according to the H-theorem. 24 The kinematic viscosity adjustment is highly localized around the critical regions to dampen the instabilities in the form of pressure waves. 24 It has been recently reported that the entropic model can be combined with the pseudopotential model for multiphase simulations. 24,25 Combination of the pseudopotential model with the newest LBM methods highlights continued interest of the LBM community in the pseudopotential model.
It is evident that a wide range of multiphase models for the LBM exists in the literature. This state of affairs exists due to a lack of consensus on what is the best strategy for simulating multiphase fluids. 2 It is the objective of this study, taking the pseudopotential model as an example, to elucidate the roles of the model schemes and key parameters in the simulation of two-phase fluids. The investigated schemes include methods of equation of state (EOS) inclusion and forcing schemes.
The criteria for the assessment of the schemes and the parameters are based on three key performance indicators (KPIs). The first KPI is the physics of real fluids dictated by the EOS. The second KPI is the interfacial interaction which is assessed based on surface tension, spurious velocities, and interfacial thickness. Finally, independent adjustment of density ratio, kinematic viscosity, and bulk viscosity is investigated.
This article is organized as follows: The LBM is described in Section 2. The collision operators are presented in more detail in Section 3. The pseudopotential models are discussed in-depth in Section 4. The results of numerical simulations are included in Section 5 in order to investigate the nature of fluid-fluid interactions in pseudopotential models. Finally, conclusions are put forward in Section 6.

LATTICE BOLTZMANN METHOD
The LBM is a discretized lattice model of the continuous Boltzmann equation with the original multi-body collision operator simplified to allow efficient simulation and analysis. The continuous Boltzmann equation takes the following where v is the fluid particle velocity, F ext is the body force, p is the momentum, f is the distribution function of fluid particles in the spatial and momentum spaces, subscripts 1 and 2 denote particles, apostrophes distinguish pre-and post-collision states, g is the relative speed, is the differential cross section and Ω is the solid angle. Discretization is carried out using the Hermite series expansion 26 and the method of characteristics. 14 LBM algorithm consists of equilibrium distribution, streaming, boundary conditions, collision operator, and a multiphase model in the case of multiphase simulations. The equilibrium distribution is a second-order truncated Maxwell-Boltzmann distribution, whose discretized version takes the following form 27 : where w i is a lattice-dependent weighting, is the density, u is the observable fluid velocity, c i is a particle's lattice velocity, and c 2 s is the lattice sound speed squared. Second-order truncation in discretization of the equilibrium distribution leads to cubic defects which result in loss of Galilean invariance. 28 However, second-order accuracy is sufficient for most applications as evinced by the widespread applicability of the LBM. A wide-range of boundary conditions for the LBM exists in the literature [29][30][31][32] ; however, this article only investigates the fluid-fluid interactions. Analysis of the LBM is carried out using the Chapman-Enskog analysis, 33 although alternative methods exist. 34

COLLISION OPERATORS
The simplest and most popular collision operator is the single-relaxation-time (SRT) or Bhatnagar-Gross-Krook (BGK) model 35 : where is the relaxation time. is a kinetic parameter and it is related to kinematic viscosity through the following relationship: where c 2 s (lattice sound speed squared) is equal to 1/3 for the D2Q9 velocity set. A natural choice of τ is 1.0. Setting to 1.0 results in full-relaxation and direct decay of distribution function to its equilibrium. 14 When τ is greater than 1.0, the system is under-relaxed and decays exponentially towards the equilibrium state. 14 Over-relaxation and oscillation of distribution function toward equilibrium occur when τ is less than 1.0. 14 Loss of stability occurs when the distribution function enters negative values, which is why loss of stability can occur at low values of viscosity. The value of has to be greater than 0.5, as apparent from Equation (4). The SRT model has the advantages of simplicity and efficiency, but it can suffer from instability and inaccuracy. When the simple bounce-back boundary condition is used with the SRT operator, the location of walls is τ-dependent, 36 which is problematic for permeability studies in porous media. The force schemes used with the SRT collision operator can introduce unphysical viscosity-dependence of physical parameters. 37 Alternatively, multiple-relaxation-time (MRT) models can be employed 38 : MRT models improve the stability and accuracy of LBM simulations. 34,39,40 MRT models work by relaxing the moments at different rates. This necessitates conversion between distribution functions and moments. The conversion is carried out using the transformation matrix and its inverse. Equilibrium moments are calculated from equations containing only density and velocity terms. The transformation matrix can be constructed using Gram-Schmidt orthogonalization or Hermite polynomials. 14 The relaxation matrix is diagonal when the MRT model is constructed using Gram-Schmidt orthogonalization and takes the following form for the D2Q9 velocity set: The first, fourth, and sixth relaxation rates correspond to mass and momentum conservation. Mass and momentum are conserved in collisions. Some authors 34 advocate setting the rates corresponding to the conserved moments to an arbitrary value, for example, zero; however, to preserve the benefits of trapezoidal integration, they should be set to a nonzero value. 41 e is responsible for bulk viscosity relaxation, corresponds to kinetic energy square, q is called the energy-flux, and v is the kinematic viscosity mode. 36 Only the bulk viscosity and kinematic viscosity relaxation rates are hydrodynamically-relevant. The remaining modes are tuned to achieve accurate simulations. e , , and v are symmetric modes and q is an antisymmetric mode. 36 Magic parameter is a convenient tool for tuning MRT models. 42 The magic parameter used in this study was calculated in the following manner: There appears to be a lack of consensus on tuning the MRT models in the literature. It has been put forward that the optimization of parameters in the MRT collision operator does not obey universality criteria, leading to stability improvement being limited to the linear regime. 43 This would be in contrast to the stability improvements in both the linear and nonlinear regimes offered by the entropic schemes.
Setting the magic parameter to 1/12 to minimize third-order errors gives the best results for multiphase simulations. 44 However, setting the magic parameter to even lower values has been reported. For example, Li and Luo 45 set the anti-symmetric relaxation rate to 1.1 at a kinematic viscosity of 0.025 resulting in the magic parameter equal to 0.0307, in order to simulate a droplet splashing on a thin liquid film at a Reynolds number of 500 and a density ratio of 500. Increasing the value of bulk viscosity to multiple times the value of shear viscosity can reduce spurious currents and increase stability by attenuating pressure waves. 14,40,46,47 In this article, the bulk viscosity was set to be five times greater than the kinematic viscosity, unless otherwise stated. The bulk viscosity range used in this paper coincides with the bulk viscosity values for water vapor. 48 The anti-symmetric mode was adjusted to give the magic parameter equal to 1/12. The kinetic energy square relaxation ( ) was set to the same value as the bulk viscosity relaxation.
The cascaded method 49 is an alternative collision operator that was reported to perform as well as the MRT model for multiphase simulations. 50 The cascaded model is not investigated in this article. Another alternative is the tworelaxation-time (TRT) model, 51 which has only two relaxation rates, that is, one for the symmetric modes and another for the anti-symmetric mode. The TRT model offers increased stability and accuracy with the simplicity of the BGK model. 42

PSEUDOPOTENTIAL LATTICE BOLTZMANN METHOD
The premise for the pseudopotential model is the observation that multiphase fluids can be simulated by introducing nearest-neighbor intermolecular force, obtained from functions of density, into the LBM 12 : where G is a parameter controlling the strength of the intermolecular interactions. Hence, phase segregation due to attraction of like molecules, in the case of single-component mixtures, can be introduced. Density is a convenient parameter since it is a moment of the distribution function and one of the main observables in LBM simulations together with velocity. The original pseudopotential function employed in the Shan-Chen model takes the following form 12 : This form of the pseudopotential is unsuitable for use in simulations involving high density ratios. A number of improvements have been developed to overcome this limitation. Higher order of isotropy models 52 are not investigated, since increasing the interface thickness was found to be more efficient than increasing isotropy. 6,52 Expanding the pseudopotential interactions beyond the nearest-neighbors also requires modification of the boundary conditions making application more difficult. 52

Methods of equation of state inclusion
Yuan and Schaefer 53 discovered that high density ratios can be simulated by inserting an EOS into the pressure term in the square root form of the pseudopotential, a form proposed originally by He and Doolen 54 and obtained from rearranging the pressure equation: In this form of the pseudopotential, G loses its physical meaning. It is simply set to −1, in order to keep the term inside the square root positive, and its role is taken over by the terms within the EOS. The two widely used cubic equations of state in pseudopotential models, that is, Carnahan-Starling (CS) 55 and Peng-Robinson (PR) 56 equations of state, take the following forms, respectively: where For example, the a parameter in the cubic equations of state represents attractive forces and b represents repulsive forces. 53 Yuan and Schaefer suggested setting a to 1.0, b to 4.0, and R to 1.0 in the CS EOS. 53 They also suggested setting a to 2/49, b to 2/21, and R to 1.0 in the PR EOS. 53 Yuan and Schaefer (YS) tested the Shan-Chen, Van der Waals, Redlich-Kwong, Redlich-Kwong Soave, PR, and CS equations of state and obtained the best results using the PR EOS. 53 In this article, the CS EOS is used in preference of the PR EOS because it has one less parameter. The a parameter in the CS EOS offers control of interface thickness. 57 It is responsible for attractive forces in the EOS and reducing its value from the 1.0 proposed by Yuan and Schaefer expands the interface. 57 Discussion of interface thickness is subjective, unless the same cut-off percentages of bulk densities are used to calculate interface thickness to compare different methods.
An alternative method for including an EOS in the pseudopotential model is a piecewise linear EOS. 58 In a piecewise linear EOS, the following equations are used to simulate non-ideal equations of state 58 : From the above equation, it can be seen that the piecewise linear method offers five free, adjustable parameters, that is, V , M , L , 1 , and 2. 58 The first step is to set the pressure slopes in the vapor branch ( V ), unstable branch ( M ), and in the liquid branch ( L ). Then the spinodal points (i.e., 1 and 2 ) can be obtained from two equations, one for the mechanical equilibrium and the other for the chemical equilibrium, given by the following respective equations 58 : Multipseudopotential interaction (MPI) 59,60 represents another class of the pseudopotential method. Multiple thermodynamically consistent pseudopotentials are introduced in order to allow simulation of high density ratios with the SRT collision operator. The Guo forcing scheme 61 for the BGK operator is modified to allow flexible selection of j for each multiple pseudopotential. 60 Multiple pseudopotentials can be described by the following equations 60 : Clearly, MPI has a significantly greater number of parameters, which are obtained from decomposition of cubic equations of state. Values of the parameters for simulation of the Van der Waals, CS, PR, and Soave-Redlich-Kwong equations of state can be found in table I in Reference 60. MPI has recently been combined with the MRT collision operator. 62 The coupling has been achieved by adapting the Li-Luo 57 forcing method.

Multiple-relaxation-time forcing schemes
The original Shan-Chen model 12 was developed for the SRT collision operator. The pseudopotential force can be introduced into LBM models using different forcing schemes. For example, it can be introduced into MRT models using forcing schemes developed specifically for the MRT collision operator. McCracken and Abraham 41 proposed a forcing scheme for the MRT collision operator which does not introduce unphysical viscosity-dependence. 63 The MRT forcing schemes discussed in this paper are based on this scheme and all of them can be described by a general equation. The individual terms in the general equation are given in Tables A1 and A2 in Appendix A and the general equation can be described as follows: Two improvements to the MRT forcing model are analyzed and discussed in the following sections. They are referred to in this paper as Li-Luo 57 and Huang-Wu 44 models. In the Li-Luo method, the value of in the mechanical stability condition is tuned using . The equation for setting takes the following form for the D2Q9 velocity set 57 : Li and Luo 64 also proposed a separate method for surface tension adjustment, which can be combined with the method for adjusting thermodynamic consistency. Surface tension is adjusted according to 1 -. Hence, setting to 0 results in unmodified surface tension and results of simulations with and without the C source term are the same. Increasing the value of reduces surface tension. Li and Luo investigated the values of ranging from 0 to 0.99. 64 The Huang-Wu 44 method, unlike the Li-Luo methods, combines thermodynamic consistency adjustment and surface tension adjustment into a single term. Thermodynamic consistency is tuned using the sum of k 1 and k 2 and surface tension is tuned using k 1 only, as follows 44 : where c is the surface tension coefficient. Huang and Wu 44 investigated setting c from 0.1 to 2.0. The Huang-Wu 44 method introduces modifications at the third-order in the Chapman-Enskog analysis, whereas the Li-Luo method is at the second-order. Recent studies illustrate the importance of third-order analysis to the pseudopotential models. 44,65,66 In the single-pseudopotential (SPI) models employing the square root form of the pseudopotential, is simply a lattice and forcing dependent numerical parameter. Its value does not have physical meaning and it needs to be adjusted on a case-by-case basis. In the MPI model, on the other hand, is an EOS parameter.

NUMERICAL SIMULATIONS
The results of numerical simulations are presented in this section in order to investigate the three KPIs mentioned in the Introduction. Only models for the D2Q9 velocity set are discussed in this paper. Extensions of models to high velocity sets and three dimensions exist in the literature. 67,68 It should be noted that extensions of the models to 3D are not trivial. D'Humières et al. 40 extended the MRT collision operator to three dimensions and benchmarked its capabilities using diagonally lid-driven cavity flow. Li et al. 69 devised, with the help of Chapman-Enskog analysis, a non-orthogonal 3D MRT model for multiphase flows. Li et al. 69 reported that the non-orthogonal approach is easier to implement and offers the same stability improvement as an orthogonal MRT model. Li et al.'s 57 forcing method of thermodynamic consistency adjustment was extended to 3D by various authors including Zhang et al., 70 Xu et al., 71 Lycett-Brown and Luo, 72 and Ammar et al. 68 Huang-Wu 44 forcing scheme does not appear to have been extended to three dimensions. Pseudopotential LBM is a diffuse interface method and in order to avoid stability issues, the phases are initialized with a prescribed interface thickness using the following formula for circular shapes 73 : where W is the interface thickness, r 0 is the spherical shape radius, and subscript c stands for center. Flat, vertical interfaces can also be initialized in a diffuse manner:

Thermodynamic consistency
Thermodynamic inconsistency impedes simulation of high density ratios using pseudopotential models. The problem and improvements designed to overcome it are investigated in this subsection. Numerical simulations were carried out in a 200 unit x 30 unit lattice surrounded by periodic boundaries in all four directions using flat interfaces at x = 50 and x = 150. Interface thickness was initially set to 3.5 units, as calculated according to Li and Luo. 57 Equations of state were included in the models using the YS method. Figure 1 illustrates the thermodynamic consistency test. The problem manifests itself in gas densities deviating from the values given by the Maxwell equal-area construction. 74 The Maxwell equal-area construction can be written in the following form 75 : Thermodynamic inconsistency stems from the fact that the mechanical stability condition has a different form than the Maxwell construction. The mechanical stability condition has the following form 76 : Figure 2 illustrates the effect of changing the value of using the Huang-Wu and Li-Luo forcing improvements. Both forcing improvements give results that are significantly better than the unmodified forcing scheme. For a given value of , both methods result in almost identical coexistence densities except for lower reduced temperatures with differences becoming noticeable around T r = 0.6. The differences at lower reduced temperatures are due to terms associated with surface tension modification by the Huang-Wu scheme.
Expanding the interface alleviates thermodynamic inconsistency and reduces spurious velocities. Interface expansion becomes a useful strategy when the reduced temperature is approximately 0.6 or less.
The limits of static high density ratio simulations are explored in Figure 3. It can be seen that the Huang-Wu forcing scheme with expanded interface allows simulations of very high density ratios. It is important to realize that the whole density range may not be suitable for practicable dynamic simulations. MPI with the SRT collision operator is also capable of achieving very high density ratios in static simulations without trial-and-error adjustments of . Thus, the advantages of using several thermodynamically consistent multipseudopotentials become apparent. Viscosity independence is an important asset of MRT models for multiphase simulations. It is desirable to have a model which can simulate high viscosity ratios. Besides stability at low values of kinematic viscosity, equilibrium densities of the phases should not be affected by the kinematic viscosity and, by extension, by v. 77 Shear viscosity is independent of the fluid density for a given temperature and is approximately equal to 78 : where m is the molecular mass, k is the Boltzmann's constant, T is the temperature, and a is the molecular diameter. Hence, changing the kinematic viscosity should not affect the equilibrium densities. Huang-Wu and Li-Luo forcing models give equilibrium densities practically unaffected by the value of the kinematic viscosity. Li and Luo 57 also reported that the stability of multiphase simulations can be increased by setting the kinematic viscosities of the gas and liquid phases to different values (i.e., setting the gas kinematic viscosity to a value multiple times higher than the liquid kinematic viscosity). 57 This can be implemented in single-component simulations by making the value of v dependent on the local value of density, with critical density acting as the dividing point between the gas viscosity and the liquid viscosity. This is permissible because can be a function of space and time. 14 Due to the fact that it is not assumed in the Chapman-Enskog analysis that t = 0 or = 0. 14 Surface tension should also be unaffected by changing the kinematic viscosity and this is illustrated in Figure 4. Surface tension of both schemes is unaffected when the kinematic viscosity is changed. The trendlines, in all the cases, have their intercepts very close to the origin, as predicted by Laplace's law. Surface tension obtained using the Huang-Wu scheme when k 1 = k 2 (−0.1132), for an epsilon value of 1.81, is 1.68 times higher than that obtained using the Li-Luo scheme with = 0. Surface tension in the Huang-Wu scheme is adjusted using the value of 1 -6k 1 . For k 1 = −0.1132, the surface tension coefficient ( c ) is 1.68. The same surface tension value can be obtained using both schemes when the k 1 coefficient in the Huang-Wu scheme is set to 0. Hence, both schemes give the same surface tension, which can be referred to as the unmodified surface tension, when k 1 = 0 in the Huang-Wu scheme and = 0 in the Li-Luo method.

Spurious velocities
This subsection investigates spurious velocities for different models. Simulations to investigate spurious velocities were carried out by initializing a liquid droplet with a diameter of 50 lattice units surrounded by vapor in a 201 units by 201 units lattice. All four boundaries were periodic. Figure 5 illustrates the simulation setup to investigate spurious velocities. Spurious velocities were recorded as the magnitude of physical velocity at a given point: Figure 6 represents the effects of changing viscosity and density ratios on spurious velocities for the Huang-Wu and the Li-Luo forcing schemes. The interface was expanded by setting the value of a in the CS EOS to 0.25, in order to achieve stability and reduce spurious velocities. It can be seen that decreasing the kinematic viscosity has an even greater effect on spurious velocities than increasing the density ratio. Spurious velocities increase by an order of magnitude for an order of magnitude decrease in viscosity. Viscosity varies by an order of magnitude between v = 1.5 (v = 0.3333) and v = 0.6 (v = 0.0333). It varies, again, by an order of magnitude between v = 1.0 (v = 0.1667) and v = 0.55 (v = 0.0167). Whereas, increasing density ratio by an order of magnitude from 100 to 1000 increases spurious velocities approximately 2.5 times. Both forcing modifications result in practically the same spurious velocities, for all the conditions tested, when the surface tension is unmodified by either scheme, that is, when k 1 = 0 and = 0. Figure 7 illustrates the effects of modifying surface tension, by changing the k 1 coefficient, in the Huang-Wu method while keeping constant. Evidently, surface tension needs to be taken into account when investigating spurious velocities, especially when comparing different methods. Increasing surface tension to 0.0118 (by setting k 1 to −0.1163, which leads to c = 1.6792) reduces spurious velocities and reducing surface tension to 0.0043 (by setting k 1 to 0.0667, which leads to c = 0.6) increases spurious velocities in comparison to the unmodified surface tension value of 0.007 (obtained by setting k 1 to 0, i.e., c = 1.0). Figure 8 illustrates the effects of modifying surface tension using the Li-Luo method. The value of is kept constant and is varied to modify surface tension. Reducing surface tension by setting to a positive value reduces spurious velocities. Li and Luo 64 did not mention the possibility of setting to a negative value to increase surface tension. Surface tension is modified by the result of 1 -; therefore setting κ to a negative value should increase surface tension. Setting to −0.4 was found to increase surface tension to 0.0098, as predicted. However, increasing surface tension was found to increase spurious velocities, in contrast to the Huang-Wu method. The changing density ratios for a given reduced temperature at different values of surface tension can also be observed in Figure 7 and in Figure 8. Up to this point, all of the results were obtained using methods directly introducing cubic equations of state into pseudopotential calculations. A piecewise linear EOS is also investigated in this article and it is combined with the Huang-Wu forcing method to allow modification of and adjustment of surface tension. Li and Luo 45 investigated the effects of changing the slopes of pressure in the vapor, unstable and liquid branches. They suggested setting V to 0.64c 2 s , L to c 2 s , and M to −0.04c 2 s . 45 Setting V to a value that is similar to L reduces the variation of vapor density with droplet size. 45 In this article, variation of density ratio with changes in surface tension and relaxation rates is examined. Interface thickness is expanded when the value of M approaches zero. 45 In order to compare the YS and the piecewise linear methods of EOS inclusion, it is necessary to set the interfacial thickness to the same value using both methods. In order to achieve this, the value of M in the piecewise linear method is varied to change the interfacial thickness. The task is complicated by the fact that the slope of the interface is different in the two methods. The YS method has a steeper interfacial slope near to the liquid phase and the piecewise linear EOS method has a steeper interface adjacent to the vapor phase. Figure 9 compares the interfacial region when the interface cut-off is taken to be 76% of the equilibrium densities. In this article, it was decided to compare both methods by setting the value of a in the CS EOS in the YS method to 0. 25  The piecewise linear EOS using the parameters listed in Table 1 has a lower surface tension than the YS method with the CS EOS (a set to 0.25). In order to reach the surface tension value of 0.007, M has to be set to −0.058c 2 s . However, this value of M would result in a thinner interface which would make comparison unfair. Therefore, the value of 1-6 k 1 used in the Huang-Wu method employed with the piecewise linear EOS was increased from 1.0 to 1.73 to obtain the same surface tension as when the YS method was used without changing the interfacial thickness. The value of surface tension was successfully set to 0.007. Figure 10 shows that spurious velocities are lower using the piecewise linear EOS than when the YS method is used. However, the piecewise linear method has problems reaching equilibrium at low viscosities. Problems reaching equilibrium become apparent around v = 0.53. No equilibrium was reached using the piecewise linear EOS at v = 0.53 for T r = 0.5. Also, no equilibrium was reached for a larger droplet (100 l.u. diameter) at v = 0.55 for T r = 0.5.
Wu et al. 79 analyzed the causes of instability of pseudopotential LBMs and proposed several strategies for extending the envelope of stable parameters. Two causes of instability are directly linked to the EOS when the YS method is used. As illustrated by Wu et al., 79 the pseudopotential becomes a complex number at high density values due to (the term inside the square root) becoming negative. A sign function is required to change the value of G from −1 to 1, in order to avoid loss of stability when high density values are encountered. Singularities of the EOS are another cause of instability that requires a limiter function restricting density to a prescribed value. 79 Singularities of the EOS are an inherent deficiency of cubic equations of state. Figure 11 investigates the pseudopotential magnitude at different values of specific volume when the piecewise linear EOS is used. Clearly, the pseudopotential behaves differently when the EOS is introduced using the piecewise linear EOS than when it is done using the YS method. The pseudopotential is asymptotic in the vapor region, that is, below the lower spinodal point 1 . It reaches a constant value in the liquid region above the second spinodal point, that is, 2 . The pseudopotential increases monotonically in the transition region. This behavior of the pseudopotential makes the piecewise linear EOS method suitable for use in simulations.
The distribution of spurious velocity is also of interest. It is well-known that spurious velocity is the highest in the vapor phase and tends to be much lower in the liquid region. Figure 12 illustrates the distributions of spurious velocity and density and compares how they are affected for the YS and piecewise linear methods by changing the surface tension using the Huang-Wu forcing scheme. Spurious velocity is found to be the highest in the vapor region. Spurious velocities are lower in the liquid region when the piecewise linear EOS is used in comparison to when the YS method is used. Density ratio generated by the piecewise linear EOS is found to be almost unaffected by changing the k 1 parameter to vary the surface tension; whereas, density is sensitive to the k 1 parameter when the YS method is employed in the pseudopotential. This phenomenon is caused by high-order error terms introduced by the modifications to the forcing scheme and it is further explained in Section 5.5. The same observation is made when the relaxation rates are varied.

Interparticle interactions
Interparticle interactions are investigated in this subsection. Unlike single-pseudopotential interaction (SPI) models, MPI consists of multiple interactions which are additive and together make up the total interaction force. In order to compare the magnitude of the different interparticle forces, the pseudopotentials should be multiplied by their respective G values.
The following graphs illustrate the interaction strength at the binodal curve starting with Figure 13, which shows how the separate multipseudopotentials generate the total MPI force. The total pseudopotential force generated by the MPI and Yuan-Schaefer (YS) models depends on the interfacial thickness. The concept of the interface thickness parameter can be used as expressed by Khajepor and Chen 60 : The comparative performance of the two schemes varies depending on the chosen EOS parameters. The following graphs illustrate examples of the different behaviors exhibited by the two schemes.

Sharp interface (CS EOS a = 1, b = 4; interface thickness parameter = 4)
The YS model is unsuitable for simulation of high density ratios with a sharp interface. The pseudopotential force becomes undefined at high values of density which are required for simulations of low reduced temperatures. On the other hand, MPI generates correct intermolecular forces in this case, as illustrated in Figure 14. a = 0.01, b = 0

.2; interface thickness parameter = 20)
In order to simulate high density ratios, the interfacial thickness should be expanded when using the YS and MPI models. Practicable interface thickness can be achieved by setting the a and b parameters in the CS EOS to, for example, 0.25 and 4, respectively, or by setting them to 0.01 and 0.2, respectively, as demonstrated in this subsection. When this interfacial F I G U R E 15 Pseudopotential forces generated by MPI and YS methods when the interfacial thickness is set to a practicable thickness [Colour figure can be viewed at wileyonlinelibrary.com] F I G U R E 16 Pseudopotential forces generated by the two models for a thick interfacial thickness when the CS EOS a parameter is set to 0.0025 [Colour figure can be viewed at wileyonlinelibrary.com] thickness is used, the pseudopotential does not become undefined for either model. However, the YS method still exhibits undesired effects at very low values of specific volume, as illustrated in Figure 15.

5.3.3
Thick interface (CS EOS a = 0.0025, b = 0.1; interface thickness parameter = 40) When this set of EOS parameters is used, MPI crashes when density reaches 20 due to negative denominator in the fourth multipseudopotential. On the other hand, the YS model generates the desired intermolecular force as illustrated in Figure 16. a = 0.005, b = 0

.2; interface thickness parameter = 40)
The EOS parameters can be changed whilst keeping the same interfacial thickness. MPI does not generate undefined interparticle forces when the EOS parameters are changed whilst keeping the interface thickness parameter constant as illustrated in Figure 17.

Interactions at a constant temperature across the phase envelope
The interparticle interactions can also be compared across the phase envelope at a constant temperature. In Figure 18, it can be seen that the only method to avoid generation of undefined pseudopotential forces across the phase envelope is the piecewise linear EOS. The other two methods generate undefined pseudopotential forces within the liquid region at densities higher than at the binodal curve due to an inherent inability of cubic equations of state to calculate pressure within that region.

Viscosity ratio between phases
Li and Luo 57 found that stability at low values of liquid kinematic viscosity can be improved by introducing a viscosity ratio between the vapor and liquid phases. This strategy is further investigated in this paper. It appears to resemble the strategy for improving stability in the entropic models. However, viscosity increase in entropic models is highly localized around the critical regions. 24 Whereas, increasing the kinematic viscosity of the vapor phase as suggested by Li and Luo affects the whole vapor region. 57 Figure 19 shows that spurious velocities are significantly reduced by increasing the kinematic viscosity ratio between the vapor and liquid phases from 1 to 20. The sensitivity of spurious velocities to density ratio increase is also almost eliminated within the examined density ratio range. Since a piecewise linear EOS was found to encounter problems reaching equilibrium at low viscosities, it was decided to investigate whether introducing a kinematic viscosity ratio between the phases would alleviate the problem. As depicted in Figure 20, the lowest stable kinematic viscosity for the density ratio range can be improved 10 times from 0.01667 at v = 0.55 to 0.00167 at v = 0.505 by introducing a viscosity ratio equal to 120.

Multiple-relaxation-time relaxation rates adjustment
Simulations were carried out in order to investigate the effects of changing the relaxation rates of the MRT collision operator. Figures 28 and 29 illustrate the fact that changing the relaxation rates has a significant impact on spurious velocities (and subsequently the stability) and thermodynamic consistency. The simulations were initialized at the same reduced temperatures (i.e. T r = 0.8, 0.65, 0.6, 0.5, 0.46), but changing the relaxation rates resulted in different density ratios. This means that the must be adjusted for each set of relaxation rates when the YS method is used. Wu et al. 80 examined this and stated that it is caused by fourth-order terms. Wu et al. 80 suggested setting the anti-symmetric mode relaxation to 1.99 to reduce the effect of changing the relaxation rates on the density ratio. This strategy equates to setting the magic parameter to a very low value, that is, 0.0013 when v = 1.0. Piecewise linear EOS is not significantly affected by changing the relaxation rates as illustrated in Figure 21. Density ratio obtained using the YS method diverges from that obtained using the piecewise linear EOS as the value of the magic parameter is increased. The parameters in the fourth-order equations responsible for the undesired effects are Q 1 and Q 8 . Attention is focused on Q 1 , because Q 8 is non-zero only when surface tension is modified, that is, when k 1 ≠ 0. Even when Q 8 is not equal to zero, its maximum value tends to be lower (between 10 and 50 times lower at the k 1 values tested) than the maximum value of Q 1 and it affects a smaller proportion of the domain than Q 1 does. Figure 22 shows two interesting observations. Firstly, the YS method generates significantly greater values of Q 1 than the piecewise linear EOS. Secondly, changing the magic parameter does not affect the value of Q 1 to a significant extent. Therefore, the change in the density ratio as the magic parameter is varied stems from the terms by which Q 1 is multiplied in the fourth-order equations derived by Wu et al., 80 that is: Similar observations are made regarding the magnitude of the thermodynamic consistency term in the Li-Luo method. The magnitude of the correction term generated by the YS method is significantly greater than generated by the piecewise linear EOS and varying the magic parameter does not have a significant influence on the magnitude of the correction term.
On the other hand, changing the surface tension by adjusting k 1 has an effect on the magnitude of Q 1 as illustrated in Figure 23 and in Figure 24. Hence, the influence of surface tension on the density ratio and on the spurious velocities. The value of Q 1 increases as k 1 is set to greater values and this is the reason why the Huang-Wu forcing scheme has negative stability effects for low values of surface tension, that is, when the result of 1 -6k 1 is less than 1.0. Figure 25 illustrates the effect of the maximum value of Q 1 on the density ratio. Density ratio does not vary significantly when k 1 is adjusted and the magic parameter is changed when the piecewise linear EOS is employed due to low values of Q 1 . The strategy proposed by Wu et al. 80 of setting the antisymmetric relaxation rate to 1.99 results in a significant density ratio change when surface tension is adjusted as illustrated in Figure 25.
Both methods generate different maximum values of Q 1 due to different maximum values of the intermolecular force. The slopes of density and the pseudopotential generated by the YS and piecewise linear methods are different in the interfacial region as illustrated in Figure 26. This means that the values of the intermolecular force (F x and F y in a 2D case), have different distributions with different maxima and this can be seen in Figure 27. The total value of the intermolecular force is almost identical for both methods with approximately 0.3% difference between them, as shown in Table 2.  Figure 28 illustrates that setting the magic parameter to 1/12 gives the best results for multiphase simulations, especially for high density ratios. Setting the magic parameter to an even lower value can result in lower spurious velocities at low density ratios. The break-even density ratio for magic parameters equal to 1/12 and 0.031 is approximately 50. Third-order spatial errors are canceled when magic parameter is set to 1/12. 14,44 Figure 29 illustrates that increasing the bulk viscosity to multiple times the kinematic viscosity increases the stability of simulations. Stability at the lowest reduced temperature tested, that is, 0.46, was achieved when bulk viscosity was higher than the kinematic viscosity. Increasing bulk viscosity was found to slightly increase the maximum spurious velocities, in the case of a static droplet suspended in a gravity-free domain. However, in dynamic cases increasing the value of bulk viscosity can significantly reduce spurious waves generated in the vapor phase as illustrated in Section 5.6.

Dynamic simulations of a droplet splashing on a thin liquid film
The findings of dynamic simulations are presented in this subsection. The well-studied problem of a droplet splashing on a thin liquid film is used as an example of the dynamic capabilities of the schemes. The investigated parameters include the kinematic viscosity, bulk viscosity, reduced temperature and the viscosity ratio between the phases. The setup consisted of a 600 units by 250 units lattice with periodic boundaries on the left and right-hand sides and a simple bounce-back boundary on the bottom wall. A 25-units deep liquid film was resting on the bottom wall and a droplet with a diameter of 100 units was initialized with a prescribed impact velocity. Figure 30 illustrates formation of satellite droplets during the impact of a droplet on a thin liquid film at Reynolds number (Re) equal to 1000 and Weber number (We) equal to 110. The dimensionless numbers are expressed as follows: Figure 31 investigates the effect of the chosen models and parameters on the achievable Reynolds number. Clearly, increasing the value of bulk viscosity to 5 or 10 times the value of kinematic viscosity significantly improves the stability of dynamic simulations. Bulk viscosity has a greater influence on the stability than the reduced temperature. Reynolds number of 300 was achieved even at high density ratios when bulk viscosity was 10 times the value of kinematic viscosity. Instability in the problem of a splashing droplet is triggered by high velocity regions in the vapor phase as illustrated in Figure 32. The formation of these high velocity regions can be inhibited by increasing the value of bulk viscosity as illustrated in Figure 33. Hence, the significant improvement of stability when bulk viscosity is increased.
The piecewise linear EOS offers improved performance allowing to increase the Reynolds number by approximately 50% when the bulk viscosity is the same as or five times the value of kinematic viscosity. MPI performs only marginally better than the YS method when interface is set to a practicable thickness as illustrated in Figure 31. The comparative performance of MPI in relation to the YS method improves when interfacial thickness is set to be sharper, as illustrated in Figure 34. These findings corroborate the nature of interparticle forces obtained using both methods depicted in Figure 14.
Introducing a viscosity ratio between the phases was found to be the most effective strategy for achieving high Reynolds numbers. Evidently, this strategy has a significant effect on the stability and allows to achieve a Reynolds number in the region of 1000 even when the bulk viscosity is equal to kinematic viscosity. However, the benefits of increasing the kinematic viscosity ratio between the phases eventually taper off and no performance improvement has been observed when the kinematic viscosity ratio was increased from 80 to 250. This tapering off of the benefits of further increasing the kinematic viscosity ratio between the phases is illustrated in Figure 35.
At a Reynolds number of 1000, the Li-Luo forcing scheme allows to achieve a Weber number value of 360 while the Huang-Wu model allows to achieve a value of 320 for the Weber number, when both methods are combined with the piecewise linear EOS. This confirms the findings that the Li-Luo method of adjusting surface tension is superior to the Huang-Wu method for low values of surface tension.

CONCLUSIONS
In this article, a number of pseudopotential models are investigated under static and dynamic conditions. From the results of the numerical simulations, it can be concluded that the MRT collision operator should be combined with the Huang-Wu 44 forcing scheme modification for high values of surface tension and with the Li-Luo 57,64 scheme for low values of surface tension. Both forcing scheme modifications perform in an equivalent manner when surface tension is unmodified, that is, k 1 = 0 for the Huang-Wu forcing method and = 0 for the Li-Luo method. Understanding the ways in which different models generate interparticle forces can be beneficial for ensuring desired simulation outcomes. Methods which directly incorporate cubic equations of state into pseudopotentials models can generate undefined forces even when the interface is smooth.
Piecewise linear EOS was found to perform better than the YS method of EOS inclusion. It not only avoids two sources of instability that can be encountered in the YS method, but the piecewise linear EOS method also helps to decouple parameters. The pseudopotential does not become a complex number at high values of density, which is a problem for the YS method. Instead, the pseudopotential reaches a constant maximum value in the liquid region denoted by the upper spinodal point. The piecewise linear EOS also does not have singularities which occur in the YS method. Hence, two of the stability strategies proposed by Wu et al. 79 become unnecessary for the piecewise linear EOS.
Combining the Huang-Wu and Li-Luo surface tension adjustment methods with the YS method causes significant changes in the density ratio, which may lead to a requirement for adjustment. Density ratio is affected by surface tension adjustment to a significantly lesser extent when the piecewise linear EOS is employed. Better performance is caused by lower values of high-order errors in the interfacial region.
At low values of kinematic viscosity, the piecewise linear method encountered problems achieving equilibrium in static tests. Difficulties achieving equilibrium at low values of kinematic viscosity using the piecewise linear EOS method can be overcome by increasing the value of vapor kinematic viscosity to multiple times the liquid kinematic viscosity. Imposing a kinematic viscosity ratio between the phases almost eliminates the sensitivity of spurious velocities to the density ratio.
The effects of changing the relaxation rates in the MRT collision operator on multiphase simulations are examined. Changing the relaxation rates influences thermodynamic consistency when the YS method is employed with the Huang-Wu or Li-Luo surface tension adjustment methods; therefore, may need to be calibrated for each set of relaxation rates. The effect of relaxation rates on thermodynamic consistency can be significantly alleviated by using the piecewise linear EOS method.
The factors affecting the stability of dynamic simulations are examined. It can be found that they have the following hierarchy from the greatest impact to the least: kinematic viscosity ratio between the phases, bulk viscosity, method of EOS inclusion, and reduced temperature/density ratio. See Tables A1 and A2. TA B L E A1 Where F α are the total forces which can include, depending on the case being simulated, fluid-fluid interactions, fluid-solid interactions, and external forces (e.g., gravity) and |F| 2 = (F int 2 x + F int 2 y ) (F α int stands for intermolecular, i.e., fluid-fluid interactions)