Cascaded lattice Boltzmann method with improved forcing scheme for large-density-ratio multiphase flow at high Reynolds and Weber numbers

A recently developed forcing scheme has allowed the pseudopotential multiphase lattice Boltzmann method to correctly reproduce coexistence curves, while expanding its range to lower surface tensions and arbitrarily high density ratios [Lycett-Brown and Luo, Phys. Rev. E 91, 023305 (2015)]. Here, a third-order Chapman-Enskog analysis is used to extend this result from the single-relaxation-time collision operator, to a multiple-relaxationtime cascaded collision operator, whose additional relaxation rates allow a significant increase in stability. Numerical results confirm that the proposed scheme enables almost independent control of density ratio, surface tension, interface width, viscosity, and the additional relaxation rates of the cascaded collision operator. This allows simulation of large density ratio flows at simultaneously high Reynolds and Weber numbers, which is demonstrated through binary collisions of water droplets in air (with density ratio up to 1000, Reynolds number 6200 and Weber number 440). This model represents a significant improvement in multiphase flow simulation by the pseudopotential lattice Boltzmann method in which real-world parameters are finally achievable.


I. INTRODUCTION
The lattice Boltzmann method (LBM) is a mesoscopic simulation approach that provides an attractive alternative to traditional computational fluid dynamics.Instead of solving the continuum macroscopic equations the LBM solves a discretized Boltzmann equation, designed to recover the Navier-Stokes equations in the macroscopic limit.Its kinetic origin makes the LBM amenable to incorporation of microscopic and mesoscopic physics.The large variety of problems to which the method has been successfully applied includes, among others, turbulent flows, microflows, flows through porous media, magnetohydrodynamics, multicomponent systems, and thermal multiphase flow (see, for example, Refs.[1][2][3][4] and references therein).In principle, multiphase behavior including interfacial phenomena can be naturally handled by the LBM by the inclusion of a forcing term in the discretised Boltzmann equation.However, the long-standing debate over the form of such a forcing term that will correctly recover a force in the Navier-Stokes equation has only recently been resolved.Lycett-Brown and Luo [5] proposed a new forcing scheme that correctly reproduces the multiphase coexistence curve at arbitrarily high-density ratios and allows for accurate control over surface tension.Here the method is generalized from a single-relaxation-time collision operator to a collision operator with multiple relaxation times.
Since the early stages of the development of the LBM, its potential to be easily extended to multiphase flow has been recognized.This has resulted in the proposal of a variety of different methods for introducing multiphase behavior.These include the free-energy models [6][7][8], those based on the kinetic theory of dense fluids [9][10][11], and the interaction potential models [12][13][14].These methods can all be formulated in the form of an additional forcing term to the collision operator.In the present paper the interaction potential model of Shan and Chen [12] is considered.The limitations of the original formulation of the Shan-Chen model to low liquid to gas density ratios and a small range of viscosities is well documented, as is the problem of the formation of spurious velocities around curved interfaces.A number of improvements to the model have therefore been proposed which can be divided into two categories: those that modify the force calculation, such as increasing the order of isotropy [15] or modifying the equation of state [16], and those that vary the incorporation of the force term into the equilibrium distribution functions.It is the latter category that is considered here.
Along with the original forcing term by Shan and Chen [12], a number of other forcing schemes have been proposed.They include the method of explicit derivatives, as used in the multiphase schemes derived from the kinetic theory of dense fluids [9][10][11], the method of Guo et al. [17], which takes into account discrete lattice effects and the exact difference method (EDM) [18].Results show that the EDM is stable over a larger range of density ratios than the method of Guo et al. [17], and gives smaller errors for high density ratio multiphase systems [18,19], when using the Carnahan-Starling equation of state.Conversely, the method of Guo et al. [17] performs better than the EDM for the exponential form of the equation of state proposed by Shan and Chen.This had previously resulted in some confusion over which forcing scheme is the most appropriate to use.Additionally, analysis has shown that while the method of Guo et al. [17] recovers the Navier-Stokes equations correctly, the EDM introduces an error into the pressure tensor, proportional to the square of the forcing term [17,20].However, the analysis used is based on a second-order expansion of the LBM, which has been shown to be insufficient to identify all relevant error terms when considering forces with large gradients, as are present at interfaces in multiphase schemes [21].Lycett-Brown and Luo [5] therefore extended the third-order truncation error analysis of Holdych et al. [22] to include forcing terms, in order to cor-rectly identify the errors in the pressure tensor.Consequently, a forcing scheme was derived accounting for both these higher-order errors, and the inherent lack of thermodynamic consistency in the Shan-Chen model.As the error terms under discussion affect the coexistence curve it is possible for errors in the pressure tensor to counteract this lack of thermodynamical consistency.The higher-order analysis also revealed how this effect can lead to some forcing schemes working better under certain circumstances, despite introducing larger errors.The outstanding debate over the correct form of the forcing term was therefore resolved with the introduction of the new forcing term based on higher-order analysis.Following Li and Luo [23], Lycett-Brown and Luo [5] also added a term to the force calculation which, when combined with the new forcing scheme, allowed for accurate control of the surface tension, independently of other model parameters.Independent control of density ratio, viscosity, surface tension, and the numerical width of the diffusive interface was shown, resolving one of the longstanding criticisms of the pseudopotential multiphase LBM.Significantly, independent control of interface width allows for systematic reduction in spurious velocities as interface width is increased.
In its simplest form, the LBM uses a single relaxation time Bhatnagar-Gross-Krook (BGK) collision operator.It has been extensively shown in the single-phase case that a significant improvement in attainable Reynolds number can be obtained through the use of a multiple-relaxationtime collision operator [24][25][26][27].Instead of relaxing particle distribution functions toward their equilibrium distribution functions, a transformation is made into moment space, where individual moments can be relaxed independently, and the relaxation rates of higher-order moments can be used to increase stability.Attempts have been made to apply this to the various multiphase schemes, leading in some cases to an increase in attainable Reynolds number [28][29][30][31].A decrease in spurious velocities has also been observed.Recently, Geier et al. [32] proposed the cascaded LBM, in which the collision operator relaxes moments in a comoving reference frame, and showed improvements over results obtained using the MRT collision operator.Lycett-Brown et al. [19,33] provided an alternative derivation of the cascaded LBM and extended it to a pseudopotential multiphase scheme, resulting in increased stability at low viscosity, and the reduction of spurious velocities by tuning the relaxation rates of the higher-order moments.The EDM was used and the interface width could not be adjusted independently of model parameters.The tuning of the relaxation rates therefore represented a compromise between reduction of spurious velocities and increased stability.Here the new forcing scheme is extended to account for the additional relaxation rates of the cascaded collision operator.A third-order Chapman-Enskog expansion is used, which is shown to give the same result as the third-order truncation error analysis [5] based on the method of Holdych et al. [22], in the limit of all relaxation parameters being equal.
The cascaded LBM with a generic forcing term and the derivation of the pressure tensor in the pseudopotential model are summarized in Sec.II.The third-order Chapman-Enskog analysis of the cascaded LBM, including a generic forcing term, is then given in Sec.III.Results for the reproduction of coexistence curves in one dimension are given in Sec.IV A, along with a study of the effect of varying model parameters and the additional relaxation rates on density ratio.Results for the two-dimensional case are then given in Secs.IV B, IV C, and IV D, including surface tension measurements, and the effect of varying surface tension and interface width on the spurious velocities found around curved phase boundaries.Section V gives results for the model applied to the problem of three-dimensional binary droplet collision at high density and viscosity ratios, and high Weber and Reynolds numbers.Finally, a summary and outlook are given in Sec.VI.

II. THE CASCADED LATTICE BOLTZMANN METHOD
The LBM with a general forcing term, F i , can be written as where f i are particle distribution functions for particles with velocities v i , x, and t are space and time coordinates, and for a single-relaxation-time collision operator, collisions are written as where ω = 1/τ is the relaxation rate and τ is the relaxation time.The number of discrete velocities, Q, depends on the lattice, where subscript i is used interchangeably throughout with subscript (i,j ) in two dimensions and (i,j,k) in three dimensions.In the former case, i = 1, . . .,Q, and the latter case refers to lattice velocities v (i,j,k) = (v (i) ,v (j ) ,v (k) ).
For the MRT and cascaded collision operators, particle distribution functions are expressed in terms of velocity moments of the distribution functions, with different relaxation rates associated with different moments.The collisions, f * i , for the cascaded collision operator are given in Appendix A. The density, ρ, and macroscopic fluid velocity, u, are found as the zeroth and first velocity moments of the distribution functions: ( Higher-order moments are defined as for the raw moments used in the MRT method, and for the central moments used in the cascaded method.The second, third, and fourth moments are written as The form of F i depends on the method, following Lycett-Brown and Luo [5] here it is defined as where F is the force acting on the fluid and γ is dependent on the method and discussed in detail in the following.c s is the speed of sound and T 0 = c 2 s is used throughout. is introduced for control of surface tension, and will also be defined in the following.For all forms of the forcing term the zeroth and first moments are With the inclusion of the term, the second moment is and we denote the general third moment of F i as .
The derivation of the Navier-Stokes equation from the LBM is deferred to the following section.For now we assume F is correctly introduced into the Navier-Stokes equation, and consider the resulting pressure tensor when F has the form of a pseudopotential force.For now the effect of is also neglected.In the Shan-Chen model the force is given by [12] where ψ is the density dependent interaction potential, G controls the strength of the interaction, and w(|v i | 2 ) are weights.The number of discrete velocities, N , used in the force calculation does not have to be equal to the number of lattice velocities.In general, weights are also different from those in the equilibrium distribution functions.The Taylor series expansion of ψ(x + v i ) gives [10] F(x) = −Gc 2 s ψ∇ψ − λGc 2 s ψ∇( ψ), (10) where λ depends on N .Using this equation and the definition the pressure tensor, P, is given by where I is the identity matrix.λ = 1/6 for a two-dimensional lattice with nine velocities (N = 9).An alternative derivation of the pressure tensor, originally proposed by Shan and Chen [13] and more recently extended by Shan [34] results in different coefficients in place of each λ in Eq. ( 12).The derivation is based on the discrete form of the integral of Eq. (11).However, as discussed in Lycett-Brown and Luo [5], these coefficients are not correct, therefore this result is not discussed further here.
The mechanical stability condition is now considered, following Shan [34].From Eq. ( 12), the pressure normal to a flat interface with gradients only in the x direction is given by where a and b are here given by a = −λ/2 and b = λ.Setting this equal to the static pressure in the bulk, p 0 , and integrating with respect to ρ (see, for example, Li et al. [35]) gives ∂ρ ∂x where subscripts l and g refer to the liquid and gas phases, respectively, and As ∂ρ/∂x = 0 in the bulk then the liquid and gas densities must satisfy From thermodynamic theory the Maxwell construction gives [13] therefore for thermodynamic consistency it is required that [36] 1 It should be noted that to be truly thermodynamically consistent, the stress tensor should match the Kortweg stress tensor.As this is not possible in the context of the Shan-Chen model, the above requirement, which allows the correct recovery of the liquid and gas densities, is termed thermodynamic consistency both in the following and in the literature.Surface tension is defined as [37] which from Eq. ( 12) is given by These expressions for the mechanical stability condition and surface tension assume that F is introduced into the Navier-Stokes equation without any error.As shown in Lycett-Brown and Luo [5] and discussed further in the following section, this is not the case.To correctly identify the errors the second-order Chapman-Enskog expansion often used in LBM analysis is insufficient, as the second-order derivatives in the pressure tensor are of higher order than the second-order in the Navier-Stokes equations [21].However, taking the expansion to higher order leads to Burnett level equations.As this is not the goal, only terms in F are needed from a third-order analysis, because F and its derivatives are large around phase boundaries.Lycett-Brown and Luo [5] derived the error terms using a thirdorder truncation error analysis following the work of Holdych et al. [22].This analysis applied only to the single-relaxationtime BGK collision operator.Here we proceed with a thirdorder Chapman-Enskog expansion of the cascaded LBM.This is shown to be equivalent to the previous third-order truncation error in the limit of all relaxation rates being equal.

III. THIRD-ORDER CHAPMAN-ENSKOG EXPANSION
The third-order Chapman-Enskog expansion of the cascaded LBM with a general forcing term is now given.The left-hand side of Eq. ( 1) can be expanded by a Taylor series expansion as where The distribution function and time derivative are expanded in powers of ε: where f 0 i = f eq i , and for all n 1.Following Buick and Greated [38], the forcing term is With the right-hand side of Eq. ( 1) expressed in terms of moments, expanding these moments individually in terms of ε is equivalent to expanding populations in terms of ε, Equation ( 1) can now be expanded to third order to give where and the notation f n i (M) is used to indicate the nth order term of the part of the distribution function, which is a function of moment M.
Terms of order ε 0 , ε 1 , and ε 2 are then equated.The zeroth to third velocity moments of the resulting zeroth-order equations, the zeroth, first, and second velocity moments of the first-order equation, and the zeroth and first velocity moments of the second-order equation are then derived.The resulting equations are then combined to derive equations for the conservation of mass and momentum.Full details of the derivation are given in Appendix B. The conservation of mass is given by the familiar form where as usual the fluid velocity, u, has been modified as The conservation of momentum gives where Tr indicates the trace.The deviatoric stress tensor is given in the usual form as where ν and ν b are the shear and bulk kinematic viscosities, given by The last three terms in Eq. (31) are errors in the Navier-Stokes equation.These errors are general and will apply to any 053313-4 force.In the case of a small external force they will have little effect, but they must be taken into account for any force with large gradients.This includes all multiphase methods that can be arranged into the form of Eq. ( 1).The last term in Eq. ( 31) is a new error term associated with the cascaded collision operator.In the BGK limit, τ b = τ , this term is zero and the error in Eq. ( 31) becomes equal to that derived in Lycett-Brown and Luo [5] for the single-relaxation-time collision operator.It should be noted that only the additional relaxation rate related to the bulk viscosity appears in this error term, the additional relaxation rates for the higher-order moments do not introduce further error at this order.Varying the additional relaxation rates of the higher-order moments in the cascaded collision operator to enhance stability should therefore have no effect on model parameters such as density ratio and surface tension.The additional error term introduced by varying bulk viscosity can be accounted for, as discussed below.
While the error terms in Eq. ( 31) are general, here we consider the case where the force is given specifically by Eq. (9).With this expression for the force, the following relationships can be derived: Following Lycett-Brown and Luo [5] we define to be where κ is a parameter for controlling surface tension.With this expression for , and using Eqs.( 35) and (36), the pressure tensor is given by The last term in this equation is the result of the additional relaxation rate related to bulk viscosity in the cascaded collision operator.Again, it should be noted that in the BGK limit, τ b = τ , this reduces to the equation for the pressure tensor derived in Lycett-Brown and Luo [5].As shown therein, it can clearly be seen that the terms from the third-order expansion must be taken into account, and that the effect of the new term should also be included.We first consider the impact of these error terms on the mechanical stability condition, before discussing their effects on surface tension.The pressure normal to a flat interface with gradients only in the x direction is now given by With the error terms included, the coefficients a and b in Eq. ( 13), corresponding to the first and second terms in square brackets in Eq. ( 39), have been changed.Consequently, is now given by where c 2 s = 1/3 and λ = 1/6 have been used.From Eq. ( 16) it can be seen that these errors affect the mechanical stability condition.However, with the correct choice of γ in Eq. ( 6) these errors can be eliminated.Rearranging Eq. ( 40) for γ gives where 0 is a free parameter that should be chosen, dependent on the equation of state, to correct for the thermodynamic inconsistency of the psuedopotential model.For example, solving Eq. ( 18) for thermodynamic consistency with = 0 gives This form for the pseudopotential is common in the literature, and its coexistence curve could be correctly reproduced by choosing 0 = 0 in Eq. ( 41), as discussed in Lycett-Brown and Luo [5].More realistic equations of state can be used by choosing ψ to be [16] and setting p 0 to the desired equation of state.Here G no longer controls the interaction strength and is set to G = −1.While it is not possible to get exact thermodynamic consistency with this form of the pseudopotential, very good approximations can be made, depending on the chosen equation of state.Solving Eq. ( 18) for a general value of and setting this equal to the pseudopotential in Eq. ( 43) gives for ρ = ρ g and ρ = ρ l .While in general this cannot be solved exactly, it can be seen that for ρc 2 s p 0 or p 0 ∝ ρ that thermodynamic consistency can be approximated by setting 0 = 2.In the results we consider the following three equations of state: (1) The Carnahan-Starling equation of state, where T is the effective temperature, and R, α, and β are constants, discussed further below; (2) An approximation about the critical point for a van der Waals equation of state [37,39], where ρ sat g and ρ sat l are the gas and liquid densities at saturation respectively, and A is a constant, also discussed further below; and (3) a piecewise linear equation of state [40], where θ g = (∂p 0 /∂ρ) g and θ l = (∂p 0 /∂ρ) l are the gradients of p 0 in the gas phase and liquid phase equal to the speeds of sound squared in each phase, and θ m = (∂p 0 /∂ρ) m is the gradient of p 0 in the mechanically unstable region of the equation of state.ρ 1 and ρ 2 are the spinodal points, which can be found by simultaneously solving the equations for mechanical and chemical equilibrium [40]: In the following results, for the Carnahan-Starling equation of state an approximation for 0 is found by an iterative procedure, the value tending toward 2 for high-density ratios.For the other two equations of state 0 = 2 is used.
The results for the coexistence curves obtained with the forcing scheme described here are given in Sec.IV.The new terms identified in the pressure tensor will also affect the surface tension.We have introduced an additional term into the forcing term, Eq. ( 6), to allow control of surface tension; however, this control will require accounting for the additional error terms.With the pressure tensor given by Eq. ( 38), and γ given by Eq. ( 41), the surface tension can be found using Eq. ( 19) to be where has been used, and the surface tension coefficient is given by This reduces to the expression given in Lycett-Brown and Luo [5] in the LBGK limit of τ b = τ .Surface tension is proportional to σ c , which can now be controlled by varying κ.The κ dependence in the expression for γ should allow surface tension to be varied without affecting the liquid and gas densities.Results for this are given in the following section.
The surface tension is also proportional to the reciprocal of the width of the diffusive interface and it is important to be able to vary this width without changing other properties of the model.
A relationship for the width of the interface can be found following Jacqmin [41].From the equation for the pressure normal to a flat interface with gradients only in the x direction, Eq. ( 39), a relationship for the interface width coefficient, W c , can be found as where again c 2 s = 1/3 and λ = 1/6 have been used.The width of the diffusive interface, in lattice units, W , is approximately proportional to W c .(As with all diffusive interface methods it should be noted that the numerical interface widths are orders of magnitude larger than the physical interfaces.)This is only an approximation due to the inexact thermodynamic consistency but, as results show, allows control of the interface width with very little error.W is also inversely proportional to the square root of a coefficient in the equation of state.This coefficient is A in Eq. ( 46), and α for the Carnahan-Starling equation of state, Eq. ( 45), with R = α.For the linear equation of state, Eq. (47), we can take this coefficient as θ m .This is only approximate if θ g and θ l are not kept proportional to θ m ; however, it is desirable to maintain independent control of the speeds of sound, and from results the approximation works well at high-density ratios.In the upcoming results the interface width, W , is measured by fitting the following curve to the density distribution across an interface: Although it should be noted that the interface does not exactly follow this function, it is a good approximation.By varying κ and the relevant parameter for each equation of state (and taking into account the required variation of 0 with varying κ for the Carnahan-Starling equation of state), it is possible to vary both surface tension and interface width independently, without affecting the liquid and gas densities (and without affecting their sound speeds, in the case of the linear equation of state).As shown in Lycett-Brown and Luo [5], this result removes the longstanding criticism of the interaction potential multiphase LBM that such independent variation of model parameters was not possible.Here, along with independent variation of shear viscosity, the extension of the forcing term for the cascaded collision operator should allow variation of bulk viscosity and the relaxation rates of the higher-order moments, independently of all these model parameters.

IV. RESULTS
In the following we present results for the new forcing scheme for the multiphase cascaded LBM described above.Lycett-Brown and Luo [5] showed that for the LBGK case the forcing term significantly improves the reproduction of the coexistence curve.For completeness we first reproduce this result for different equations of state, before presenting results for the additional correction derived for the cascaded LBM.The effects of varying the surface tension parameter and interface width, in combination with additional relaxation rates of the cascaded LBM, are then discussed in one dimension.Droplets in two dimensions are then considered.Lycett-Brown and Luo [5] showed that for the LBGK case significant reductions in surface tension could be achieved, and that independent control of the width of the diffusive interface allowed for systematic control of spurious velocities.Here further results are presented for measured surface tension, spurious velocities, and droplet isotropy, for varying each of the parameters of the model, including the additional relaxation rates of the cascaded collision operator.

A. Coexistence curve
All existing forcing schemes can be rearranged into the form of Eq. ( 1), therefore γ can then be derived.It is well known that the original forcing scheme in the Shan-Chen method gave liquid and gas densities dependent on τ , this is due to a τ 2 term appearing in γ , and it is not discussed further here.The exact difference method has this cancels the τ dependence in Eq. ( 40) (the Shan-Chen method coincides with the EDM at τ = 1), and the method of Guo et al. [17] has this removes both the τ and ψ 2 /ρ dependencies in .
As an example, Fig. 1 shows the vapor branch of the coexistence curves obtained with the method of Guo et al. [17] and the EDM, for the equation of state given by Eq. ( 46), compared with the gas density input into the equation of state, ρ sat g .Here, and in the following, gas density is measured from simulation of a flat interface on a double periodic 450 × 5 grid, initialized with half of the domain as the predicted gas density and half as the predicted liquid density and then run to equilibrium.This domain size was found to be large enough that the measured liquid and gas densities far from the interface had converged even for the largest interface widths (i.e., further increasing the domain size would not change the densities measured at the midpoint between the interfaces).Unless otherwise stated, ω = 1 is used, as very little variation in the error in gas density with ω was seen (for example, at T = 0.05 in Eq. ( 45), giving a predicted density ratio of about 400, the error in gas density with the present method for ω = 1 is 7.38% and 7.46% for ω = 1.8, this can also be seen in Fig. 6).
Both the EDM and the method of Guo show significant errors in the gas density, with the EDM performing slightly better than the Guo method.Both schemes are shown using A = 1/10, and increasing the interface widths by decreasing A to 1/40 was found to give worse results in both cases.This agrees with the results for the Carnahan-Starling equation of state given in Lycett-Brown and Luo [5].The results for the present method, using γ given by Eq. ( 41) and 0 = 2 are shown in Fig. 2. Three values of A are given, A = 1/10, A = 1/40, and A = 1/160, giving increasingly wider interfaces.It can be seen that the vapor branch of the coexistence curve using the present method is in much better agreement with the expected gas density and that the error in the gas density reduces systematically with increasing interface width.As a further example, Fig. 3 shows the results for the vapour branch of the coexistence curve for the equation of state given by Eq. (47).Results are shown for θ m = −1/64, which gives interface widths of approximately four lattice units (relatively independent of density ratio; for the present method we found W = 4.2 at ρ sat g = 0.1 and W = 3.9 at ρ sat g = 0.0001).Decreasing θ m increases interface width, and results are also shown for θ m = −1/256 (giving an interface width of approximately 8 lattice units).The method of Guo was unstable for the smaller interface width, and only stable at low density ratios for the larger width.The EDM gave better results than those of the Guo method, with the present method again giving the lowest errors.Figure 4 shows the numerical convergence in terms of the error reduction in the gas density with increasing interface width for the present method, for density ratios from 10 to 10 000.Again all simulations were FIG. 2. Simulation results for gas densities in one dimension, for the equation of state given by Eq. ( 46) with ρ sat l = 1.The present method (in the LBGK limit) with increasing interface width (decreasing A) using A = 1/10 (triangles), A = 1/40 (crosses), and A = 1/160 (circles), is shown compared with the expected density, ρ sat g (black line).Significant improvement over the existing methods (shown in Fig. 1) is observed.
run in a periodic domain of size 450 × 5, this is large enough that even for the largest interface width the midpoint between the two interfaces gives the same density as would be measured at "infinity."Therefore, only resolution in the interface is increased with increasing W , not the whole domain resolution.Convergence is seen to be of second order, a vital improvement over the previous methods whose results did not converge with increasing interface width.Further data for the improvements of the present method in the LBGK limit, over the method of Guo and the EDM are given in Lycett-Brown and Luo [5], here we will focus on the additional correction derived here for the cascaded LBM.
Figure 5 shows the results for the Carnahan-Starling equation of state for the method of Lycett-Brown and Luo [5],  which does not include the additional correction to the pressure tensor derived here for the additional relaxation rates of the cascaded collision operator [the last term of Eq. ( 38  Maxwell construction, the gas densities vary significantly for ω b = ω.This is expected from the pressure tensor derived here, which clearly shows dependence on ω b .To correct for the additional term in the pressure tensor, γ should be defined as in Eq. ( 41). Figure 6 shows the same results as Fig. 5, but with this additional correction included.For all four cases the gas densities are found to be in very close agreement with the Maxwell construction, even at very high density ratios.They are also in close agreement with each other, suggesting that the present method correctly accounts for the additional higher-order errors for the general case, where ω b = ω.It should also be noted that the additional correction to the pressure tensor is not dependent on the relaxation rates of the higher-order moments, ω 3 and ω 4 .We found that varying these relaxation rates resulted in variations in the gas density that were significantly less than the error between the gas density and the Maxwell construction.We therefore will not discuss these relaxation rates further and set ω 3 = ω 4 = 1 for maximum stability in the following results, unless otherwise stated.
The width of the diffusive interface is predicted to depend on ω b (and again to be independent of ω 3 and ω 4 ), in addition to κ. Figure 7 shows results for interface widths measured by fitting Eq. ( 53) to steady-state density profiles.The results are for the Carnahan-Starling equation of state, Eq. ( 45), with α = 1/32 and T = 0.0455 giving a density ratio between the liquid and gas phases of approximately 1000, and ω = 1.The interface width for κ = 0 and ω b = 1 was first measured to be W = 9.46.From this result Eq. (52) was then used to predict interface widths for varying ω b at different κ.These predictions are shown as solid lines on Fig. 7. Also shown are measured interface widths for varying ω b for κ = 0, 1, 2, and 3.All measured values are found to be in very good agreement with the theoretical prediction, lending further evidence that the newly derived correction to the pressure tensor is correct, and that Eq. ( 52) provides a very good approximation for the interface width coefficient W c .The equation used to measure the interface width, Eq. ( 53), is the solution to the interface profile for the equation of state given by Eq. ( 46), near to the critical point.This is used as a further test of the thermodynamic consistency of the present model.Figure 8 shows the profiles across a flat interface for the present method, the EDM and the method of Guo, compared with Eq. (53), for ρ sat l = 1,ρ sat g = 0.2 and A = 1/10.For this a periodic domain of length 200 lattice units was initialized with liquid between x = 50 and x = 150 and gas elsewhere, with the predicted hyperbolic tangent interface profile, and run to equilibrium.The present method clearly provides the best match for the interface profile.
Figures 9 and 10 show results for error in measured gas density for varying the surface tension parameter, σ c over a large range, for the Carnahan-Starling equation of state, Eq. ( 45), and the linear equation of state, Eq. (47), respectively.FIG. 8. Density profiles across a flat interface for the present method (crosses), the EDM (circles), and the method of Guo (triangles), for the equation of state given by Eq. ( 46) with ρ sat l = 1,ρ sat g = 0.2 and A = 1/10, compared with the theoretical solution, Eq. (53).Both are for a density ratio of approximately 1000, using T = 0.0455 in the former case and ρ sat g = 0.001 in the latter.For the linear equation of state θ l = 1/3 and θ g = 1/6.Using Eq. ( 51) for the surface tension coefficient and Eq. ( 52 10.Error in gas density (compared with ρ sat g ) with varying surface tension parameter, σ c , for the present forcing scheme.The results are for the linear equation of state, Eq. (47), with ρ sat g giving an expected density ratio of 1000, and θ l = 1/3 and θ g = 1/6.σ 0 is σ c with κ = 0. Interface width is kept approximately constant at W = 10.11.Measured surface tensions, σ m , against theoretical surface tension σ th (both normalized by surface tension at κ = 0,σ 0 ).Results are for the linear equation of state, Eq. ( 47), with ρ sat g = 0.001 (density ratio of 1000), and θ l = 1/3 and θ g = 1/6.Interface width is kept approximately constant at W = 11.Surface tension measurements are made using the Laplace law.Results for ω = ω b = 1 (crosses), ω = 1 and ω b = 1.6 (circles), and ω = 1.8 and ω b = 1.2 (diamonds) are almost indistinguishable, and in good agreement with the theoretical prediction.slightly in most cases as the surface tension parameter is reduced (measurements of surface tension with varying surface tension parameter are given in the following section).The interface widths also show very little variation, from 9.34 to 9.61 lattice units in the first case, and from 10.61 to 10.73 in the second.Even for this high-density ratio case, the surface tension parameter can be varied with negligible effects on both gas density and interface width.ω b can be varied independently of ω with negligible effect on this result.Significantly, this will allow its use for enhanced stability without impact on model parameters.

B. Surface tension
We have shown that the surface tension parameter, σ c , can be varied with only slight effects on the gas density, and that with the newly introduced correction accounting for the cascaded collision operator, this result is unaffected by varying the associated additional relaxation rates.We proceed by looking at varying surface tension at fixed interface width and then by varying interface width at fixed surface tension, again focusing on the effects of varying the additional relaxation rates.In the following results, in addition to measuring interface width, using Eq. ( 53), and surface tension, the average spurious velocity magnitude and the isotropy were   56), and average spurious velocity, ū, with ω 3 and ω 4 .Results are for the linear equation of state, Eq. ( 47), with ρ sat g = 0.001 (density ratio of 1000), and θ l = 1/3 and θ g = 1/6.Surface tension and interface width are kept approximately constant at σ m = 0.014 and W = 10.5, respectively.Results are for varying ω 3 (with ω 4 = 1) ( ū, black crosses; I , gray triangles) and varying ω 4 (with ω 3 = 1) ( ū, black diamonds; I , gray circles).also measured.Isotropy is defined as where r is the radius of the droplet, and subscripts indicate the angle to the x axis of the line along which radius is measured (radii being measured at the density midpoint in the diffusive interface).Surface tension was measured using the Laplace law which in 2D is given by 16. Variation in isotropy, I , defined in Eq. (56), and average spurious velocity, ū, with ω 3 and ω 4 .Results are for the linear equation of state, Eq. ( 47), with ρ sat g = 0.001 (density ratio of 1000), and θ l = 1/3 and θ g = 1/6.Surface tension and interface width are kept approximately constant at σ m = 0.00073 and W = 11.3,respectively.Results are for varying ω 3 (with ω 4 = 1) ( ū, black crosses; I , gray triangles) and varying ω 4 (with ω 3 = 1) ( ū, black diamonds; I , gray circles).
where P is the pressure difference between the inside and outside of the droplet.For this, droplets of radius 40 were initiated in the center of a 250 × 250 domain.Once equilibrium was reached the liquid and gas densities and radius were measured.Pressures can then be calculated from the equation of state and the surface tension then derived from Eq. (57).
surface tension reduced, by Eq. ( 51) and using σ ∝ σ c /W .Figure 11 shows the reduction in measured surface tension compared to the theoretically predicted surface tension, as surface tension is reduced at fixed interface width.Three cases are shown with ω = ω b = 1, ω = 1, and ω b = 1.6, and ω = 1.8 and ω b = 1.2.Results are normalized by the surface tension at κ = 0 (and ω = ω b = 1), σ 0 .Across all three cases the maximum and minimum measured interface widths were 11.5 and 10.5, respectively.Measured surface tension is seen to be in very good agreement with the theoretical prediction across the two orders of magnitude reduction over the κ = 0 case.The result is almost identical for the three combinations of ω and ω b .Figure 12 shows the measured average spurious velocity and droplet isotropy for the same cases.In general, spurious velocities are seen to decrease with decreasing surface tension, while anisotropy increases.This result is again unaffected by the choice of ω and ω b , with  20.Variation in isotropy, I , defined in Eq. (56), and average spurious velocity, ū, with sound speed c s = √ θ using the linear equation of state, Eq. (47), with ρ sat g = 0.001 (density ratio of 1000), and θ m = −1/512.Results are for varying θ g with θ l fixed at 1/3 (σ , black crosses; W , gray circles), and varying θ l with θ g fixed at 1/768 (σ , black diamonds; W , gray pluses).

C. Relaxation rates
To further quantify the results in the previous section, Figs. 13 and 14 show variations in spurious velocities and isotropy for varying ω b , for two different surface tensions, σ m = 0.014 in Fig. 13 and σ m = 0.00073 in Fig. 14.In both figures results are given for two different values of ω.All results show that ω b can be varied across a significant range with only minimal effect on the spurious velocities and isotropy.For the full range of results shown, the largest variation in measured surface tension and interface width was 3.4% and 3.0%, respectively (both in the σ m = 0.00073 case, at the lowest stable ω b ), with most variations being significantly smaller.
All results for surface tension and interface width have so far been given at ω 3 = ω 4 = 1.It is predicted from the third-order analysis that varying these higher-order relaxation rates will have no effect on model parameters.In one dimension almost no effect on density ratio and interface width was observed for varying these relaxation rates.To further investigate this we chose two cases: σ m = 0.014 with ω = ω b = 1, and σ m = 0.00073 with ω = 1.6 and ω b = 1.2.In both cases, ω 3 was varied between 0.2 and 1.8 with ω 4 = 1, and then ω 4 was varied between 0.2 and 1.8 with ω 3 = 1.The largest variation in measured surface tension and interface width across both of these ranges was 2.0% and 0.9%, respectively (both in the σ m = 0.00073 case, at ω 3 = 0.2).As predicted, the higher-order relaxation rates can therefore be varied with almost no effect on surface tension and interface width.However, they do impact the spurious velocities and isotropy, as shown in Figs. 15 and 16, which show measured spurious velocities and isotropy with varying ω 3 and ω 4 in each of the two cases.The large variations in spurious velocity are worth noting, as tuning the value of the higher-order relaxation rates can be an additional method for reducing spurious velocity.While there will be a tradeoff with stability, this is useful as in some cases it will allow the use of thinner interface widths (and therefore lower resolution simulations with reduced computational requirements).With the third-order correction derived above ω b can also be varied with minimal effect on model parameters such as density ratio and surface tension, and can therefore also be used to give slight reductions in spurious velocities, or more significantly to enhance stability.

D. Interface width
We now consider varying interface width at fixed surface tension.As interface width can be controlled independently, this can be used to systematically reduce spurious velocities [5].Here we show that this is also unaffected by the additional relaxation rates in the cascaded collision operator.
Figures 17 and 18 show the variation in spurious velocities and isotropy with measured interface width for the linear equation of state, Eq. ( 47 for ω = ω b = 1,ω = 1.0 and ω b = 1.6, and ω = 1.8 and ω b = 1.2.Surface tension is kept approximately constant at σ = 0.0128 in Fig. 17 and σ = 0.00074 in Fig. 18 (the largest variations in surface tension across the whole range of interface widths and relaxation rates were 1.7% for σ = 0.0128 and 7.8% for σ = 0.00074).In all cases the average spurious velocity is found to reduce with increasing interface width.This relationship is unaffected by changing ω and ω b .Having shown that interface width can be varied without significantly affecting the density ratio and surface tension, this allows systematic control over spurious velocities.ω b can be tuned for stability without changing this result.Anisotropy is found to be small in all cases, although in general is higher for the lower surface tension case.Again varying the relaxation rates does not have a significant effect.Some of the interface widths explored here are quite large.In practice, for simulating for example a droplet of radius R, it is required that R W .While it could therefore be impractical to use wide interfaces without the use of mesh refinement around the interface, the systematic decrease in FIG.24.Experimental results from Pan et al. [43] for the head-on splattering of equal-sized droplets at Re = 6210 and We = 440, as simulated in Fig. 26.
spurious velocity with interface width shown here provides this option, which was previously unavailable with pseudopotential methods.In any diffusive interface method such as the LBM, the numerical interfaces are already much larger than physical interface widths.While the results presented here are a significant advance for the pseudopotential multiphase LBM, it should be noted that other methods such as that by Lee and Fischer [42] have successfully reduced spurious velocities for comparable interface widths and density ratios.

E. Sound speeds
Finally, we briefly consider the possibility of tuning the sound speeds using the parameters θ g and θ l in the linear equation of state, Eq. ( 47).This is included for completeness, and no attempt is made to keep other model parameters constant while the sound speeds are varied.Figure 19 shows the variation in measured surface tension and interface width for two cases: varying the speed of sound in the gas (while the speed of sound in the liquid is fixed by setting θ l = 1/3) and varying the speed of sound in the liquid (while the speed of sound in the gas is fixed by setting θ g = 1/768), in both cases θ m is fixed at −1/512.The corresponding changes in spurious velocity and isotropy are shown in Fig. 20.It is important to note that variations in spurious velocity and isotropy are small, and that interface width and surface tension tend to constant as θ is decreased.If large values of θ were desired in simulation then surface tension and interface width can be tuned independently as discussed above.The theoretical sound speeds were compared with sound speeds measured from the one-dimensional decay of a density perturbation.Errors in measured sound speed were found to increase approximately exponentially as sound speed was decreased, for example, ranging from 0.1% at θ g = 1/6 to nearly 10% at θ g = 1/384 in the gas phase.Importantly, little or no dependance on θ m , κ or any of the relaxation rates, was observed.
Having established that all model parameters, including density ratio, surface tension, viscosity, interface width, and the additional relaxation rates of the cascaded collision operator can be independently varied, we now apply the model to the simulation of binary droplet collisions.

V. BINARY DROPLET COLLISIONS
Here the increased parameter range afforded by the current method is demonstrated through simulation of binary droplet collisions in three dimensions.The study of binary droplet collisions has many important applications across different scientific disciplines, from understanding cloud formation in climate theory, to engineering applications, such as turbine blade cooling, ink-jet printing, spray coating, and spray combustion in diesel internal combustion engines.Consequently, there have been a number of attempts to simulate such droplet collisions.In all methods, attaining high density ratio, low viscosity, and low surface tension simultaneously is a considerable challenge.The head-on collision of two droplets can be described by four dimensionless parameters, the density ratio between liquid and gas phase, ρ r = ρ l /ρ g , the viscosity ratio, ν r = ν l /ν g , and the Reynolds (Re) and Weber (We) numbers, defined as where R is the droplet radius and U their relative velocity.
Using the LBM Inamuro et al. [44] achieved droplet collision at Re = 2000 and We = 100 using a free-energy-based model, but at a density ratio of only 50.Similarly, Monaco et al. [45] achieved a Weber number of 760, at a higher density ratio of 150, but a lower Reynolds number of 200, using a Shan-Chenbased scheme.Lycett-Brown et al. [33] used the cascaded LBM with the Shan-Chen scheme at a similar density ratio of 120, with a higher maximum Reynolds number of 1200 but a lower maximum Weber number of 100.
For each of the simulations shown here, density ratio is fixed to 1000, to match that of water droplets colliding in air.For this the Carnahan-Starling equation of state is used with T = 0.056 and α = 0.02.The viscosity ratio is also set equal to that of water to air, ν r = 0.06, using a variable ω given by where ω l and ω g are the relaxation rate in the liquid and gas phases respectively.For all simulations we set ω b = ω 3 = ω 4 = 1.0.The interface width, W , is kept at approximately 5 lattice units, and the maximum recorded spurious velocity at any point in the domain (measured around the droplets at equilibrium) was 0.003.Droplet radius was set to R = 60, and simulations were run in periodic domains of dimensions around 10R, using an MPI parallelised code on up to 7000 cores.For the code used we found that the present method requires approximately 60% more computational time than the cascaded model with the EDM [the additional time being partly associated with the additional gradient calculations for the surface tension in Eq. ( 37)], which itself requires approximately 65% more computational time than the BGK model with the EDM.However the code is unoptimized and these figures are provided as a rough guideline only.
Figure 21 shows droplets coalescing with Re = 1380 and We = 27, and Fig. 22 shows droplet separation and production of one satellite droplet with Re = 1720 and We = 58.Both cases have been simulated successfully in previous studies, here we note that the relative size of the satellite droplet is in good agreement with the experimental result of Pan et al. [43], as shown in Fig. 23.Pan et al. [43] also produced experimental results for higher Weber number collisions.They showed a transition from separation along the axis of collisions, to breakup of the expanding disk in the axis perpendicular to the axis of collision as Weber number is increased.This results in satellite droplets from the breakup of the rim surrounding the central droplet.Further increase in Weber number then results in complete breakup of the colliding droplets into many smaller droplets, termed splattering, as shown in Fig. 24. Figure 25 shows our simulation results with Re = 4690 and We = 280, in which breakup of the expanding disk is captured.Additional breakup of the lamella is observed during contraction due to under-resolution.Galilean invariance is also observed in the rim breakup.The fingering on the interface observed in experiment is not seen here and should be the subject of future investigation.Figure 26 shows the case of Re = 6210 and We = 440, this was achieved with κ = 3.7, ω l = 1.982, and ω g = 1.75.An initial random oscillation and rotation were given to the droplets to break symmetry.Simulation of the splattering regime is achieved.
The purpose here is not to give a detailed analysis of droplet breakup, but to demonstrate the enhancements in stability of the present model.The current method has allowed significant increase in all parameters simultaneously.It is particularly significant that despite increasing the density ratio by an order of magnitude over previous LBM simulations, the Weber number and Reynolds number have both also been increased.This provides the ability to simulate real world droplet collisions, a more detailed study of which should be the focus of future work.

VI. CONCLUSION
A third-order analysis based on the Chapman-Enskog expansion has identified the errors in the pressure tensor in the pseudopotential cascaded lattice Boltzmann method (LBM).This has extended the recent analysis of Lycett-Brown and Luo [5] from the single-relaxation-time LBGK collision operator to the multiple-relaxation-time (MRT) cascaded collision operator.For the LBGK case the authors have shown that a third-order analysis is necessary to correctly identify the errors in the pressure tensor when large gradients in the forcing term are present at phase boundaries.The forcing term derived therein enabled the full range of coexistence curves to be accurately reproduced, even at arbitrarily high density ratios.When combined with the addition of a new term to the pseudopotential method, which allows variation of surface tension over a wide range, the resulting scheme allowed independent variation of density ratio, surface tension, interface width, and viscosity.Importantly, with this forcing scheme, errors in the coexistence curve decrease with increasing interface width (unlike in the exact difference method [18] and the scheme of Guo et al. [17], which give increasing errors as interface width is increased).While the new method does not significantly affect the spurious velocities found around curved interfaces when compared with the previous forcing schemes, it does allow them to be controlled (spurious velocities decreasing with increasing interface width).This opens the door to useful mesh refinement around interfaces in the pseudopotential multiphase LBM.
The previous analysis was limited to the LBGK collision operator, which is well known to be unstable at low viscosities, and has bulk viscosity fixed to be equal to shear viscosity.Introducing additional relaxation rates to allow separate control of bulk and shear viscosity and independent relaxation of higher-order moments can significantly enhance stability at low viscosities.Here the cascaded collision operator, which relaxes higher-order moments defined in a reference frame co-moving with the fluids is used, as this has shown further improvement over the standard MRT, which relaxes higher-order moments in the reference frame of the underlying lattice.Consequently, a further correction to the forcing term, dependent on the relaxation rate associated with the bulk viscosity (but not those of the third and fourth-order moments), has been derived.With this further correction to the forcing term it has been shown that the additional relaxation rates of the cascaded collision operator can be varied without impact on the previous results.Results are given for various equations of state, including a linear equation of state which further allows independent variation of liquid and gas sound speeds.
In summary, the pseudopotential cascaded LBM derived here is capable of varying all model and physical parameters (density ratio, surface tension, interface width, shear viscosity, bulk viscosity, the higher-order relaxation rates associated with enhanced stability, and liquid and gas sound speeds using the linear equation of state) virtually independently and over a significantly wider range than previously achieved.While the full exploration of this large parameter space should be the subject of future investigation, it is clear that this development will allow the pseudopotential multiphase LBM to tackle a wide range of multiphase flow problems in real-world applications.To illustrate this, we have presented results for binary droplet collisions.Density and viscosity ratio are matched to that of water droplets in air at room temperature and pressure (density ratio of 1000, kinematic viscosity ratio of 0.064), and collisions are given for Reynolds and Weber number up to Re = 6210 and We = 440, respectively.This is a significant improvement over the previous achievements of the pseudopotential LBM, of collisions at a density ratio of 120, equal kinematic viscosities, Re = 1190 and We = 100 [33].The results also bring the multiphase LBM in line with the cutting edge in CFD; recently Kuan et al. [46] achieved the first high Reynolds and Weber droplet collision simulation using an immersed boundary method, composed of a Eularian solver with adaptive mesh refinement for the liquid and gas phases, combined with a Lagrangian solver for the moving interface.They achieved the high parameters using elaborate procedures in contrast with the algorithmic simplicity of the present LBM.

ACKNOWLEDGMENT
Support is acknowledged from the UK Consortium on Mesoscale Engineering Sciences (UKCOMES) under the EPSRC Grant No. EP/L00030X/1.

APPENDIX A: CASCADED COLLISIONS
The collisions, f * i , for the 2D, 9-velocity, cascaded collision operator are derived in Lycett-Brown and Luo [19], and are given by where The central moments used in the cascaded collision operator, indicated by tildes, are related to the raw moments used in the MRT method by For the 3D, 27-velocity lattice the collisions, f * i , of the cascaded collision operator are derived and given in Lycett-Brown et al. [33].

APPENDIX B: CHAPMAN-ENSKOG EXPANSION
Starting from Eq. ( 27), the equations for conservation of mass and momentum, Eqs. ( 29) and (31), are derived.Equating terms of order ε 0 , ε 1 , and ε 2 in Eq. (27) gives When taking the velocity moments of these equations the following relations are used: for n 1, and for n 0. Additionally, at the order considered here, the following terms, which appear throughout the derivation, are equal to zero and therefore neglected: The zeroth, first, second, and third velocity moments of the zeroth order equation, Eq. (B1), give where in 2D ˜ n s = The zeroth-order velocity moment of the first-order equation, Eq. (B2), gives which using the derivative with respect to t 0 of Eq. (B7) and the divergence of Eq. (B8) becomes The first-order velocity moment of the first-order equation gives which using the derivative with respect to t 0 of Eq. (B8) and the divergence of Eq. (B9) becomes The second-order velocity moment of the first order equation gives which using the divergence of Eq. (B10) becomes The last term on the left-hand side comes from the cascaded expression of the third-order moment, given in Eq. (A3).As each of the terms within that last term are products of velocity and a second-order moment, we assume (for the low velocities required by the LBM) they can be neglected here.It should be noted apart from this term the derivation given here would apply to the MRT collision operator (removing tildes from moments), and the result should also therefore apply to the MRT method.
The zeroth-order velocity moment of the second-order term, Eq. (B3), gives which using the second derivative with respect to t 0 and the first derivative with respect to t 1 of Eq. (B7), the first derivative with respect to t 0 of Eq. (B8) and the divergence of Eq. (B9) becomes The first-order velocity moment of the second-order term gives which using the second derivative with respect to t 0 and the first derivative with respect to t 1 of Eq. (B8) and the divergence of Eq. (B18) becomes Combining Eqs.(B7), (B15), and (B21), we derive the familiar equation for the conservation of mass where the fluid velocity is now defined to take into account the additional terms in the conservation of mass equation that result from the third-order expansion, Combining Eqs.(B8), (B17), and (B23) and using this definition of velocity we can derive an equation for the conservation of momentum, Using the expression, where Tr indicates the trace, and assuming derivatives in time of F are smaller than spatial derivatives, then Eqs.(B9) and (B19), and their respective traces, can be used to rewrite the momentum equation as Taking into account the relations where T 0 ∇ρ = F has been used, and deriving the following for the trace terms, which is the momentum conservation equation given in Eq. ( 31).

40 FIG. 1 .
FIG.1.Simulation results for gas densities in one dimension, for the equation of state given by Eq. (46) with ρ sat l = 1.The method of Guo with A = 1/10 (crosses) and A = 1/40 (pluses), and the EDM with A = 1/10 (triangles) and A = 1/40 (diamonds) are shown compared with the expected density, ρ sat g (black line).

FIG. 4 .
FIG.4.Reduction in error in gas density with interface width, W , for the present forcing scheme with the equation of state given by Eq. (47).Gas density is measured far from a flat interface with ω = 1 in all cases.Four different density ratios are shown: ρ r = 10 (diamonds), ρ r = 100 (circles), ρ r = 1000 (crosses), and ρ r = 10000 (triangles).The solid black line shows second-order convergence.

FIG. 5 .ω = 1 , 6 FIG. 6 .
FIG.5.Simulation results for gas densities in one dimension, for the Carnahan-Starling equation of state with α = 1/16.The results are for the method of Lycett-Brown and Luo[5], which does not include the additional correction derived here for the cascaded collision operator, in the LBGK limit with ω b = ω = 1 (crosses), and for ω = 1,ω b = 1.6 (triangles), ω = 1,ω b = 0.4 (circles), ω = 1.8,ω b = 1.2 (diamonds), and ω = 1.8,ω b = 0.6 (squares).They are shown compared with the Maxwell construction (black line).As in Lycett-Brown and Luo[5] the LBGK results are in very good agreement with the Maxwell construction (greatly improving on the previous methods, as shown therein), but as expected there are significant errors in the gas density for the cascaded LBM with ω b = ω.(Liquid densities are all in close agreement with the Maxwell construction and are therefore not shown.)

3 FIG. 7 .
FIG. 7. Measured interface widths, W , for varying ω b , for the Carnahan-Starling equation of state with α = 1/32 and T = 0.0455 giving a density ratio of approximately 1000.Results are shown for κ = 0 (circles), κ = 1 (crosses), κ = 2 (triangles), and κ = 3 (diamonds).The solid line shows the theoretical prediction for the interface based on Eq. (52) and a measured interface width of 9.46 for the case ω b = 1 and κ = 0 (filled circle).All measured interface widths are in very good agreement with the predicted values.
) for the interface width coefficient, surface tension can be varied at constant interface width by varying κ and either α for the Carnahan-Starling equation of state or θ m for the linear equation of state.Both σ c and W c also depend on ω and ω b , therefore for both equations of state four different cases are shown: ω = ω b = 1, ω = 1 and ω b = 0.6, ω = 1 and ω b = 1.6, and ω = 1.8 and ω b = 1.2.In all cases very little change in gas density is observed, with the error reducing

T
is the trace of the momentum flux tensor, T = xx + yy in 2D, and T = xx + yy + zz in 3D.N indicates the normal stress differences, in 2D N = xx − yy is used, and in 3D we use N xz = xx − zz and N yz = yy − zz .From Eqs. (B6) and (B9), it can be shown that∂ ∂t 0 ∇ • ˜ 1 s = ∂ ∂t 0 ∇ • ˜ 1 b = 0, (B13)therefore these terms are neglected in the following derivation.