A three-dimensional non-orthogonal multiple-relaxation-time phase-field lattice Boltzmann model for multiphase flows at large density ratios and high Reynolds numbers

This study proposes a three-dimensional non-orthogonal multiple-relaxation-time (NMRT) phase-field multi-phase lattice Boltzmann (PFLB) model within a recently established unified lattice Boltzmann model (ULBM) framework [Luo et al., Phil. Trans. R. Soc. A 379, 20200397, 2021]. The conservative Allen-Cahn equation and the incompressible Navier-Stokes (NS) equations are solved. In addition, a local gradient calculation scheme for the order parameter of the Allen-Cahn equation is constructed with the non-equilibrium part of the distribution function. A series of benchmark cases are conducted to validate the proposed model, including the two-phase Poiseuille flow, Rayleigh-Taylor instability, binary liquid/metal droplet collision, and a bubble rise in water. The present simulation results are in good agreement with existing simulation and experimental data. In the simulation of the co-current two-phase Poiseuille flow, the present model is proven to resolve the discontinuity at the phase interface and provide accurate results at extremely high density ratios (i.e., up to 10 6 ). Finally, the proposed model is adopted to simulate two challenging cases: (1) water droplet splashing during its impacting on a thin liquid film and (2) liquid jet breakup. The simulation results demonstrate an excellent agreement with previous experimental results, both qualitatively and quantitatively. In these simulations, the Weber number and Reynolds number reach 10 5 and 6000, respectively, and the viscosity can be as low as ∼ 10 (cid:0) 4 , in the lattice unit.


Introduction
Multiphase fluid flow is ubiquitous in science and engineering and has received considerable research attention in the past few decades.A wide range of applications can benefit from a thorough understanding of multiphase phenomena, such as the anti-icing of the turbine blade (Cao et al., 2009;Sarkar and Farzaneh, 2009), inkjet printing (Castrejón-Pita et al., 2008;Gan et al., 2009) and spray cooling (Kim, 2007;Liang and Mudawar, 2017;Liang and Mudawar, 2017).Multiphase interfacial dynamics is a typical nonlinear, nonequilibrium and multiscale phenomenon.The complexity of multiphase flow presents a significant challenge to both analytical and experimental approaches due to the high dimensionality of the physical parameter space, hindering attempts to uncover the underlying mechanisms.Advanced numerical methods are therefore desired to supplement the analysis and experiment (Brennen and Brennen, 2005).Owing to the mesoscopic nature, high computational efficiency and simple algorithm, the lattice Boltzmann (LB) method has been regarded as a promising numerical approach to model multiphase fluid flow (Benzi et al., 1992;Qian et al., 1995).One significant feature of LB method is that it is fully explicit.In the simulation of incompressible multiphase flow.It does not require solving the pressure Poisson equation, which is usually computationally expensive.Numerous multiphase LB models have been proposed for multiphase flow, such as the colour gradients model, the free energy model, the pseudopotential model, and the phase field model (Huang et al., 2015;Li et al., 2016).Amongst them, the phase field model has attracted much interest because of its potential to capture large interface topological deformations and to model complex instability dynamics (Wang et al., 2019).
In the phase field LB model, two evolution equations of distribution functions are incorporated to solve the interface tracking equation and the incompressible Navier-Stokes (N-S) equations.To obtain accurate solutions of these two equations, improved LB models have been continuously proposed.On the one hand, for the interface tracking equation, there are two widely-used governing equations, namely the Cahn-Hilliard (CH) equation (Cahn, 1959) and the Allen-Cahn (AC) equation (Allen and Cahn, 1976).The CH equation has been commonly used in the early-stage phase-field LB models (He et al., 1999;Lee and Lin, 2005) because it satisfies the mass conservation by nature.Nevertheless, the CH equation involves a fourth-order derivative of the phase indicator, the solution of which requires extra efforts using the LB method (Liang et al., 2023).In contrast to the CH model, the AC equation does not include high order derivatives, which makes it more compatible with the LB model (Wang et al., 2019;Mitchell et al., 2018).Initially, the use of AC interface equation is limited because it does not satisfy the mass conservation.This problem is solved later by reconstructing the equation to be conservative (Chiu and Lin, 2011;Sun and Beckermann, 2007).As pointed out in Refs.(Liang et al., 2023;Wang et al., 2016), the AC equation is superior in computational efficiency, numerical stability, and simplicity over the traditional CH interface tracking equation.Therefore, the AC equation becomes popular in the phase field LB model.To solve the conservative AC equation, Geier et al. (Geier et al., 2015) proposed a central moment-based LB model.However, both Ren et al. (Ren et al., 2016) and Wang et al. (Wang et al., 2016) claimed that Geier's model cannot fully recover the target AC equation, where the time differentiation should be included in the source term.Accordingly, an improved multiple-relaxation-time (MRT) model is developed in Ref. (Ren et al., 2016).In some recent studies, the efficiency and the accuracy of the LB model for the AC equation are further improved using, for example, a nondiagonal relaxation matrix (Zhang et al., 2022;Xu et al., 2021), the local gradient calculation scheme (Zhang et al., 2022;Liang et al., 2018), or improved collision operators (De Rosis and Enan, 2021).However, challenges remain in the LB model of the AC equation, especially in achieving high-order accuracy and simulating the narrow interface.
For the incompressible NS equations, two LB models are proposed, i. e., the momentum-based model and the velocity-based one.For the momentum-based model, a modified equilibrium function is proposed, which integrates the hydrodynamic pressure as the primary quantity of interest (Lee and Lin, 2005;He et al., 1999).By contrast, in the velocity-based LB model proposed by Zu and He (Zu and He, 2013), the equilibrium distribution functions recover fluid velocity-based macro dynamic equations rather than the fluid momentum-based ones.Recently, a central moment-based cascaded LB model is developed for the velocity-based model (De Rosis and Enan, 2021;Gruszczyński et al., 2020).It has been pointed out that the momentum-based model is more efficient in parallel computing than the velocity-based one because it does not require a prediction-correction scheme to update the pressure and it has no velocity derivatives in the interface forces (Fakhari et al., 2017).Besides, the momentum-based model has better potential in simulating large viscosity ratio and density ratio as it follows the H-theorem (Zu and He, 2013).Moreover, the momentum-based model has been improved by using a weighted MRT collision operator (Fakhari et al., 2017), coupling large eddy simulations (An et al., 2021), and proposing a simplified force term (Liang et al., 2018).Nevertheless, discontinuity at the phase interface is produced in the momentum-based model, which causes velocity fluctuations and hampers its applications in simulating large-density ratio multiphase flow (Zu and He, 2013).
Despite continued progress in the development of phase field LB models, there is still room for further improvements in order to simulate challenging multiphase flow problems, especially realistic multiphase flows at large density ratios and high Reynolds numbers.For high density ratio (~10 3 ) multiphase flow, the achievable Reynolds number in the existing phase field LB model is typically around (~10 2 ).Additionally, comprehensive comparisons between numerical and experimental results are still lacking.In this work, we proposed a threedimensional phase-field LB model within a recently proposed unified lattice Boltzmann model (ULBM) framework (Luo et al., 2021).In our model, the conservative AC equation (Chiu and Lin, 2011;Sun and Beckermann, 2007) is used for the interface tracking, coupling with the incompressible NS equations for the fluid flow.A non-orthogonal moment set (Fei et al., 2018;Fei et al., 2019) is used to construct the non-orthogonal MRT LB model (NMRT) to solve the target governing equations.The proposed model is then used to reproduce various complex multiphase experiments, including binary liquid/metal droplet collision, bubble raising, water droplet splashing and liquid jet break.
The outline of this paper is as follows: we begin with an introduction of the conservative phase field multiphase model and the ULBM framework in Section 2.1.Then, the constructed ULBM (NMRT) model for the interface tracking equation and incompressible NS equations are given in Sections 2.2 and 2.3, respectively.In Section 3, we provide comprehensive validations and assessments for the proposed ULBM (NMRT-PFLB) PF model.The conclusions of this work are given in Section 4.

The ULBM (NMRT) and the conservative phase field multiphase model
In the present phase field multiphase LB model (i.e., ULBM (NMRT) PF ), a conservative AC equation is used for interface tracking (Chiu and Lin, 2011): where ϕ is the order parameter (phase indicator).In the two-phase system, the heavy phase and the light phase are indicated by the order parameter ϕ h and ϕ l , respectively.In this work, the value of ϕ h for pure heavy phase is set as 1 and that of ϕ l for pure light phase is set as 0, with ϕ 0 = (ϕ h +ϕ l )/2 = 0.5 indicating the position of the interface.n is the unit vector normal to the interface, which is defined as n = ∇ϕ/|∇ϕ|.M ϕ is the mobility and w stands for the interface thickness.u = [u x , u y , u z ] is the fluid velocity.In the phase field multiphase model, the conservative AC equation is coupled with incompressible NS equations, and the continuity and momentum equations for incompressible multiphase flows can be expressed as (Unverdi and Tryggvason, 1992): ρ indicates the fluid density, ν is the fluid kinematic viscosity, ν b stands for the non-hydrodynamic bulk viscosity and P is the pressure; F s and F b are the surface tension and body force, respectively.The surface tension F s is determined by: (3) (2) G. Wang et al. tial, and the parameters β = 12σ/w and k = 3σw/2 are related to the interface thickness w and the surface tension σ.In this work, we adopt the ULBM (NMRT) framework for the phase field multiphase model.The general evolution equation of the collision operator with forcing term in the ULBM framework can be written as (Luo et al., 2021): f i and f * i are the pre-collision and post-collision distribution functions, respectively.f eq i is the equilibrium distribution function and C i is the forcing term.e i and Δt = 1 are the discrete velocities and the time step, respectively.I is the unit matrix, and S represents the relaxation matrix.It should be noticed that e i and S are varied accordingly based on the adopted lattice model (i.e., D3Q19 model and D3Q27 model in this work).The details can be found in Appendixes 2 and 3, respectively.
In the ULBM framework, the transformation matrix M is adopted to transform the distribution functions (f i ) to their raw moments (m), and the shift matrix N is used to convert the raw moments into the central moments ( m), the transformation/shift can be expressed as: (5) As pointed out in the recent work of Mitchell et al., (2021), for the phase-field LB model, the central moment-based collision operator increases the numerical dispersion.Therefore, the raw moment-based collision operator is adopted, which corresponds to the shift matrix N = I.Accordingly, the collision step in the raw moment space can be re-written as: m eq = Mf eq i and R = MC i are the discrete equilibrium moment and discrete forcing term in the raw moment space, respectively.In this study, we adopt a simplified non-orthogonal moment set which is originally proposed in Refs.(Fei et al., 2018;Fei et al., 2019) to construct the ULBM (NMRT) model.We first define the raw moments as:

∑
f i e p ix e q iy e n iz , p, q, n ∈ {0, 1, 2}, where u x , u y , and u z are the velocity components in x, y and z directions, respectively.The non-orthogonal moment set (m) in raw moment corresponding to D3Q19 and D3Q27 lattice models can be found in Appendixes 2 and 3, respectively.Then, by substituting Eq. ( 7) into m, the transformation matrix M can be obtained, and its explicit expression can be found in the Appendix.
It is noted that compared with the traditional orthogonal moment set, there are more zero terms in M and M − 1 , and the sub-lattice model (e.g., D3Q19) can be deduced from the full-lattice model (D3Q27) directly (as shown in Appendix 3).Furthermore, in traditional MRT-LBM phase field models with conventional orthogonal moment sets, the discrete equilibrium moments and forcing terms are often very complex.In particular, in the three-dimensional case, it may contain hundreds of terms.In contrast, our ULBM (NMRT) model leads to a concise expression.Therefore, our proposed model exhibits favourable features such as a simplified implementation and improved computational efficiency.The above ULBM (NMRT) model has been widely adopted in the pseudopotential multiphase model.It has been demonstrated that the ULBM (NMRT) model has advantages in terms of computational efficiency, flexibility and numerical robustness (Wang et al., 2020;Wang et al., 2021;Wang et al., 2020).

ULBM (NMRT) model for the conservative Allen-Cahn equation
In this section, we integrate the conservative AC equation into the ULBM (NMRT) framework.The collision operator in raw moment space is written as: The subscript ϕ represents the phase indicator.To recover the target AC Eq. (1), differing from Refs.(Fakhari et al., 2017;Fakhari et al., 2017), we choose the equilibrium distribution function as: where ω(|e i | 2 ) is the weight, and C s = 1/ ̅̅̅ 3 √ stands for the lattice sound speed.The discrete equilibrium moment in the raw moment space (m eq ϕ ) can be achieved by multiplying the transformation matrix M with f eq ϕ,i .The forcing term is written as (Ren et al., 2016): where F ϕ,i is given by: Similarly, by multiplying the transformation matrix M, the discrete forcing term in the raw moment space is R ϕ = MC ϕ,i .The details of the D3Q19 and D3Q27 ULBM (NMRT) model for the AC equation can be found in Appendix.The order parameter can be calculated by f ϕ,i is reconstructed by f ϕ,i = M − 1 m * ϕ after the collision in the raw moment space.The relaxation parameters are s ϕ,1 = s ϕ,2 = s ϕ,3 = 1/(M ϕ /c 2 s + 0.5) and the other relaxation parameters can be chosen freely.To save the computational cost, we set the other relaxation parameters as 1.Applying the Chapman-Enskog (CE) analysis in Appendix 1, it can be proven that the above LB model can recover the target AC equation Eq. (1) without high-order error terms.It is worth mentioning that the gradient terms in Eq. ( 11) can be calculated by the second-order lattice-based finite difference (FD) scheme: And the Laplace operator in Eq. ( 3) is given by: The time derivative in Eq. ( 25) is then calculated by the Eulerian scheme, i.e.,

ULBM (NMRT) model for the incompressible Navier-Stokes equations
In this study, the momentum-based LB model is used for incompressible flow.Motivated by Refs.(Lee and Lin, 2005;He et al., 1999), when solving the incompressible NS equations in the phase-field multiphase model, a modified distribution function is now introduced: where P 0 = ρc 2 s stands for the ideal gas equation of state.The equilib-G.Wang et al. rium distribution function is where Γ i is defined as: substituting Eq. ( 16) into the discrete Boltzmann equation for the particle distribution function: D/Dt = ∂ t + e i .∇ is the total derivative and Ω g stands for the collision term for g i .F is the total force, which can be expressed as (He et al., 1999;Fakhari et al., 2017;Fakhari et al., 2016): P is constant in the nearly incompressible limit, and P 0 = ρc 2 s is also constant as the conservative AC equation is adopted.Thus (Fakhari et al., 2017;Fakhari et al., 2016): F g = F s + F b , and therefore the corresponding collision operator of the ULBM (NMRT) model for incompressible NS equation is given as: where the discrete equilibrium moments in the raw moment space can be calculated as m eq g = Mg eq i .Similarly, the discrete forcing terms in the raw moment space can be given by Eq. ( 21) and M, where R g = MC g,i .Notably, we simplify the discrete forcing terms by ignoring the high order terms (~ O(Ma 3 ), related to u 3 and uuF g ).The explicit expression of m eq g and R g for the D3Q19 lattice model and D3Q27 lattice model can be found in Appendixes 2 and 3, respectively.
In this study, the fluid density is calculated by interpolating the order parameter as, where ρ h and ρ l stand for the density of the heavy phase and the light phase, respectively.According to Eq. ( 23), the derivatives of P 0 can be determined by ϕ, where: In most previous phase-field LBM models, the gradient in Eq. ( 24) is given by the lattice-based FD scheme (i.e., Eq. ( 13)) (Wang et al., 2019).Alternatively, some recent works (Wang et al., 2016;Zhang et al., 2022;Liang et al., 2018) pointed out that the gradient term can be deduced locally from the non-equilibrium part of the distribution functions.Through the CE analysis in Appendix 1, the gradient of ϕ in the current ULBM (NMRT) PF model relates to the non-equilibrium part of the m ϕ : ) ) .
(25) This is the local scheme for the gradients of ϕ in Eq. ( 24), which is designated as the Local scheme in this paper.Although the above equation needs implicit terms to calculate F ϕ , it can be made explicit by algebraic derivation (e.g., Eq. ( 20) in Ref. (Zhang et al., 2022)).Geier et al. (Geier et al., 2015) pointed out that undesired grid-scale oscillations are generated when using the Local scheme for the interface tracking equation.Consequently, in our work, the Local scheme Eq. ( 25) is only adopted when solving the incompressible NS equations (i.e., in Eq. ( 24)).The gradient terms involved in the interfacing tracking equation (i.e., in Eq. ( 11)) are still calculated by the FD scheme.In the following Section 3.1.1,we will show that this setup can solve the interface discontinuity caused by using the FD scheme when solving the NS equation.
Following the collision step in raw moment space (Eq.( 22)), the distribution function can be reconstructed by g i = M − 1 m * g .The macroscopic variables are calculated as: ) . ( With the CE analysis in Appendix 1, the above ULBM (NMRT) PF model can recover the target NS equations in the nearly incompressible limit.The relaxation parameters s v = [s g,4 , s g,5 , s g,6 , s g,8 , s g,9 ] are related to the fluid kinematic viscosity, where ν = (1 /s v − 0.5)C 2 S Δt.Besides, the bulk viscosity ν b depends on s g,7 , and ν b = 2/3(1 /s g,7 − 0.5)C 2 S Δt.All the other relaxation parameters can be chosen freely.In the following simulations, without any statement otherwise, we set s g,0 = s g,1 = s g,2 = s g,3 = s g,7 = 1 and the rest of high-order (i > 9) relaxation parameters as 1.25.To avoid the implicit velocity, the gradient terms involved in Eq. ( 26) are also calculated by Eq. ( 13).

Model validations
In this section, we validate the accuracy and stability of the proposed ULBM(NMRT) PF model through several benchmark cases.In subsection 3.1.1,we conduct the Laplace test of a static droplet to assess the model stability.Then in subsection 3.1.2,we validate the accuracy of our model by simulating a co-current two-phase Poiseuille flow.Finally, in subsection 3.1.3,the Rayleigh-Taylor instability (RTI) is simulated to verify the ability of the proposed model for modelling multiphase flows with complex interfacial deformations.

Laplace verification and static droplet test at large density ratios
In the Laplace test, a static droplet with an initial radius R 0 is placed at the centre of a 4R 0 × 4R 0 × 4R 0 box, with periodic boundaries in all directions.The simulation is first conducted using the D3Q19 ULBM (NMRT) PF model with the Local scheme.The initial density field is described by the following function: where X stands for the distance from the droplet centre.In this case, we choose w = 5, M ϕ = 0.01.The density ratio and dynamic viscosity ratio are set as 1000 (ρ h = 1000, ρ l = 1) and 1 (ν h = ν l = 0.1), respectively.
In this case, R 0 varies from 12 to 50, for three different surface tensions ρ ana.|/ρ ana.× 100%) between the simulation results and the analytical solution for the density distribution.It can be found that for narrower interface cases, the maximum relative errors are larger than those for wider interface cases, while the differences between the FD scheme and the Local scheme are small.The average relative errors for the cases of w = 3 and w = 5 are 2.0% and 1.2%, respectively.We then conduct a static droplet test to compare the numerical performance of different LB models, including the D3Q19 model with FD scheme (D3Q19 FD), the D3Q19 model with Local scheme (D3Q19 Local), the D3Q27 model with Local scheme (D3Q27 Local), and the single relaxation time D3Q19 model with Local scheme (D3Q19 SRT).The difference between the Local scheme and the FD scheme is reflected in the calculation of Eq. ( 24).Similar to the setup in the Laplace verification, a static droplet with R 0 = 30 is placed in the centre of a 4R 0 × 4R 0 × 4R 0 box, and the density ratio is set as 1000.Table 1 records the maximum spurious velocities for different simulation parameters.The

Table 1
Comparison of maximum spurious velocities in the gas phase under different simulation parameters, by using D3Q27 and D3Q19 ULBM (NMRT) PF.The last row stands for the total CPU time of 20,000Δt for various cases.that the D3Q27 MRT model can remain stable when ν l,g is as low as 0.002, In the meantime, it has an additional cost of 50% CPU time compared with the D3Q19 model.The Local scheme requires less CPU time compared with the FD gradient scheme because it does not calculate the integration of the non-local variables.

Simulation of co-current two-phase Poiseuille flow
To verify the numerical accuracy of our ULBM (NMRT) PF model, a co-current two-phase Poiseuille flow is investigated.In this simulation, the computational domain in x, y and z directions is set as 10 × 100 × 1, where the periodic boundaries are set in x and z directions and the nonslip boundary in y direction.The initial density field in the y direction is given by: where the bottom part is occupied by the heavy-phase fluid and the rest by the light-phase fluid.In the simulation, the fluid is driven by a constant body force in the x direction, where F b = ρgx.Three different density ratios are tested as, ρ h /ρ l = 1000/1, 100/1 and = 10/1.In this test, we set σ = 0.001, M ϕ = 0.1 and g = 10 − 8 ; the viscosity ratioμ h /μ l = ρ h ν h /(ρ l ν l ) is fixed as 100 for all cases.It is worth nothing that, similar to the density profile, the viscosity profile is also described by a linear interpolation relation, where: It is known that for the co-current Poiseuille flow driven by the constant body force, the N-S equation can be simplified as (Fakhari et al., 2017): For comparison, we calculate the analytical solution of u x by using a second-order FD method.
The D3Q27 ULBM (NMRT) PF model with FD scheme and Local scheme is adopted in the simulation, with two different interface widths w = 2.5 and w = 5.0.Fig. 2(a) presents a comparison of the order parameter gradient achieved by the simulations and the analytical solutions based on Eq. ( 30).As shown in the figure, the FD scheme and the Local scheme both achieve reasonable results compared with the analytical solution at different values of w.Nevertheless, the Local G. Wang et al. scheme provides more accurate results for the interfacial points (also highlighted in the inserted figure of Fig. 2(a)).Fig. 2 (b ~ d) presents the comparison of simulation results and the analytical solutions for u x at different density ratios.It can be found that, for the low density ratio ρ h /ρ l = 10/1 case, results from both the FD and the Local schemes are almost consistent with the analytical solutions.However, with the increase of the density ratio, the simulation results obtained by the FD scheme clearly deviate from the analytical solutions at the large density ratios (e.g., ρ h /ρ l = 1000/1, ρ h /ρ l = 100/1).This finding is consistent with results in Zhang et al., (2022), which can be explained by the discontinuity at the interface.As indicated in Figs.2(c) and (d), when the interface is thinner (w = 2.5), the simulation results obtained by the FD scheme become worse in the case of ρ h /ρ l = 100/1 and diverge when ρ h /ρ l = 1000/1.Table 2 shows the relative errors between the simulation results and the analytical results for the cases in Fig. 2 (b ~ c).As can be observed, the relative error is increased with the density ratio for the FD scheme while it remains small for all the cases by using the Local scheme.The maximum relative error for the Local scheme is 4.07%, which occurs for the case with large density ratios (~1000).It can be seen that the local scheme improves the continuity across the phase interface and therefore leads to more accurate results and better algorithmic stability, especially for the narrow interface thickness and large density ratio.Now we give a theoretical analysis of the interface discontinuity phenomenon.According to the CE analysis of the ULBM model for the NS equations at O(ε 1 ) level (in Appendix 1), it can be found that the derivative ∇P 0 in Eq. ( 24) is used to offset the derivative caused by the additional term P 0 = c 2 s ρ in the distribution function.For example, for the m (1) g,4 , according to Eq. (A5) in Appendix 1 we have: based on Eq. (A14), we have: (32) considering the nearly incompressible limit, which implies ∂ t (P) ≈ 0, Eq. ( 31) can be re-written as: where the last term on the right-hand side of Eq. ( 33) represents the offsetting term.From this offsetting term, it can be found that the difference between ∇P 0 and c 2 s ∇ρ leads to the discontinuity at the interface.Ideally, the term c 2 s ∇ρ which is implicitly generated in the LB evolution at O(ε 1 ) level should be eliminated by the terms ∇P 0 in R g .Nevertheless, as indicated in Fig. 2.(a), when adopting the FD scheme, the errors of the order parameter gradient can be generated at the interface, thus leading to the deviations between ∇P 0 and c 2 s ∇ρ.In addition, the deviations are increased with the density ratio ρ h /ρ l (related to the ∇P 0 ) and viscosity ratio v h /v l (link to the 1/s v ).Therefore, the interface discontinuity is more significant for higher density and viscosity ratios.
To validate the above analysis, we simulate the cases in Fig. 2 in a wider range of operating parameters.The relative errors between the simulation results and the analytical solutions are listed in Table 3.The relative errors for the FD scheme increase with ν h /ν l and ρ h /ρ l and are always one magnitude higher than those for the Local scheme.For Local schemes, its accuracy is more dependant on the viscosity ratio than on the density ratio.This finding is consistent with the above analysis.It can be found that the Local scheme can provide an accurate result at an extremely high-density ratio (ρ h /ρ l = 10 6 ), which demonstrates the power of our current model in simulating the large density ratio multiphase problems.

Simulation of the Rayleigh-Taylor instability
Finally, we validate our model by simulating the RTI in twodimension (2D) and three-dimension (3D).The RTI is an instability phenomenon with strong interfacial deformation, which has a wide range of applications such as astrophysics and confinement fusion (Instability of Liquid, 1950;Dimonte et al., 2005).Firstly, we simulate a quasi-2D case by using D3Q27 ULBM (NMRT) PF model with both the Local scheme and the FD scheme.The size of the computation domain in x, y and z directions is L × 1 × 4L, where L stands for the reference length.The top and bottom walls are set as non-slip boundaries and the side walls are set as periodic boundaries.In the RTI simulation, the upper part of the simulation domain in z direction is set as the heavy phase fluid, and the fluid flow is driven by gravity with Initially, a small deformation is added on the phase interface, for the two-dimensional case, the phase-field is described as: To compare the current simulation results with the existing results in the literature, we set the values of the relevant non-dimensional numbers the same as, the Atwood number At = (ρ h + ρ l )/(ρ h − ρ l ) = 0.5, the Reynolds number Re = ρ h U 0 L/μ h = 3000, and the viscosity ratio μ h /μ L = 1, where the reference velocity is . Other parameters are also set the same: the capillary number Ca = μ h U 0 /σ = 0.26, Péclet number Pe = LU 0 /w =1000.The reference length is L = 256 and the interface thickness w = 5.We choose the reference time t 0 = ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ L/(gAt) √ as 16,000, and the dimensionless time is defined as T * = t/t 0 .The simulation results for the 2D case are shown in Fig. 3.As presented in Fig. 3(a), during the evolution, the front interface of the heavy phase falls and rolls up along the flank of the spike.With the continuous fall of the heavy phase, the rolled-up heavy phase fluid is constantly stretched and finally breaks up into small dissociative droplets.The temporal evolution of the fluid is consistent with pervious simulations, e.g., in Refs.(Ren et al., 2016;Fakhari et al., 2017;Dinesh Kumar et al., 2019).We also record the evolution of the bubble front and liquid front, and compare our simulation results with previous simulation results from Ren et al. (Ren et al., 2016) in Fig. 3(b).Both the results from the Local scheme and the FD scheme can achieve excellent agreement with previous data.
We then extend the simulation to 3D using the D3Q27 ULBM (NMRT) PF with the Local scheme.The simulation domain is set as L × L × 4L,

Table 2
The relative errors between the simulation results and analytical results.A comparison of the relative errors between the simulation results and analytical solution for the lawyer two-phase Poiseuille flow at various viscosity ratios and density ratios.and the initial phase-field is described as: For the 3D case, we set At = 0.5, Re = 128, μ h /μ L = 3, Ca = 9.1 and Pe = 744, which are the same as set by Zu&He (Zu and He, 2013).L = 128 and t 0 = 6000 are set in this simulation case.Fig. 4(a) presents a comparison between the previous simulation results (Zu and He, 2013) and our current simulation results for the time evolution of the positions of the bubble, spike and saddle (marked in Fig. 4(b)).Our simulation results agree well with Zu&He's simulation results, except for some small discrepancies at the later stage of the simulation.Fig. 4(b) demonstrates the evolution of the phase interface, which is in good agreement with the previous simulation results for the same case (He et al., 1999;Zu and He, 2013).Based on the simulation results presented in this section, the accuracy of the proposed ULBM (NMRT) PF model can be validated.In the following section, we will adopt this model to investigate the realistic multiphase problem with a large density ratio.Based on the above test results, we would like to conclude this subsection by highlighting the advantages of our proposed ULBM (NMRT) PF model: (1) The numerical stability and computational efficiency are improved by using non-orthogonal moment sets to construct MRT collision operators.(2) The proposed Local scheme resolves the interface discontinuity encountered in the conventional FD scheme.(3) The current model has the potential to simulate multiphase flows with a large density ratio, up to 10 6 .

Simulation of droplet and bubble dynamics with large density ratios and high Reynolds numbers
In this section, we demonstrate the capability of the proposed ULBM (NMRT-PFLB) model in reproducing multiphase experiments with realistic physical parameters.In subsection 3.2.1,we reproduce the experiment of a bubble rising in a water-air system.In subsection 3.2.2,we first simulate the collision of binary tetradecane droplets, and then extend the simulation to the collision of binary liquid metal droplets with extremely high density ratios.

Simulation of a bubble rising
We first consider the case of a bubble rise driven by buoyancy, which is a common phenomenon in a wide range of engineering applications.The bubble rising phenomenon has been widely investigated, but most of the existing works used small density ratios (10 1 ∼ 10 2 ) to maintain numerical stability (Wang et al., 2019).By contrast, we calculate a realistic water-air system with a density ratio of 1000 using the D3Q27 ULBM (NMRT) PF model with the Local scheme.The computational domain is a L × L × 4L tunnel, with the reference length L = 256.The bottom and top walls are set as non-slip boundaries, and the other directions are periodic boundaries.A bubble with a radius R 0 = 0.125L is placed at a distance of 0.25L from the bottom of the tunnel.The density ratio of the surrounding liquid and bubble is ρ h /ρ l = 1000, with the viscosity ratio μ h /μ l = 100, which corresponds to the water-air system.
Due to the buoyancy effect, the bubble rises from the bottom of the tunnel and eventually leads to a steady-state terminal shape and velocity (U t ) when the buoyancy force is balanced by the resistance.In this system, the terminal shape and velocity of the bubble are governed by two critical dimensionless numbers, the Bond number Bo = ρ h (2R 0 ) 2 g /σ, and Morton number Mo = μ 4 h g/(ρ h σ 3 ).To ensure the droplet terminal velocities satisfy the incompressible limitation, the gravity is set as 5 × 10 − 6 .In this section, we consider three different cases with various Bo and Mo, and the simulation parameters can be found in Table 4.
Our simulation results are compared with the experimental data in Ref. (Bhaga and Weber, 1981) and the previous simulation results from the velocity-based phase field LB model (Dinesh Kumar et al., 2019), Smoothed Particle Hydrodynamics(SPH) simulation results (Zhang et al., 2015) and front tracking method (FTM) (Hua et al., 2008) simulation results.We first quantitatively compare the droplet terminal Reynolds numbers (Re = 2ρ h U t R 0 /μ h ).As indicated in Table 4, our proposed model demonstrates better performance over other models for the large-Re cases.The maximum relative error between our simulations and experiments are 4.47%, which is lower than the simulation results obtained by the previous LBM model in Ref. (Dinesh Kumar et al., 2019).The qualitative comparison of the bubble's terminal shapes can be found in Fig. 5.As can be seen, the results of ULBM (NMRT) PF model achieve good agreement with the experimental results (Bhaga and Weber, 1981) for the bubble's terminal morphologies.Furthermore, similar to the previous simulations (Tripathi et al., 2015), at a large terminal Re, small satellite bubbles peel off from large bubbles.Thanks to the full 3D simulation with realistic physical parameters, we can capture the small fingering shapes at the edge of large bubbles, which is less obvious in previous LB simulations (Hua et al., 2008).

Simulation of the binary droplet collisions
In this subsection, we first simulate a head-on collision of tetradecane droplets using the D3Q27 ULBM (NMRT) PF model with the Local scheme and compare the simulation results with the experimental data (Tang et al., 2012).The computational domain is set as a 250 ×250 × 500 box, with periodic boundaries in all directions.Two cases are considered, namely binary equal-size and unequal-size droplet collisions.In the simulation, we set the density ratio as ρ h /ρ l = 1000.During the collisions, the droplet dynamics is governed by two dimensionless numbers, which are the Weber number We s = 2ρ h R s,0 U 2 0 /σ (R s,0 and U 0 stand for the radius of the small droplet and relative speed, respectively) and size difference Δ = R l,0 /R s,0 (R l,0 = 50 lattices stand for the larger droplet radius).Following the suggestion of Fakhari et al. (Fakhari et al., 2017), we choose M ϕ = 0.05 and w = 5 for better stability and accuracy.
For the case of equal-size droplet collision, we have We s = 40, Δ = 1, Reynolds number Re = 2ρ h U 0 R s,0 L/μ h = 500, and the viscosity ratio μ h / μ l = 166.7,which are comparable to the experiment setup (Tang et al., 2012).The qualitative comparison between the numerical and experimental results can be found in Fig. 6(a).In agreement with the experimental results, a liquid disk is created due to the coalescence and collision of the binary droplets.Then, a reflective separation of the two droplets can be observed and a liquid bridge is formed between the two separated droplets.Finally, the neck of the liquid bridge breaks and the liquid bridge retracts to form the third satellite droplet, which eventually leads to the reflective separation of the three satellite droplets.Notably, owing to the good numerical robustness, a very small satellite droplet (shown in the red braked in Fig. 6(a)) with a radius around 2 lattices (5 μm in physical unit) can be captured in our simulation.Then we simulate an unequal droplet collision, with the simulation parameter We s = 62.4,Δ = 1.78,Re = ρ h U 0 R s,0 L/μ h = 375 and the viscosity ratio μ h /μ l = 166.7.As shown in Fig. 6(b), similar to the experimental results, the droplet presents a penetration-separation process after the collision.The droplet morphologies in the simulation results are in excellent agreement with the snapshots of the experiment.
In addition to the head-on tetradecane droplet collision, we also simulate off-centre collision of binary liquid metal droplets.To satisfy the liquid metalair system, the density ratio in this simulation is set as ρ h /ρ l = 5000 and viscosity ratio μ h /μ l = 117.8.Besides, we choose We s = 61, Δ = 1 Re = 4000 and the off centre distance Δx = X/D s,0 = 0.7 (X indicates the separation between the vertical centre lines of droplets) to match the experimental conditions in Ref. (Jia et al., 2019).The comparison between the experimental results and our simulation results is shown in Fig. 6(c).It can be seen that our ULBM (NMRT) PF   Kumar et al., 2019) and the experimental results [Bhaga and Weber, 1981).
model reproduces this extremely large density ratio multiphase flow successfully, as the droplet exhibits stretching separation after collision with satellite droplets formation.The numerical droplet evolution qualitatively agrees with the experiment snapshots, except for a difference in the number of generated satellite droplets (3 for the experiment and 4 for our simulation), which is caused by the reverse breakup of the centre stretching jet.Considering this dynamic is very sensitive to the off-centre distance (Qian and Law, 1997) and the initial shape of the metal droplet, this discrepancy is acceptable.To the best of our knowledge, this is the first time that the LB method achieves this kind of extreme density ratio metal droplet dynamics with a large Re number.

Simulation of complex multiphase flow
In this section, we adopt the ULBM (NMRT) PF model to simulate the multiphase phenomena with extreme operation conditions (e.g., extreme density ratios, Weber numbers, Reynolds numbers).To test the capability of our proposed model in capturing complex multiphase instability phenomena, two experiments are reproduced, that is, the water droplet crown splashing and the liquid jet breakup.

Simulation of droplet splashing
Firstly, we simulate a water droplet splashing as it impacts a thin liquid film.In this simulation, the density ratio of the heavy phase and  et al., 2019), We s = 40, Δ = 1 and Δx = 0.7, at ρ h /ρ l = 5000 and μ h /μ l = 117.8.
G. Wang et al. light phase is set as ρ h /ρ l = 1000 with a viscosity ratio μ h /μ l = 67, which corresponds to the realistic water-air system.The same case has been studied previously using the phase field LB models, but most studies used unrealistic parameters (e.g., a small Reynolds number, usually lower than 2000 and 2-D simulation (Zhang et al., 2022;Liang et al., 2018;Fakhari et al., 2017;An et al., 2021).In our simulation, we reproduce experimental cases for a water droplet impacting a thin water film, by using the D3Q27 ULBM (NMRT) PF model with the Local scheme.The computational domain is set as a 14R 0 × 14R 0 × 6R 0 box, with periodic boundaries for the side boundaries and non-slip boundaries for bottom and top walls, the initial droplet radius is R 0 = 60.Following Fakhari et al. (Fakhari et al., 2017), we set the interface thickness w = 5 and mobility M ϕ = 0.05 to maintain numerical stability and accuracy.In the simulation, consistent with the experiment condition, we keep the droplet Weber number We = 2ρ h R 0 U 2 0 /σ = 380 and Reynolds number Re = 2ρ h U 0 R 0 /μ h = 6000 (Zhu et al., 2021).We simulate two cases with different liquid film thicknesses, h * = h /(2R 0 ) = 0.15 and h * = 0.35 (h is the liquid film thickness) and compare with the experimental cases h * = 0.07 and h * = 0.16 with the same Weber number We = 380.Notably, considering the influence of the interface thickness, the liquid film thickness in our simulation is thicker than the experiment.
Initially, the droplet is placed 10 lattices away from the liquid film, with a vertical velocity U 0 = 0.045 in lattice unit.Similar to the previous simulation, we introduce a random disturbance with a magnitude of 0.05U 0 (Ming and Jing, 2014;Liang et al., 2019) into the droplet velocity from t * = 2R 0 /U 0 = 0.1 and withdraw it when t * = 0.5.Besides the Local scheme, we also attempted the FD scheme for the same case.The simulation diverges as the droplet contacts with the liquid film with the FD scheme, indicating Local scheme have a better numerical stability in the practice.The comparison between the numerical and experimental results is shown in Fig. 7.In our simulation, We = 380 0044 , which corresponds to the breakup conditions proposed in Ref. (Burzynski et al., 2020).As shown in Fig. 7(a) and Fig. 7(b), similar to the experiment observation (Zhu et al., 2021), both cases capture the droplet prompt splash (shown in the red bucket) which occurs a short period after the droplet contact with the liquid film.The radiuses of the satellite droplets in the prompt splash range from 1~3 lattice (1.7% ~ 5% of the R 0 ), which are in good agreement with the results of the experiments in Ref. (Burzynski et al., 2020).
Besides the prompt splash, the simulation results also reproduce the key dynamics features observed in the experiment.As shown in Fig. 7 (a), when the droplet impacts a thicker liquid film, rim instability can be observed at the edge of the liquid crown.Liquid fingers can be found as the crown is retracted while there is no breakup of liquid fingers, which is consistent with the experimental observation (Zhu et al., 2021).We also observe the central jet at the later stage of the droplet evolution, which is formed by the rebound of the central liquid sheet.This phenomenon is also observed in the experiment.For the case of droplets impinging on the thinner liquid film, we can find liquid fingers generated with the liquid crown developing, and secondary droplets escape from the liquid fingers as the liquid crown still expands.Moreover, we observe the fingering jet formed from the retraction of the liquid crown and eventually breakup into satellite droplets.These numerical results are consistent with the experimental phenomena (Wu et al., 2021;Lu et al., 2020).It should be noticed that the fingering breakup of the liquid crown is less obvious in previous 3D LB simulations (Sitompul and Aoki, 2019;Schwarzmeier et al., 2023), which proves our new model has a better ability to reproduce the experimental phenomenon.
In additional to the qualitative comparison of the droplet morphologies, we also quantitatively compare the evolution of the droplet crown diameter with the experimental data (Zhu et al., 2021) and the power law fitting function r/D 0 = C(Ut/D 0 ) 0.5 (Fakhari et al., 2017).As shown in Fig. 7(c), our simulation results agree well with the experiment, for both the diameter evolution and breakup point (marked in the figure).By fitting simulation results based on the power law fitting functions, the fitting prefactor C is obtained as 0.8 and 1 for the cases of h * = 0.07 and h * = 0.16, respectively.The fitting values are close to previous experiment findings (Yarin and Weiss, 1995) and the famous Wagner's theorem (Wagner, 1932).

Fig. 7.
Comparison between the experimental snapshots (grey) [Zhu et al., 2021) and the simulation results (colour, coloured by the normalized velocity) of a water droplet impacting a thin water film at We = 380 and Re =6000.The dimensionless liquid film thickness (a) h * = 0.16 and (b) h * = 0.07; (c) the time evolution of droplet crown diameters and a comparison with the experiment data (symbols) and power law fitting (lines).

Simulation of liquid jet breakup
In this subsection, we simulate the breakup of a liquid jet when injected into another liquid.In this simulation, the density ratio of the heavy phase is set as ρ h /ρ l = 1.7 and viscosity ratio μ h /μ l = 7.14 (ν h /ν l = 1/4.2).These setups correspond to the simulation setup (Saito et al., 2017) and experimental conditions (Saito et al., 2016) of Saito et al.'s works.The computational domain is a cuboid box, and the top wall is set as the non-slip boundary of a circular inlet with a diameter d 0 at the central, where the heavy phase is injected with an initial speed u 0 .The bottom wall is set up as a Neumann outflow boundary, and the side walls are set as periodic boundaries.Following the setup of Saito et al. (Saito et al., 2017), in this case, the body force is written as which indicates the gravity only acts on the heavy phase.The interface thickness is set as w = 5 and mobility M ϕ = 0.05.The dynamics of the injected liquid jet is governed by three critical dimensionless numbers, ).Three typical jet dynamics are simulated in the present work, including dripping, sinuous and atomization.Firstly, we consider the situation of the dripping jet.In this simulation, we set d 0 = 20 and the computational domain is 140 × 140 × 350.The characteristic numbers are Fr = 8.5, We = 5 and Re = 300, corresponding to Oh = 0.0075.The evolution of the liquid jet is shown in Fig. 8, in which we can find a swollen part generated at T * = tu 0 / d 0 = 7.5 and pinch-off from the liquid jet at T * = 30.The jet tip performs a long-wave mode primary breakup, and the short waves generated after the pinch-off of the primary droplet propagate upstream and lead to a later short-wave mode breakup.These qualitative findings are consistent with previous simulation results for the small Oh liquid jet breakup mechanism (Shinjo and Umemura, 2010).Notably, benefiting from the high accuracy of the ULBM (NMRT) model, we reproduce the tiny satellite droplets (pointed out in the figure) when the primary droplet breaks up.This phenomenon is also observed in the previous experiments of the same liquid-liquid system (An et al., 2021;Saito et al., 2017), while it is hardly observed in previous simulations by using the colour gradient LBM (An et al., 2021).
We then simulate the liquid jet breakup in the sinuous regime, with Fr = 8.5, We = 220 and Re = 2200, corresponding to Oh = 0.0075.As shown in Fig. 9, a mushroom-shape head appears after the heavy-phase liquid injection.With the development of the liquid jet, similar to the experimental observations (Saito et al., 2017) and previous simulation (An et al., 2021), we can find the return flow of a mushroom head causes Kelvin-Helmholtz instability waves (i.e., T * = 15) on the liquid column.Afterwards, the pinch-off of the instability waves causes droplet entrainment (i.e., T * = 28).Since there is no artificial perturbation, the instability waves are caused naturally owing to the nonlinear hydrodynamics.In addition, we quantitatively compare the liquid jet penetration length with the experimental and numerical results.As shown in Fig. 11(a), our simulation shows an excellent agreement with Saito et al.'s experiments (Saito et al., 2016).Both our current simulation results and their latest simulation by using the D3Q27 MRT colour gradient LB model (Saito et al., 2017) present better accuracy than their previous simulation by using the D3Q19 MRT colour gradient LB model (Saito et al., 2016).
Finally, we simulate the jet atomization at large We and Re.In this simulation, Fr = 8.5, We = 10 5 and Re = 3400, corresponding to Oh = 0.03.It is worth mentioning that the experimental viscosity ratio (ν h /ν l = 1/4.2) (Saito et al., 2016) is adopted in this simulation, which leads to the kinematic viscosity of the light phase being as lower as υ l = 3.5 × 10 − 4 .As indicated in the Fig. 10, similar to the sinuous case, a mushroom shape head can be observed after the liquid jet injection (i.e., T * = 7.5).On the other hand, due to the strong Kelvin-Helmholtz instability, we can find the mushroom-shape head breaks up soon after the jet penetration and tiny droplets are generated by the pinch-off of the instability waves (as shown when T * = 15).With the continuous development of the liquid jet, numerous tiny droplets are detached from the liquid column.The qualitative evolution of the liquid jet also agrees well with previous simulations (Saito et al., 2017;Saito et al., 2016).
We conclude by plotting all of the simulated cases in a regime map (Saito et al., 2017;Saito et al., 2017) and qualitatively comparing the liquid jet morphologies with the experimental results in Ref. (Saito et al., 2017).As presented in Fig. 11(b), our simulation results are consistent with the experiment snapshots, and our simulation parameters match the corresponding regimes.The simulation results in section 3.3 demonstrate the excellent stability and accuracy of the ULBM (NMRT) PF model since our simulations not only reproduce the key instability phenomena involved in various complex multiphase experiments but also provide accurate and stable results with the demanding physical conditions (extreme lower viscosities, large density ratios, high Weber numbers).The results in this section prove the power of the proposed model in modelling the complex multiphase problem.

Conclusion
In this paper, we propose a three-dimensional phase-field lattice Boltzmann model (PF LBM) for multiphase flows at large density ratios and high Reynolds numbers.The proposed model is built within a recently proposed unified lattice Boltzmann model (ULBM) framework, with a non-orthogonal moment set being used to construct the multiplerelaxation-time (NMRT) collision operator.In addition, we provide a local gradient calculation scheme which is based on the non-equilibrium part of the distribution function in raw moment space.Benefiting from the non-orthogonal moment set, our proposed model leads to a simplified implementation and enhanced computational efficiency.According to a detailed Chapman-Enskog analysis, our proposed model can accurately recover the target conservative Allen-Cahn equation for interface tracking and the target Navier-Stokes equations in the nearly incompressible limit.
The proposed model is comprehensively validated and assessed by comparing our simulation results with experimental results, analytical results, and previous simulation results.It is shown that the constructed ULBM(NMRT) PF model can provide better numerical stability compared with the single-relaxation-time collision operator.Besides, the local gradient calculation scheme is found to reduce the computational cost and provide better numerical accuracy compared with the finite difference scheme for gradient term calculation.We also find that the Local scheme prevents the interface discontinuity that occurs in the simulation of co-current Poiseuille flow, which can be explained by the fact that the derivative terms calculated by the Local scheme in the discrete forcing eliminate the derivative terms generated by the modified distribution function.
After validations, the proposed model is adopted to simulate the Rayleigh-Taylor instability and the realistic multiphase problems with large density ratios, i.e., binary droplet collision and bubble rising.The current simulation results achieve good agreement with previous simulation results and experiment results, both qualitatively and quantitatively.Notably, our simulation accurately reproduces the binary collision of liquid metal droplets at an extremely high density ratio (up  Finally, the proposed model is adopted to simulate challenging multiphase problems under extreme conditions, including droplet splash on impacting liquid film and liquid jet spray.For the case of droplet splash, the simulated Reynolds number reaches 6000, at a density ratio of 1000.For the case of a liquid jet spray, our simulation can reproduce various liquid jet breakup regimes for a variety of operating parameters, with the Weber number reaching 10 4 and the viscosity as low as ∼ 10 − 4 in lattice unit.The simulation results agreed well with the previous experimental results in both cases, demonstrating unprecedented capabilities of the proposed NMRT-PFLB model within the ULBM framework.Considering the generality and versatility of the ULBM framework, the new NMRT-PFLB model can be used to tackle a wide range of realistic multiphase problems, such as jet spray in engines and boiling in powerplants.In addition, the ULBM framework is able to incorporate any new and advanced multiphase flow models to further extend its capability. where , e 1a ……, e 18a )]M − 1 for a = x, y, z. by introducing the following relation: where ε stand for an order parameter, Eq. ( A2) can be written in the consecutive orders of ε: Firstly, we provide a CE analysis for the ULBM (NMRT) PF model in Section 2.2 for the AC equation.For the distribution function f ϕ,i , it has Consequently, according to Eqs. ( A24) and (A25) and assuming m (1) Substituting Eqs.(A25) and (A26) into Eq.(A5), we have the relation at ε 1 level for i = 0: It can be found that the Local gradient calculation scheme Eq. ( 25), named as Local scheme) is achieved by substituting the relation m (1) ϕ = m ϕ − m eq ϕ into Eq.(A8).Liang et al. (Liang et al., 2018) proved that the Local scheme achieves second-order accuracy in space.By Eq. (A6), we can get the following relation at ε 2 level for i = 0: Substituting Eq. (A26) and Eq.(A8) into Eq.(A9), we can achieve the following equation at O(ε 2 ) level: where s ϕ,1 = s ϕ,2 = s ϕ,3 = 1/(M ϕ /c 2 s +0.5).Combining Eq. (A7) at O(ε 1 ) and Eq.(A10) at O(ε 2 ) level and substituting Eq.( 11) for F ϕ , we can get: ∂ϕ ∂t which is the target AC equation (Eq.( 1)).Therefore, it shows that the current ULBM (NMRT) model in Section 2.2 can recover the target interface tracking equation without high order error terms involved.Then, we provide a CE analysis for the ULBM (NMRT) model for the incompressible NS equation.Based on Eqs.(A24, A27) and (A28), and assuming m (1) g = m g − m eq g , m (n) g = 0 (n ≥ 2), we have : (A12) substituting Eqs.(A27, A28) and (A12) into Eq.(A5), for i = 0 ∼ 3, the following continuity and momentum equations at ε 1 level can be achieved: where i = 0…26 and |⋅〉 donates a 27-column vector.The superscript T is the transposition symbol.The relaxation matrix S is: For the D3Q27 lattice model, the non-orthogonal raw moment set is designed as: Where the subscript of moments is in ascending order of (p + q + n).Accordingly, the transformation matrix is: Following the introductions in Section 2.2, the discrete equilibrium moment and forcing terms in the raw moment space for D3Q27 lattice model can be written as: m eq ϕ = Mf eq ϕ,i = ⎡ ⎣ ϕ, ϕu x , ϕu y , ϕu z , 0, 0, 0, ϕ, 0, 0, ϕu x c In solving the NS equation, similarly, the discrete equilibrium moment and forcing terms in the raw moment space can be written as:
(i.e., σ = 0.1, 0.05, 0.01).The pressure difference ΔP between different phases and the analytical results based on the Laplace law ΔP = 2σ /R 0 are plotted in Fig. 1(a).The numerical results show good agreement with the analytical results, with the maximum relative error being less than 3%.Fig. 1(b) demonstrates the density profiles for the case of R 0 = 30, with two different interface widths, w = 3 and w = 5.For comparison, we also plot the numerical results of the FD scheme and the analytical solutions obtained by Eq. (27).As indicated in Fig. 1(b), both the Local scheme and the FD scheme achieve excellent agreements with the analytical results.Fig. 1(c) presents the relative error (ϵ = |ρ sim.−

Fig. 1 .
Fig. 1.(a) Laplace law verification for a static droplet, where the different symbols stand for different surface tension, the dashed lines stand for the analytical results; (b) the density profiles of the droplet interface at the density ratio of 1000 (ρ h = 1000, ρ l = 1), obtained by using FD scheme and Local scheme, the lines in the figure represent the analytical solution; (c) relative errors between simulation results and analytical solutions of density distribution.
et al. maximum spurious velocities are similar between the FD and the Local schemes, and the spurious velocities of the MRT model are always one order of magnitude lower than the SRT model.It is worth mentioning

Fig. 2 .
Fig. 2. Simulation of a two-phase co-current Poiseuille flow at viscosity ratio μ h /μ l = 100; (a) the gradient of the order parameter in the y direction; the x direction velocity at the density ratio equals (b)ρ h /ρ l = 10, (c) ρ h /ρ l = 100, and (d) ρ h /ρ l = 1000; the simulation results in the figure are obtained by FD ULBM (NMRT) PF (blue) and Local ULBM (NMRT) PF (red) models, and the lines in the figure stand for analytical results.

Table 4
Comparison of simulated Re with the experiment result and previous simulations.