General mechanisms for stabilizing weakly compressible multiphase models and their applications

The present paper analyzes three representative weakly compressible multiphase models. It is found that these models contain some identical numerical dissipation terms, pressure diffusion, and bulk viscosity terms. Numerical investigations show that these identical numerical dissipation terms interpret general mechanisms for stabilizing computations. The generality of mechanisms is reflected in (a) they are likely to present in many weakly compressible multiphase models and (b) they represent interpretable physical mechanisms. Based on those general mechanisms, many weakly compressible multiphase models can be incorporated into a general theoretical framework. And a general weakly compressible solver for multiphase flows (GWCS-MF) is proposed. It is derived from standard governing equations and can easily be applied to nonuniform meshes. Detailed numerical tests demonstrate that it can achieve good numerical stability and accuracy for challenging conditions (inviscid fluids, large density ratios, high Weber numbers, and high Reynolds numbers) and can simulate complex interface evolution well. These good performances exhibit the advantages of GWCS-MF and further validate those general mechanisms.


Introduction
Incompressible multiphase flow is a common phenomenon in both nature and engineering.Compared with experiments, numerical simulation is a popular method to investigate complex phenomena in this field owing to its low cost, high efficiency, and generality.Numerical algorithms for incompressible multiphase flows can be roughly categorized into two groups.The first group is conventional solvers that directly solve macroscopic governing equations.The second group is mesoscopic lattice Boltzmann models describing discrete distribution functions' evolution.
Conventional solvers for incompressible multiphase flows fall into two categories according to approaches of updating the pressure.One is the exactly incompressible model [1,2] that updates pressure implicitly through the continuity equation.Without compressibility, it has no acoustic oscillation and, thus, is more likely to achieve good numerical stability.On the other hand, it mostly requires the solution of a Poisson equation for pressure or a pressure-correction equation through complicated iteration steps.Another one is the weakly compressible model [3][4][5] that updates pressure through an explicit pressure equation.However, it often requires additional treatments to stabilize computations.By combining them with methods to track the interface, such as volume of fluid [6], level-set method [2], front tracking algorithm [3], and diffuse interface method [4], conventional solvers can simulate nearly incompressible multiphase flows.
Multiphase lattice Boltzmann models are mesoscopic models based on the Boltzmann equation.They have been successfully applied to investigate practical engineering problems [7,8].These models can be roughly divided into four categories, color-gradient models [9], pseudopotential multiphase models [10,11], free energy models [12], and phase-field-based models [13].
These original models are unstable for multiphase flows with large density ratios (like water-air flow with a density ratio of about 1000) [14].Therefore, improved models based on different perspectives have been proposed to enhance the numerical stability at large-density ratios.
The first approach is the modification of these original models.For color-gradient models, the isotropic color gradient model [15] and the multiple-relaxation-time colorgradient lattice Boltzmann model [16] were proposed.For pseudopotential models, an improved forcing scheme for the multiple-relaxation-time model [17] was proposed to achieve thermodynamic consistency.It can handle large density ratios consequently.For free energy models, large density ratios beyond 1000 can be simulated by improving the Galilean invariance of higher-order terms [18].For interface tracking models, a stable discretization scheme based on directional derivatives [19] was proposed.All mentioned measures can be interpreted as modifications for higher-order terms.
The second approach is the simplification of lattice Boltzmann models.Two models, the multiphase lattice Boltzmann flux solver (MLBFS), and the simplified multiphase lattice Boltzmann model (SMLBM) [20], were proposed to improve the numerical stability at large density ratios.MLBFS [21] is a finite volume solver that reconstructs a simplified lattice Boltzmann model at the cell face to calculate pressure and momentum fluxes.The phase interface is captured by directly solving the macroscopic Cahn-Hilliard equation.Consequently, it can be applied to nonuniform meshes and was found to have good numerical stability for multiphase flows with high density ratios.The improved model [22] with a simplified computation of interface fluxes was proposed to remove the complex calculation of the compensation tensor and maintain numerical stability at a relatively thinner interface.Aiming to unify the different computational frameworks for velocity and phase fields, an interfacial lattice Boltzmann flux solver (ILBFS) [23] was proposed to solve the Cahn-Hilliard equation.Based on ILBFS, a simplified MLBFS with slight simplification in interface fluxes [24] and an improved MLBFS (IMLBFS) [25] based on the original phase-field-based multiphase lattice Boltzmann model (MLBM) [26] were proposed.
As SMLBM, it is a second-order approximation of the phase-field-based MLBM [26].The non-equilibrium distribution functions are approximated by interpolations of equilibrium distribution functions at different positions and moments.Therefore, its computational process involves only equilibrium distribution functions determined by macroscopic variables, which significantly decreases the memory size.SMLBM has also been proven to have good numerical stability at large density ratios, high Weber numbers, and high Reynolds numbers for 2D multiphase flows [20].
Although these MLBFS models and SMLBM have been successfully applied to simulate multiphase flows with large density ratios, the mechanisms for good numerical stability at large density ratios have not been well explained.The reason for the good numerical stability of the original MLBFS was thought to be using the lattice Boltzmann model [21].IMLBFS is believed to be more stable than the original MLBFS because it is derived from the multiphase lattice Boltzmann model and thus has a more robust physical basis [25].The reason for the good numerical stability of SMLBM is recognized as that SMLBM inherits the good stability feature of the reconstruction strategy [20].These hypotheses appear unclear and empirical, and clear evidence for the good numerical stability of these models has not been well established yet.
It should be noticed that the computational procedures of these MLBFS models and SMLBM involve only macroscopic variables, which indicates that they are intrinsically macroscopic models.Furthermore, the recovered governing equations of these models imply that they are weakly compressible models.From a macroscopic perspective, for both macroscopic weakly compressible models for single-phase [27] and multiphase flows [5], additional treatments [28,29] are needed, in general, to improve numerical stability.Therefore, it is reasonable to believe that mechanisms for the good numerical stability of MLBFS models and SMLBM can be established on the macroscopic scale.The present paper aims to reveal the general mechanisms that can incorporate these models into a general theoretical framework and then construct a general weakly compressible solver for multiphase flows (GWCS-MF).
In the present paper, macroscopic equations of a phase-field-based MLBM, IMLBFS, and SMLBM are derived first by approximating their actual computational process.Unlike continuous macroscopic equations recovered by using Chapman-Enskog expansion analysis, the key point of present studies is to recover the timediscretized macroscopic equations with numerical dissipation terms.Through detailed analyses, it is found that there are some identical numerical dissipation terms in these models, the pressure diffusion and bulk viscosity terms.The effect of these numerical dissipation terms on stabilizing computation is confirmed by numerical investigations, which indicate that these numerical dissipation terms indeed interpret general mechanisms for stabilizing computations.
GWCS-MF derived from standard governing equations is proposed based on the general mechanisms.Detailed numerical investigations confirm its good numerical stability and accuracy for multiphase flows with large density ratios, zero viscosities, high Reynolds numbers, and high Weber numbers.These results validate the general mechanisms further.
The remainder of this paper is organized as follows.Sections 2, 3, and 4 analyze the macroscopic equations of MLBM, IMLBFS, and SMLM, respectively.Section 5 summarizes these general mechanisms, proposes GWCS-MF, and makes numerical investigations to validate the general mechanisms.In Section 6, seven benchmark tests are simulated further to evaluate the numerical stability and accuracy of GWCS-MF and validate the general mechanisms.Finally, conclusions are given in Section 7.

Governing equations of the weakly compressible multiphase model
In the phase-field-based weakly compressible multiphase model, the following governing equations are adopted for simulations: where p is the pressure, ρ is the density, s c is the sound speed, u α is the velocity, µ is the dynamic viscosity, F α is the forcing term including the surface tension and body force.The popular phase field equation includes the Cahn-Hilliard equation, the Allen-Cahn equation, and their conservative forms [30][31][32].In the present paper, the Cahn-Hilliard equation is adopted, and it can be written as where C is the phase fraction of the heavier fluid, Once the interface width ξ and the surface tension σ are given, parameters λ and κ can be determined by The surface tension force can be evaluated by For a solid-fluid boundary [33], firstly, to ensure the mass conservation law, the boundary condition for C µ at a wall is wall 0 where n α is the unit outer normal vector.Secondly, to minimize the total free energy contributed by the specified wall free energy, the boundary condition for C is ( ) where ε is related to the equilibrium contact angle eq θ through For the interfacial area, the density is determined by ( ) Without a particular illustration, the kinematic viscosity is determined ( ) The subscripts H and L denote the parameters of the heavier fluid and the lighter fluid, respectively.

The phase-field-based MLBM
The phase-field-based MLBM has two evolution equations that describe the velocity and phase fields, respectively.The following model [26] is used for analysis, and the evolution equations are where i f and i g are pressure and phase fraction distribution functions of discrete velocity i e , respectively; eq i f and eq i g are the corresponding equilibrium distribution functions; f τ and g τ are dimensionless relaxation times for i f and i g , respectively; i H is the source term; x is the location; t is time; t δ is the time interval.For 2D and 3D situations, the D2Q9 model and the D3Q19 model are adopted, respectively.The discrete velocities of the D2Q9 and D3Q19 models are The equilibrium distribution functions eq i f and eq i g are ( ) ( ) ( ) where i w is the weight coefficient, p is the pressure, s c is the sound speed, u α is the velocity, and A is an adjustable parameter.The i w and s c of the D2Q9 and D3Q19 models are The equilibrium distribution functions satisfy The source term i H is where where Ma is the Mach number defined as s uc α .
The relaxation parameters f τ and g τ are respectively related to the kinematic viscosity υ and mobility Macroscopic variables p , u α , and C are respectively recovered by 2 0.5 0.5 Using the Chapman-Enskog expansion analysis, the phase-field-based MLBM can recover the governing equations [26], Eqs.(1) to (3), with second-order accuracy.

MEs-MLBM
Since MLBM is a discretized algorithm, the time-discretized macroscopic equations rather than continuous governing equations need to be investigated to explain mechanisms for good numerical stability.By using a second-order Taylor series expansion, ( ) , and i H can be expanded as , 0.5 , 0.5 , 0.5 Substituting Eqs. ( 28), ( 29), (30) to Eq. ( 13), we can obtain According to Eq. ( 13), we have ( ) ( ) which can also be rewritten as   (35) In Eq. ( 35), the last three lines are small deviation terms, and the term Finally, the viscous terms in the recovered momentum equations are 2 22 0.5 33 It can be seen that it has a bulk viscosity term 2 2 0.5 3 By using a second-order Taylor series expansion, ( ) , , 0.5 , 0.5 Substituting Eqs. ( 38), (39) to Eq. ( 14), we can obtain According to Eq. ( 14), we have ( ) ( ) which can also be rewritten as With the aid of Eqs. ( 20), (27), and (42), and taking summation of the zeroth-order moment of Eqs.(40), we can recover the macroscopic equation Eqs. ( 34), (35), and ( 43) are labeled by MEs-MLBM.
On the other hand, the time-discretized governing equations with the explicit firstorder scheme are It can be seen that compared with the discretized governing equations, MEs-MLBM contain a pressure diffusion term in the pressure equation and an additional bulk viscosity term in the momentum equation.

Analysis of the IMLBFS
In this Section, IMLBFS is introduced first, and then MEs-IMLBFS is derived.

IMLBFS
Using a Chapman-Enskog expansion analysis, the MLBM mentioned above can be rewritten as a finite volume scheme.By using some approximations to calculate nonequilibrium distribution functions, IMLBFS [25] was constructed.
The space-discretized equations of IMLBFS are ( ) where R α , αβ  (b) Calculation of the distribution functions at ( ) at the left cell at the right cell 0.5 at the interface where φ denotes an arbitrary variable, subscripts S , L , and R denote the variables at the cell face, left cell, and right cell, respectively.The partial derivatives at cell centers are calculated by the discrete Gauss theorem The macroscopic variables at the cell face in Eq. ( 51) are obtained by a linearized interpolation Once the macroscopic variables at  16) and (17), respectively.
(c) Calculation of the predicted variables at ( ) . The predicted variables at ( ) , , , Substituting Eq. (60) into Eq.(58) leads to the momentum flux Finally, substituting Eqs.(61) to (66) into Eqs.(47) to (49) leads to ( ) ( ) 1.5 0.5 where the viscous terms can be given as 2 22 33 It can be seen that MEs-IMLBFS also introduce a pressure diffusion term to the pressure equation and an additional bulk viscosity term to the momentum equation.

Analysis of SMLBM
SMLBM is a second-order approximation of MLBM.It has been proven to have good numerical stability for multiphase flow with large density ratios and high Reynolds numbers.

MEs-SMLBM
By using a second-order Taylor series expansion, ( ) ( ) ( ) As to macroscopic equations of the corrector step, using a second-order Taylor expansion leads to f Eqs. ( 89) to (94) are labeled as MEs-SMLBM.They contain a pressure diffusion term in the pressure equation.By combining Eqs. ( 90) and ( 93), all viscous-related terms can be summarized as It can be seen that MEs-SMLBM have a bulk viscosity coefficient

General mechanisms for stabilizing computations
Analyses in Sections 2 to 4 show that the two identical dissipation terms, the pressure diffusion and bulk viscosity terms, exist in these weakly compressible models for multiphase flows.It has been proven that the pressure diffusion term is efficient in stabilizing computations of both single-phase [34] and multiphase flows [5].The viscous terms also contribute to damping the pressure oscillations.Numerical experiments [35] confirm that sound waves experience the proper dissipation due to the intended bulk and shear viscosities.Therefore, it is thought that the two general dissipation terms provide general mechanisms for stabilizing computations of weakly compressible multiphase models.Generalities of mechanisms include: (a) they are commen for weakly compressible models for multiphase flows; (b) they represent physical mechanisms that can be applied for different discretization schemes, such as finite difference and finite volume schemes.

GWCS-MF
GWCS-MF adopts the computational procedures of IMLBFS, where numerical dissipation terms are introduced through a predictor-corrector step.The spacediscretized equations of GWCS-MF are ( ) The detailed procedures to calculate interface fluxes through reconstructing a local lattice Boltzmann model are introduced as follows:  method [36] to obtain the predicted variables at ( ) Then Eqs. ( 103) to ( 105) can be solved by different time discretization schemes.The explicit first-order scheme is adopted in the present paper to clarify the general mechanisms.Substituting Eqs. ( 106), ( 107), (109), and (110) to Eqs. ( 103) and ( 104), the macroscopic equations of GWCS-MF are ( ) It can be seen that it introduces a pressure diffusion term and an additional bulk viscosity term implicitly.

Dissipative and non-dissipative models based on the finite difference scheme
The dissipative model is constructed here to show that generality applies to different discretization schemes.The governing equations of the dissipative model are ( ) , respectively.
For simplicity, the governing equations are discretized by an explicit first-order scheme, and the discretized governing equations are ( ) ( ) where the subscripts n and 1 n + denote variables at time steps n and 1 n + , respectively.Partial derivatives in Eqs. ( 117) to (119) are discretized by the LSFD method [36].More details can be seen in Ref. [37].Note that the weight coefficients in the LSFD method are set to be proportional to 4  r − and 6 r − for 2D and 3D situations, respectively, where r is the distance between the neighboring nodes to the local node.

Numerical stability test
This section investigates the effect of these dissipation terms on stabilizing computation.The numerical stability of GWCS-MF, the dissipative model, the nondissipative model, MLBM, SMLBM, MLBFS, and IMLBFS is tested by simulating a steady 2D droplet immersed in the gas.Because the spurious velocity in the interfacial area is very small for this case at the steady state, it decreases the influence of different solution methods for the phase field equation on numerical stability.
In the computational domain , a liquid droplet of radius 0 R is placed in the center, while the surrounding is gas.Initially, the computational domain has a uniform pressure and zero velocity.The initial phase field is The convergence criterion is while the time steps of the other three models are 1.The Neumann boundary condition of zero gradients is adopted for all variables on all boundaries, and a uniform mesh of size 200×200 is adopted for all test cases.
Figure 3 exhibits the convergence situations of the four models.It can be seen that without the general dissipation terms, the non-dissipative model is divergent in all test cases.On the contrary, with these general dissipation terms, the dissipative model has a noticeable improvement in numerical stability.The comparison proves that these general dissipation terms are effective in stabilizing computation.Like the dissipative model, SMLBM is divergent at large kinematic viscosities while convergent at small kinematic viscosities.As for MLBM, it has good numerical stability at large kinematic viscosities while limited numerical stability at small kinematic viscosities.Second, the numerical stability of GWCS-MF, MLBFS, and IMLBFS is compared.The MLBFS analyzed here combines the improved model [22] and ILBFS [23].The predictor and corrector time steps are set as 0.5 t δ = , and 0.5 t ∆= for all models.Other parameters are assigned the same as above.As shown in Fig. 4, the three solvers can achieve good numerical stability at small kinematic viscosities with those general numerical dissipation terms.It means that the good numerical stability of MLBFS and IMLBFS can be well explained by those general numerical dissipation terms rather than the mesoscopic theory of the lattice Boltzmann method.Note that a smaller t ∆ can be adopted for the three solvers to obtain convergent results at large kinematic viscosities.In summary, the general numerical dissipation terms effectively stabilize computations and can well explain the good numerical stability of SMLBM, MLBFS, and IMLBFS.As for MLBM, it has other intrinsic mechanisms to achieve good numerical stability at large kinematic viscosities.However, it is unstable at small kinematic viscosities.It has been proven that the numerical stability of MLBM can be improved by using the multi-relaxation time model [38] that aims to adjust the higherorder unphysical moments and bulk viscosity.Thus, the instability of MLBM at small kinematic viscosities may result from the complex deviation terms shown in MEs-MLBM.

Numerical validations
Seven benchmarks are simulated in this section to validate GWCS-MF.

Laplace's law in inviscid fluids
To further investigate the performance of GWCS-MF at small kinematic viscosities, this section tests Laplace's law in inviscid fluids.The model has been depicted in Section 5.3.The liquid and gas are inviscid, and the density ratio is 1000.Other settings are the same as those in Section 5.3.The velocity and pressure contours are shown in Fig. 5.It can be seen that these contours are smooth even at such a large density ratio of 1000, which implies that GWCS-MF effectively suppresses the acoustic oscillation.The pressure differences at different surface tensions are compared to verify the correctness of GWCS-MF. Figure 6 exhibits that the present results are in good accordance with the analytical solutions determined by Laplace's law 0 pR σ

∆=
. The good numerical stability and accuracy for inviscid fluids imply the potential of GWCS-MF in simulating high Reynolds number multiphase flows, which can be seen in Section 6.3.On the contrary, the problem is a big challenge for multiphase lattice Boltzmann models, and no successful simulations have been reported to the best of our knowledge.
Fig. 6.Comparison of pressure differences at different surface tensions.

Two-phase Taylor-Couette flow in an annular area
The two-phase Taylor-Couette flow in an annular area is simulated to show the flexibility of GWCS-MF in handling curved boundaries.The schematic diagram is depicted in Fig. 7.The area 01 R rR ≤≤ is filled with lighter fluid, and the remaining space is filled with heavier fluid.The inner boundary has a fixed angular velocity ω , while the outer boundary is stationary.The steady velocity field is determined by the dynamic viscosity ratio The corresponding analytical solution is 1 An O-type body-fitted mesh is adopted for simulations, and its size is 160×80.The fixed computational parameters include 0 1 R = , 1 15 R = .

2D Rayleigh-Taylor instability
The 2D Rayleigh-Taylor instability problem is simulated to validate GWCS-MF for multiphase flows with complex interface evolution.In this problem, the heavier

GWCS-MF
upward, and two secondary roll-ups appear at the tails.The interface droplet configurations at different moments agree well with those in Refs.[13,21,23].
The positions and velocities of the spike tip and bubble front obtained by GWCS-MF are compared to those given by Wang et al. [21] using MLBFS and He et al. [13] using a phase-field-based lattice Boltzmann model to quantify and validate the present results.Note that the time and interface velocity here are nondimensionalized by dg and dg , respectively.As shown in Fig. 10, good agreement can be observed between the present results and reference results [13,21].
As to the situation of Re=2048, Fig. 11 exhibits the transient phase fields at different moments.At the early stage, phenomena like primary and secondary roll-ups are similar to those at Re=256.After 2.5 T = , the interface exhibits more complex evolution, and interface breakup occurs.The transient interface configurations also match those given by a phase-field-based lattice Boltzmann method [13].Comparisons of positions and velocities of the spike tip and bubble front are shown in Fig. 12.It can be seen that GWCS-MF results are in good accordance with those given by a phasefield-based lattice Boltzmann method [13].Taylor instability problem at Re=2048, together with those given by He et al. [13] using a phasefield-based lattice Boltzmann method.

2D bubble rising with a large density ratio
The example of a bubble rising is simulated to further validate GWCS-MF for tracking complex phase interface evolution with a large density ratio.The schematic diagram of the problem is shown in Fig. 13 , where g is the gravitational acceleration.The computational parameters follow the settings in Ref. [30] for comparisons.
The fixed parameters are 100 The buoyance is ( ) , where ρ is the local density.  .For a low Eo (Eo=10), which means a relatively larger surface tension, a smooth bubble droplet configuration like a semicircle is retained.As Eo increases, the constraint from the surface tension decreases.Consequently, at Eo=50 and 125, the bubble deformation becomes significant, and two symmetric tails form later.All bubble droplet configurations at different moments are in good accordance with those given by SMLBM [20].
As shown in Fig. 17, the present curves given by GWCS-MF agree well with the reference data [25,39].

Droplet splashing on a thin liquid film
In weakly compressible models, simulations of multiphase flows with large Weber numbers, large density ratios, and high Reynolds numbers are challenging owing to numerical instability issues [19].Therefore, the problem of a droplet splashing on a thin liquid film with a large Weber number of 8000, a large density ratio of 1000, and high Reynolds numbers up to 8000 is simulated to evaluate the performance of GWCS-MF in these challenging situations.The schematic diagram of this problem is depicted in Fig. 15.Initially, the liquid drop of diameter D is tangential to the liquid film surface.
It will impact the liquid film with a velocity ( ) Since the case is symmetric, only half the domain needs to be simulated.A uniform mesh of size 1000×500 is adopted for all simulations.Following the setting in Ref. the linear mean is adopted for the kinematic viscosity in the interfacial area.19 to 23, respectively.It can be seen that as the droplet hits the liquid film, the droplet deforms and tends to spread horizontally.The liquid film hinders the tendency.Thus, the drop periphery extends outward and is pushed upward by the surrounding liquid film simultaneously.The large viscous force prevents splashing when the Reynolds number is relatively low (Re=20).And an outward-moving surface wave can be observed in Fig. 19.For large Reynolds numbers (Re=100, 500, 2000, and 8000), the relatively small viscous forces cannot prevent the liquid from moving upward, and obvious splashing can be observed.These observations are consistent with those reported in Ref. [20].The transient dimensionless impact radius ( ) 2 rR is investigated to quantify the present results.The impact radius r is defined as the distance from the intersection of the unperturbed droplet and the original liquid film surface zH = [40].Theoretical analysis [40] indicates that for large Reynold numbers, the transient dimensionless impact radius approximately satisfies the power law ( ) ( ) , where A is a constant of about 1.0 [40].Many studies have also confirmed this regularity [19,20].As shown in Fig. 24, the present numerical results roughly satisfy the prediction of the power law, which verifies the correctness of GWCS-MF for multiphase flows with large Weber numbers, large density ratios, and high Reynolds numbers.

3D head-on collisions of binary microdroplets
Head-on collisions of binary microdroplets at different parameters are simulated to test GWCS-TF for 3D multiphase flows with large density ratios.The simulations mimic the experiments of tetradecane droplet collision in a nitrogen environment at 1 atm pressure conducted by Qian and Law [41] Only one-eighth of the computational domain is simulated owing to the model's symmetry.A uniform mesh of size 320×160×160 is adopted.Two cases are simulated, which correspond to We=32.8,Oh=0.615, and We=61.4,Oh=0.598.The liquid-gas density ratio is 666, and the dynamic viscosity ratio is 119 [42].Other parameters are set as 1 The transient phase fields of the two cases are depicted in Figs.(25) and (26).The droplet behavior is similar at the early stages of the two cases.The two droplets collide and merge into a larger one.And then, the larger droplet stretches in the x direction and forms a long liquid cylinder with two rounded ends.The later droplet behavior is different.For the case with We=32.8, which corresponds to a relatively large surface tension, the long droplet retracts inward and finally remains one droplet.In the case with We=61.4,the long liquid cylinder finally breaks into three droplets.The droplet evolution observed here matches the corresponding experimental results [41] well.The transient droplet configurations also agree with the numerical results [42] given by an axisymmetric multiphase lattice Boltzmann method.

Micro-droplet impacting a dry surface
Micro-drop impacting on surfaces exhibits various behaviors determined by many factors, including properties of the liquid, size and speed of the drop, and dropletsurface interaction characterized by static or dynamic contact angle(s).The complex outcomes and abundant experimental and numerical results of this problem make it a good benchmark test.In this section, GWCS-MF is applied to simulate it for further validation.
The present simulation mimics an experimental case made by Dong et al.HL µµ = , respectively.For the droplet impacting with such a high velocity, the gravity can be ignored in a short period.In the simulation, the computational domain is 22 Dx D − ≤≤ , 22 Dy D − ≤≤ , 02 zD ≤≤ .Since the model is symmetric, only a quarter of the domain needs to be simulated, and a uniform mesh of size 160×160×160 is adopted.Figure 27 exhibits the deformation process of the droplet.As the droplet hits the dry surface, it spreads to form a thin plate.At about t=10μs, the droplet deformation reaches the maximum, and then the droplet is pulled back by the surface tension and droplet-surface interaction.With time evolution, the droplet retracts to a cylinder and tends to jump.Furthermore, the dimensionless impact diameter D * and height H * given by GWCS-MF are compared with the experimental and numerical results.The impact diameter and height are respectively defined as the droplet diameter on the dry surface and the maximum height of the droplet, and they are nondimensionalized by D . Figure 28 shows that the present results match the experimental data [43] well and show better agreement than other numerical results [22,33].FIG.28.Comparison of impact diameter and height for the 3D droplet impacting a dry surface with Re=685 and We=103, θ eq =107°.

Conclusion
Many weakly compressible models based on lattice Boltzmann models have been proposed to achieve good numerical stability for multiphase flows with large density ratios.However, their mechanisms for stabilizing computation have not been well understood.The present paper aims to establish the general mechanisms that incorporate many weakly compressible multiphase models into a unified theoretical framework.
The present paper first derived the macroscopic equations of MLBM, IMLBFS, and SMLBM.Unlike the continuous equations recovered in previous references, the current derivation recovers the time-discretized macroscopic equations with numerical dissipation terms.It was found that the recovered macroscopic equations contain some common numerical dissipation terms.Numerical investigations prove that they provide general mechanisms for stabilizing computations.The generality of mechanisms applies to many weakly compressible multiphase models.They represent interpretable physical mechanisms that do not rely on specific discretization schemes.
GWCS-MF based on the finite volume scheme has been proposed based on the general mechanisms.The stable droplet immersed in gas was simulated by the dissipative model, the non-dissipative model, GWCS-MF, MLBM, SMLBM, MLBFS, and IMLBFS to compare their numerical stability.The results imply that the general mechanisms can well explain the good numerical stability of SMLBM, MLBFS, and IMLBFS.As to MLBM, it has other intrinsic mechanisms to achieve good numerical stability for large kinematic viscosities.However, it is unstable for small kinematic viscosities, which may result from complex deviation terms shown in MEs-MLBM.
Seven benchmark cases were simulated to evaluate GWCS-MF in detail.The results prove that it can achieve good numerical stability and accuracy for multiphase flows at challenging conditions, such as inviscid fluids, large density ratios, high Weber numbers, and high Reynolds numbers.It is also flexible in nonuniform meshes and can tackle complex interface evolutions well.These good performances prove the advantages of GWCS-MF and further validate the general mechanisms.

CM
is the constant mobility, and C µ is the chemical potential determined by the total free energy of fluid-fluid or fluid- wall interfaces.The total free energy is

( ) 0 E−
C is the bulk free energy defined as , λ and κ are two fixed parameters; ( )S C ϕis the wall free energy per unit area.For the equilibrium state, the total free energy is minimized, and the chemical potential C µ of the inside fluid can be determined by

Π
and Q α are interface fluxes, s c is the sound speed, V ∆ is the control volume, subscript k denotes the th k interface of the control volume, S α ∆ and k n α are the area and the outer normal vector of the th k interface, respectively.The detailed procedures to calculate interface fluxes through reconstructing a local lattice Boltzmann model are as follows: (a) Reconstruction of the local lattice Boltzmann model.As shown in Fig. 1, a unit lattice of the D2Q9 model is constructed at the finite-volume cell face, and S x is the position of the face center.The discrete velocities of the D2Q9 model and the D3Q19 model are given in Eq. (15).

Fig. 1 .
Fig. 1.The reconstructed unit lattice at the cell face.
respectively calculated by Eqs. ( Eq. (61) into Eq.(59) leads to the volume fraction flux )Eqs.(62) to (65) and Eqs.(68) to (70) are labeled as MEs-IMLBFS.To show the numerical dissipation terms included in MEs-IMLBFS, Eqs.(62) and (63) are substituted into Eqs.(68) and (69), and we can obtain (a) Reconstruction of the local lattice Boltzmann model.As shown in Fig. 3, a unit lattice is constructed at the cell face and S x is the position of the face center.In the present paper, the lattice size is set as the minimum value among SL − xx , SR − xx , SU − xx , and SD − xx .

Fig. 2 .
Fig. 2. The reconstructed unit lattice at the interface.
µ are the adjustable pressure diffusion coefficient, and bulk viscosity coefficient, respectively.Referring to Eqs. (112) and (113), p D and b The non-dissipative model is similar to the dissipative model, except that the p D numerical stability of four models based on the finite difference scheme, the non-dissipative model, the dissipative model, SMLBM, and MLBM, are investigated.Subscripts H and L represent the parameters of the liquid and gas, respectively.The fixed parameters include 0an extensive range from 1 to 1000.The predictor and corrector steps of the dissipative model are set as 0

Fig. 3
Fig. 3 Stability of the non-dissipative model (a), the dissipative model (b), SMLBM (c), and MLBM.The blue circle indicates stable, and the red cross indicates unstable solutions.
considered.The nonslip boundary condition is imposed on the inside and outside boundaries.

Figure 8 Fig. 8 .
Fig. 8.Comparison of the azimuthal velocity profiles at different dynamic viscosity ratios for the two-phase Taylor-Couette flows in an annular area.

1 L
fluid of density H ρ floats on a lighter fluid of density L ρ .Once a slight pressure or interface curvature perturbation occurs, the interface becomes unstable, and the two fluids would reverse. is the characteristic length.The nonslip boundary condition is imposed on the upper and lower boundaries, and the periodic boundary condition is set on the left and right boundaries.The Reynolds number Re and the Atwood number At are used to characterizing the problem.They are respectively defined as of size 200×800 is adopted.The Atwood number is set as 0.5, and Re has two values, 256 and 2048.Other simulation parameters are set as Re=256 are discussed first.Figure 9 shows the transient phase fields at different moments.At the early stage, the central heavy fluid falls to form a spike, and the lighter fluid at the two sides rises to form bubbles.At 2 T = , two symmetric roll-ups occur at the spike end.As the evolution continues, the two roll-ups stretch

Fig. 9 .
Fig. 9. Evolution of the fluid interface for the 2D Rayleigh-Taylor instability problem at Re=256.

Fig. 10 .
Fig. 10.Transient positions and velocities of the spike tip and bubble front for the 2D Rayleigh-Taylor instability problem at Re=256, together with those given by Wang et al. [21] using MLBFS and He et al. [13] using a phase-field-based lattice Boltzmann method.
. The computational domain is DxD −≤≤ , 3 Dy D −≤≤ .A stationary bubble of diameter D is initially immersed in the heavier fluid, and the bubble center is ( ) 0, 0 .The upper and lower boundaries are set as the nonslip boundary condition with zero velocity, while the left and right boundaries are set as the periodic boundary condition.The problem is characterized by three dimensionless parameters: the density ratio HL ρρ ,

Figures 14 to 16
Figures 14 to 16 depict the evolution of the bubble interface defined as the contour 0.5 C =.For a low Eo (Eo=10), which means a relatively larger surface tension, a smooth bubble droplet configuration like a semicircle is retained.As Eo increases, the constraint from the surface tension decreases.Consequently, at Eo=50 and 125, the bubble deformation becomes significant, and two symmetric tails form later.All bubble droplet configurations at different moments are in good accordance with those given by SMLBM[20].
center at Eo=125 are compared with the reference results[25,39] to validate the present results in quantity.The vertical position and velocity of the bubble mass center are defined as

Fig. 17 .
Fig. 17.Transient vertical position and velocity of the bubble mass center at Eo=125, togetherwith the reference data given by Aland and Voigt[39], and Yang et al.[25].

[ 19 ]
, the computational parameters are set as 0, 100, 500, 2000, and 8000, are chosen.The dynamic viscosity of the lighter fluid is fixed, and the viscosity ratio HL µµ is 40 at Re = 500.Referring to Ref.[19],

Fig. 18 .
Fig. 18.Schematic diagram of the droplet splashing on a thin film.

Fig. 19 .
Fig. 19.Evolution of the phase field for droplet splashing on a thin liquid film at Re=20.

Fig. 20 .
Fig. 20.Evolution of the phase field for droplet splashing on a thin liquid film at Re=100.

Fig. 21 .
Fig. 21.Evolution of the phase field for droplet splashing on a thin liquid film at Re=500.

Fig. 22 .
Fig. 22. Evolution of the phase field for droplet splashing on a thin liquid film at Re=2000.

Fig. 23 .
Fig. 23.Evolution of the phase field for droplet splashing on a thin liquid film at Re=8000.

Fig. 24 .
Fig. 24.The transient impact radiuses at different Reynolds numbers for droplet splashing on a thin liquid film.
boundaries, the Neumann boundary condition of zero gradients is adopted for all variables.Initially, the two droplets of Diameter D are placed at ( ) The problem is characterized by the Weber number We and the Ohnesorge number Oh defined as
[43] for comparison.Initially, the droplet of diameter 0 50.5μmD = is tangential to the dry surface with an equilibrium contact angle 107 eq θ =  and will impact the dry surface with a vertical speed of 12.2m/s.The case corresponds to Re=685 and We=103, where Re and We are defined in Eq. (126).The density and viscosity ratio of the droplet and

FIG. 27 .
FIG. 27.Evolution of interface morphology for the 3D micro-droplet impacting on a dry surface with Re=685 and We=103,