A Momentum-Conserving Weakly Compressible Navier–Stokes Solver for Simulation of Violent Two-Phase Flows with High Density Ratio

A consistent and conservative formulation for mass and momentum transport is proposed in the context of simulating incompressible two-phase flows by using weakly compressible method. Combined with the evolving pressure projection method to prevent oscillation of the solution induced by the acoustic wave, this solver aims at a robust and accurate computation of violent two-phase flows with a high density ratio, while taking advantage of fully explicit time integration of the weakly compressible Navier–Stokes equations. Coupled with the volume of fluid method for capturing interfaces, the mass and momentum fluxes are evaluated in a consistent manner using the finite volume method. In addition, a special implementation of the pressure projection is devised to avoid velocity-pressure decoupling on a collocated grid. The solver's accuracy and stability are demonstrated through various two-phase flow simulations, including dam break and liquid jet atomization scenarios, emphasizing its momentum-conserving properties.


Introduction
Among the various numerical modelings and simulations of two-phase flows, the direct approach that solves the macroscopic Navier-Stokes equations in combination with a two-phase interface capturing method is widely used in fundamental research and industrial applications (Mirjalili, Jain, and Dodd 2017;Zaleski and Xiao 2020). The one-fluid model framework allows for the straightforward extension of numerical methods initially designed for single-phase flow to two-phase flow computation. However, significant challenges emerge when dealing with high density ratio between two phases, where the order of 10 3 is the upper limit for many existing solvers, particularly those utilising the non-conservative form of Navier-Stokes equations. The situation deteriorates when the motion of the interface is violent, as steep changes in density and momentum across the interface result in pressure gradient discontinuity, which can pose problems for unspecialised numerical discretisation.
CONTACT Kai Yang yang@sim.gsic.titech.ac.jp To achieve a robust and accurate simulation of violent two-phase flows with a high density ratio, a consistent formulation of mass and momentum transport is essential (Rudman 1998). While a staggered grid system is commonly used for the spatial discretisation of governing equations, maintaining consistency becomes particularly complex due to the different locations of mass and momentum control volumes. A classical solution is to compute the volume fraction of fluids on a mesh twice as fine as the velocity-pressure staggered mesh. This makes the mass flux through a small control volume directly available for calculating the momentum flux through a large control volume, even if they do not fully coincide. To reduce the computational overhead of solving an advection equation for the volume fraction on sub-cells, some approaches manage to construct auxiliary volume fractions in respective face-centred control volumes of the momentum, either by summing the PLICreconstructed values in the sub-cells (Zuzio et al. 2020) or by shifting the original volume fractions of two adjacent cells sharing the face by using Weymouth and Yue's advection scheme (Arrufat et al. 2021). While momentum transport can be made consistent with the auxiliary mass, it is not conservative because of the discrepancy in the reconstructed interface. Instead of using the volume fraction equation, another method synchronises the face-centred density with the level-set function after linearly interpolating the level-set function to the face (Nangia et al. 2019). Recently, a phase-field model based transport scheme is developed by algebraically modifying the momentum equation taking into consideration the artificial flux term in the phase-field equation (Mirjalili and Mani 2021). A different realisation of consistent transport is presented in Patel and Natarajan (2017) where a hybrid staggered/non-staggered grid is employed, and the control volume of the momentum actually occupies two cells associated with the mass. And an identical discrete operator is required for convection of both mass and momentum.
In contrast, a collocated grid is inherently suitable for evaluating convective fluxes of momentum consistently with the mass, because both the density and velocities are defined at the cell centre. Since the pioneering work of Bussmann,Kothe, and Sicilian (2002), new numerical methods for two-phase flow simulations on a collocated grid, where the consistent and conservative momentum transport can be simply implemented, have been continuously emerging (Popinet 2009;Balcázar et al. 2014;Fuster and Popinet 2018;O'Brien and Bussmann 2020). Unlike in the case of a staggered grid, a crucial consideration for ensuring a stable computation on a collocated grid is how to avoid velocity-pressure decoupling. Various techniques, such as the projection filter (Rider 1998) and double projection (Aprovitola and Denaro 2007), have been proposed to enhance the coupling and eliminate numerical modes in the solution. A comparative study has also been reported on the stability and accuracy of pressure projection methods on the collocated grid (Abbasi, Ashrafizadeh, and Shadaram 2013).
The two-phase flow solvers discussed above rely primarily on the incompressible Navier-Stokes equations, which require an implicit solution to a pressure Poisson equation. However, iterative methods such as Jacobi, Gauss-Seidel and SOR are known for their slow convergence rates when applied to two-phase flows with a high density contrast. To address this challenge, Krylov methods and multigrid preconditioners are used to efficiently solve the Poisson equation with variable coefficients (Nangia et al. 2019;El Ouafa, Vincent, and Le Chenadec 2021). Despite their effectiveness, the complexity of these algorithms presents significant implementation challenges, particularly on largescale multi-node clusters (Zaspel and Griebel 2013;Ha,Park, and You 2018;. Alternatively, weakly compressible Navier-Stokes equations are gaining popularity for simulating incompressible fluid flows. Under low Mach number and isothermal conditions, an independent pressure evolution equation can be derived (Toutant 2017) and explicitly integrated over time. And its spatial discretisation depends only on a local and compact stencil. The hyperbolic-parabolic system of weakly compressible Navier-Stokes equations offers excellent performance and scalability for large-scale computing on multinode supercomputers Yang and Aoki 2023). The simple numerical scheme used for calculating the pressure evolution equation is flexible and can be readily applied to both unstructured and adaptive Cartesian grids, as demonstrated in Matsushita and Aoki (2021). This stands in contrast to multigridbased Poisson solvers, which may be less versatile in their implementation. A comparison of various pressure evolution equations has been conducted to verify their accuracy (Toutant 2018).
Although weakly compressible Navier-Stokes equations have been successfully applied to two-phase flow simulations (Matsushita and Aoki 2019), several major issues persist. One is the oscillation of the velocity and pressure induced by the acoustic wave in weakly compressible flows. This phenomenon has also been observed in simulations based on the lattice Boltzmann method (Inamuro et al. 2016;Sitompul and Aoki 2019), which is a well-known weakly compressible method (Lallemand et al. 2021). We previously proposed an evolving pressure projection method to mitigate compressibility effects and enhance the accuracy of two-phase flow simulations (Yang and Aoki 2021). However, there has been little research aimed at achieving momentum conservation and consistent transport in weakly compressible method for simulating incompressible two-phase flows. A recent study based on the entropically damped artificial compressibility (EDAC) model solves the conservativeform momentum equation in combination with an artificial mass-conservation equation derived from the phase-field model (Kajzer and Pozorski 2020).
Nevertheless, the two-phase flow simulations in that study are restricted to moderate density ratios and low Reynolds numbers.
The objective of this study is to develop a weakly compressible Navier-Stokes solver that features a consistent and conservative formulation for momentum transport, an aspect that was lacking in our previous work (Yang and Aoki 2021) and, to the best of our knowledge, in the entire research field of weakly compressible methods for simulating incompressible two-phase flows. We also propose a novel implementation of our evolving pressure projection method to dampen acoustic waves while maintaining the coupling between pressure and velocities on the collocated grid. The proposed solver aims to provide a robust and accurate computation for violent two-phase flows with high density ratios, leveraging the hyperbolic-parabolic characteristics of weakly compressible Navier-Stokes equations for fully explicit time integration.
The remainder of this paper is organised as follows. In Section 2, we introduce the equations governing the dynamics of weakly compressible two-phase flows. In Section 3, we propose numerical methods for solving the governing equations in a conservative and coupled manner. Section 4 presents simulation results to validate the momentum-conserving weakly compressible solver and demonstrate its benefits. Finally, we provide concluding remarks in Section 5.

Mathematical Model
The Navier-Stokes equations for weakly compressible fluid under isothermal condition can be written as: where u = (u, v, w) is the velocity vector, ρ the density, p the pressure, μ the dynamic viscous coefficient and F the external body force. c s is an artificial speed of sound, which must be large enough to satisfy the low Mach number condition Ma = |u| max /c s 1. The conservative form of the convection term in the momentum equation is crucial in the present model.
We simulate an immiscible two-phase fluid system as a one-fluid model. An indicator function φ, which represents the volume fraction of the heavy fluid, is used to distinguish the different fluid regions in the domain and the motion of the interface between two fluids is governed by the volume of fluid (VOF) equation (Hirt and Nichols 1981): Note that the right-hand side of Equation (3) accounts for the change in volume fraction due to fluid compressibility, and it vanishes when the fluid flow is incompressible. We neglect the right-hand side term for two reasons. One is that compressibility is small enough at low Mach number to consider the fluid as being nearly incompressible. In addition, since the continuity equation is not explicitly included in this model, the VOF equation without the right-hand side also represents conservation of mass. Fluid properties such as density and the viscosity coefficient can be evaluated on the basis of the volume fraction function, where the physical properties of the heavy and light fluids are indicated by the subscripts h and l, respectively.
The surface tension force in the two-phase flow is formulated by a density-scaled continuum surface force (CSF) model (Yokoi 2014), F s has the following form: where σ is the coefficient of surface tension, and κ is the curvature of the interface.

Consistent Transport of Mass and Momentum
A collocated grid system is used for spatial discretisation, which means all primitive variables are defined at the centre of each cell. The finite volume method (FVM) is used to ensure numerical conservation of mass and momentum when calculating the convective parts of Equations (1) and (3). In the present method, the mass transport is first calculated by solving the VOF equation. The THINC/WLIC scheme (Yokoi 2007), which has both good accuracy and simplicity of implementation, is employed to algebraically evaluate the numerical flux of the volume fraction function. The THINC (tangent of hyperbola for interface capturing) scheme (Xiao, Honma, and Kono 2005) uses a piecewise hyperbolic tangent function to represent the jump of volume fraction across the interface in one dimension: where α = 1 for φ i−1 < φ i+1 and α = −1 for φ i−1 > φ i+1 , the parameter β which controls the steepness and thickness of the interface is chosen to be 2.5 here. x i represents the distance from the jump centre of the interface, and it is determined by solving To extend the THINC scheme to multiple dimensions, the WLIC (weighted line interface calculation) method reconstructs the interface by using a weighted average of one-dimensional line interfaces along each axis of space. The weights depend on the interface normal. When calculating the flux of the volume fraction F V in a certain axis direction, only the line interface along this axis is reconstructed by the THINC scheme and its flux is a spatial integration of Equation (7). The fluxes of the line interfaces perpendicular to this axis can be simply evaluated by using the 1st-order upwind scheme. For example, the flux in the x-direction can be expressed as where iup = i for u i+ 1 2 > 0, and iup = i + 1 for u i+ 1 2 ≤ 0. The weight of the line interface along the x-axis is calculated using the interface normal n = ∇φ/|∇φ|: where n x , n y and n z are three components of the normal vector in the x, y and z directions respectively. As illustrated in Figure  crossing the cell face during t are related as follows: where the subscript f denotes an arbitrary cell face. u f is a face-centred advection velocity, and its construction will be described in the next subsection. The total mass flux F Mass is obtained by summing the volume fluxes of the two-phase fluids multiplied by their corresponding densities: Finally, the convection term of the momentum equation is calculated using Gauss's divergence theorem: where A f represents the area of the cell face. U f is the velocity at the cell face that is reconstructed from the cell-centred velocities by a 3rd-order MUSCL scheme (Van Leer 1979). Since the volume flux in the VOF equation is included in the momentum flux, the mass and momentum transports are coupled by the above consistent formulation. Only by this way, the solution of velocity is correct when extracting the velocity from the momentum after solving the transport equations, For comparison, we provide an inconsistent formulation. The convective part of the momentum equation is replaced by a non-conservative form, and solved by the finite difference method with the 3rd-order WENO scheme (Jiang and Peng 2000).
Since the time integration is directly implemented on the velocity rather than the momentum, there is no need to worry about the abnormal solution caused by extracting the velocity from the momentum via Equation (13). In this approach, the mass transport and momentum transport are independent at each time step.

Evolving Pressure Projection Method
Following Chorin's projection method for solving the incompressible Navier-Stokes equations, an intermediate velocity u * is updated from u n at time step n by solving the momentum equation, Equation (1), without the pressure gradient term. The intermediate velocity needs to be corrected by a pressure projection step, which gives the solution at the next time step For weakly compressible flow simulations, we have proposed an evolving pressure projection method (Yang and Aoki 2021) to alleviate the undesirable effect of the acoustic wave and consequently diminish the fluctuations in the velocity and pressure fields. First, the pressure can be easily updated from p n using Equation (2): By substituting p n+1 in Equation (15) with p * ,0 , the intermediate velocity u * is corrected to u * * ,0 . Next, the following iteration process is conducted for N steps: (i) compute the increment of pressure δp and update the temporary pressure p * ,i to p * ,i+1 , (ii) update the velocity field u * * ,i to u * * ,i+1 by projecting the increment of pressure, where the superscript i indicates the current step of the iteration and starts from 0. The effectiveness and computational efficiency of different number of iterations have been studied in our previous work (Yang and Aoki 2021). Considering the severe restriction of the time step imposed by the speed of sound, there is a trade-off between the Mach number and the number of iterations in order to accelerate the time-to-solution.
Here, we find that 20 iterations in every time step are sufficient if the Mach number is smaller than 0.05, which can be realised by dynamically adjusting the speed of sound as the simulation goes on.
The above method with two iterative steps can be simplified to a new pressure evolution equation as follows: This hyperbolic equation for pressure can be explicitly advanced in time by using a local spatial stencil. It significantly benefits scalability of computation in spite of the speed-of-sound restriction on the time step. It is worth noting that the Equation (20) reduces to the Poisson equation when the speed of sound is infinite or when the evolution of pressure converges with time. This means that in the limit case, the proposed method will converge to the incompressible fluid flow.
To avoid pressure-velocity decoupling when the pressure projection is performed on the collocated grid, we follow the concept of approximate projection (Popinet 2003). Imitating the staggered grid system, a face-centred velocity field u f is constructed by linear interpolation of the corresponding velocity components defined at the cell centre, as illustrated in Figure 2. The external forces coming from the surface tension and gravity are first applied to the face-centred velocities, where u * c− and u * c+ are the intermediate cell-centred velocities at either side of the face. The density ρ f on the cell face is also linearly interpolated from the densities at the cell centres. Then, the staggered pressure and velocities are iteratively corrected by using the evolving pressure projection method. Finally, the cell-centred velocity u c is corrected by averaging the pressure gradient and body forces from the cell faces to the cell centre: Additionally, the corrected face-centred velocity will be used as the advection velocity, which also appears in Equation (10), in the next time step.

Time Integration
We advance the numerical solutions from time step t n to t n+1 = t n + t by using the fractional-step method, which proceeds as follows: (1) Compute the consistent transport of mass and momentum by simultaneously solving the convective parts of the VOF equation and the momentum equation based on Equations (7) - (12): We use a three-stage third-order strongstability-preserving Runge-Kutta scheme (Gottlieb, Shu, and Tadmor 2001) for the explicit time integration of the transport equations.
(2) Compute the average density and viscosity with the updated VOF function by Equations (4) and (5).
(3) Explicitly solve the viscous part of Equation (1) by using a first-order Euler forward scheme for the time integration: (4) Construct the face-centred velocity and apply the external forces to the cell face by using Equation (21). (5) Iteratively correct the face-centred velocity and pressure by using the evolving pressure projection method, i.e. Equation (16)-(19). (6) Approximate projection and forcing for the cell-centred velocity by using Equation (22).
The algorithm described above can be easily switched to an inconsistent formulation by replacing Equation (23b) with Equation (14). In this case, we still use the Runge-Kutta scheme to integrate the nonconservative momentum convection equation in time, although it is no longer coupled with the transport of the volume fraction, Equation (23a).
The time step t in our weakly compressible Navier-Stokes solver for two-phase flow simulations is constrained by advection, propagation of the acoustic wave, viscous stress and surface tension. Thus, the following conditions should be satisfied (Matsushita and Aoki 2019): where N d is the number of spatial dimensions. Safe coefficients CFL adv = CFL acs = 0.4, s visc = s sf = 0.1 are used in the simulations described below.
In our computation, we dynamically adjust the time step based on changes in flow velocities. As the Mach number is fixed at 0.05, we adjust the artificial speed of sound accordingly to 20 times the maximum flow velocity. Since advection dominates the following twophase flow simulations, the time step size is restricted by the speed of sound. Consequently, our time step is limited to one-twentieth that of an incompressible solver. Despite the limitations on time step size, the weakly compressible flow solver still outperforms an incompressible solver that employs an SOR method to solve the Poisson equation in terms of the time to solution, as demonstrated in our previous study (Yang and Aoki 2021).

Single Vortex
The first benchmark test is the single vortex of Rider and Kothe (1998) for evaluating the accuracy of the interface capturing method. A circular disk of radius 0.15 is centred at (0.5, 0.75) in a unit square domain. The time-dependent velocity field for advection is given by u = 2 sin 2 (π x) sin πy cos πy cos πt T , (26a) v = −2 sin (πx) cos (πx) sin 2 π y cos πt T , where the period T = 8. The circle is stretched to a thin and long filament at t = 0.5T, then it returns to its original location at t = T. To quantify the numerical errors of the THINC/WLIC method used in this paper, the L1 error is defined as where the exact solution φ exact is computed by Equation (7). And the final solution of the VOF function at t = T is used to evaluate the accuracy. The evolutions of interface in single vortex flow on 64 × 64, 128 × 128, 256 × 256 and 512 × 512 are depicted in Figure (3) at t = 0.5T and Figure (4) at t = T, respectively. As the mesh resolution increases, the thin tail and final shape can be well resolved. A reasonable convergence order of the L1 error is also shown in the Table 1.

Two-Dimensional Rising Bubble
For the validation of two-phase flow, we first solve a non-violent problem in order to compare the numerical results with other benchmark data. The accuracy of the present method is evaluated in a two-dimensional rising bubble simulation. At the beginning, a circular bubble with diameter 0.5 is surrounded by still liquid in a 1 × 2 rectangular domain. The centre of the bubble is located at (0.5, 0.5). No-slip conditions are imposed on the top and bottom enclosing walls, and slip conditions are imposed on the left and right sides. The computation is performed on a 128 × 256 uniform grid. The physical properties of fluids are initialised by the two groups of parameters in Table 2 for two cases respectively. The evolution of the shape of the bubble until the ending time t = 3 is depicted in Figure 5. Since the present method is mass-conservative, and accurate enough for interface reconstruction, the bubble trail in the case 2 is well captured in our results, which is difficult for the level-set-based method (Hysing et al. 2009). Next, we compare the rising velocity of the bubble with the results of Hysing et al. (2009) and Aland and Voigt (2012). The latter employs the finite element method for solving the Navier-Stokes and diffusive-interface phase-field system. The rising velocity plotted in Figure 6 shows good agreement in both cases. Different from our previous study, the iteration number is fixed and the speed of sound is adjusted dynamically, as described in Section 3.2. There is almost no oscillation in the velocity, which indicates that the proposed approach successfully minimises the effects of weak compressibility and the acoustic wave. Moreover, as both references solve the incompressible Navier-Stokes equation, these comparative results prove that our solver can be used to accurately simulate an incompressible two-phase flow.

Transport of a Heavy Droplet
The transport of a heavy droplet is a test that is often conducted in studies on the momentum-conserving or consistent transport methods. Here, a circular droplet of diameter D = 0.25 is placed at the centre of a periodic unit square domain. The density of the liquid droplet is ρ h = 10 6 , and the density of the surrounding gas is ρ l = 1. Both the droplet and the gas are inviscid, i.e. μ h = μ l = 0. This problem is free of the gravity and surface tension force. In a disk area of diameter D + 2 x concentric with the droplet, the velocities are initialised to u = 1, v = 1, while the fluid in the other region is at rest. The droplet velocity is initialised in an area slightly larger than the droplet in order to improve the results in the early stage of the simulation (Arrufat et al. 2021). The computation is performed on a 128 × 128 uniform grid, correspondingly, 32 meshes are assigned to the diameter of the droplet.
The droplet goes through a cycle of motion at t = 1 and returns to its original position. Figure 7 compares the droplet shapes obtained by the consistent and inconsistent methods described in Section 3.1. Although the density ratio 10 6 is very high and this problem is dominated by convection of the heavy fluid, the droplet undergoes unphysical distortion if the mass and momentum transports are inconsistent. In contrast, the consistent method retains the circular shape. There is no significant deformation of the droplet. Hence, it is necessary to use the consistent method when the density ratio between the two phases is high. Or rather, the consistent method is always preferable for the two-phase simulation.

Drop Oscillation with Surface Tension Force
Next, we investigate a drop oscillation problem driven by surface tension force. This problem is widely used to verify surface force models. However, here we focus more on the loss of kinetic energy during the oscillation. We use the same computational setting as in Torres and Brackbill (2000), Xie et al. (2020), in which a 20 2 square region is discretised by a 128 × 128 uniform grid. The initial drop interface is an ellipse defined by the equation, (x − 10) 2 /9 + y − 10 2 /4 = 1.
The density and viscosity of the drop are 1 and 0.01 respectively, and the fluid outside the drop has density  0.01 and viscosity 5 × 10 −5 . This problem is gravityfree, only surface tension force is applied. The surface tension coefficient is σ = 1. At the beginning, the fluid in the whole domain is at rest. Slip conditions are imposed on the enclosing walls.
The evolution of the total kinetic energy 1 2 ρ|u| 2 dV of the computational domain is plotted in Figure 8. The frequencies of the oscillation obtained by the consistent and inconsistent methods are basically the same and are close to the values reported in the references   (Torres and Brackbill 2000;Xie et al. 2020), which indicates that the accuracy of the surface tension force in the present solver is sufficient. However, in the drop oscillation process, loss of fluid momentum in the inconsistent method is much greater than that in the consistent method. This shows the significant advantages of the consistent method by stably solving the conservative-form momentum equation. Not only is the momentum conserved during the transport process, but the consistency with the mass transport also prevents non-physical transfer of momentum from the heavy phase to the light phase.

Collapse of Water Column
The two-dimensional collapse of a water column, which is also known as dam break, is simulated to demonstrate the importance of consistent transport on the stability of the two-phase interface at high Reynolds number. We initialise the computation by following the experimental setup of Lobovskỳ et al. (2014). A reasonable simplification to a twodimensional simulation is made based on the fact that the water surface profile is almost the same in the breadth direction of the tank during the collapse process. The initial water column has a width 0.6 m and height H = 0.3 m, located at the right and bottom of a rectangular tank of size 1.61m × 0.6m. A schematic representation of the water and domain is illustrated in Figure 9. On the downstream vertical wall, we select four pressure monitoring points at the same height as in the experiment: 3, 15, 30 and 80 mm. The physical properties of the water and air are ρ h = 997 kg/m 3 , ρ l = 1.2 kg/m 3 , μ h = 8.8733 × 10 −4 Pa·s, μ l = 1.8 × 10 −5 Pa s, and surface tension coefficient σ = 0.072 N/m. The gravitational acceleration is g = 9.8 m/s 2 . Because of the high Reynolds number 3.8 × 10 6 reported in the experiment, we discretise the computational domain on a relatively high-resolution grid that has 1610 × 600 cells in order to extract the data more accurately in the flow field. No-slip conditions are imposed on the enclosing walls. In this simulation, the downstream distance of the domain is sufficiently long to imitate a strong shear flow situation. The interface is unstable when there is a large velocity difference across it with small perturbations between two fluids. This phenomenon is known as Kelvin-Helmholtz instability (KHI). According to an analysis using linear theory (Chandrasekhar 1961), gravity and surface tension force play a role in stabilising the interface. Substituting typical parameters in the case of a collapsing water column into the formula for the linear growth rate of KHI: where k is the wave number of the interface perturbation, yields a negative value within the square root. Therefore, the Kelvin-Helmholtz instability shown in Figure 10(a) shouldn't occur. In contrast, the consistent method more realistically reproduces the profile of the water front, which compares favourably to the experiment and another numerical simulation using the SPH method (Rezavand, Zhang, and Hu 2020). It reveals that solving the non-conservative momentum equation will result in that if the low-density fluid has a large velocity, it will have a large impact on the movement of the high-density fluid. The time evolution of the water surface obtained by the consistent method is compared with the experimental snapshots in Figure 11. Except in the early stage of collapse, which is under the influence of the gate release, the surface profiles, water heights and front positions are in good agreement. Figure 12 shows the time variation of the impact pressure at four monitoring points on the downstream wall in comparison with the sensor data in the experiment. Non-dimensional time t * = t g/H and pressure p/(ρgH) are used, where ρ is the density of water and H is the initial height of the water column. The colour band in the figure represents the experimental data between the 2.5% and 97.5% percentile levels of 100 tests. In Figure  12(a), the rising point of the impact pressure measured in the numerical simulation is a little ahead of that of the experiment. The peak pressure is also higher. This can be explained by the fact that the water column is not subject to resistance from the side wall in the twodimensional simulation, and the water front hits the lower corner of the tank earlier with a slightly higher speed. Nevertheless, the simulation is in satisfactory agreement with the impact pressures measured by the four sensors.

Three-Dimensional Rising Bubble
To validate our conservative and consistent solver for realistic three-dimensional two-phase flow simulations, we simulate the rising of a single bubble in a liquid under various conditions and compare the results with the experiment conducted by Bhaga and Weber (1981). The final bubble shapes are categorised based on the non-dimensional Eötvös number, also known as the Bond number, and the Morton number: In our simulation, the density ratio and the viscosity ratio are set to 1000 and 100, respectively. The liquid density is ρ h = 1000 and the bubble density is ρ l = 1. The viscosities of the two phases are μ h = 10 −3 and μ l = 10 −5 . The bubble diameter is D = 0.01. Once the Eötvös and Morton numbers are known, the gravitational acceleration g and the surface tension coefficient σ can be determined accordingly. Reference numerical simulations (van Sint Annaland, Deen, and Kuipers 2005;Hua and Lou 2007) have shown that a computational domain with lateral dimensions of four Before comparing our results with the experiment, we examine grid independence. We compute the case with Eo = 116 and Mo = 41.1 on several different uniform grids, where the bubble diameter is resolved by 15 ∼ 35 grid cells. Figure 13 shows the contour of the bubble shape on the central slice of the domain at the non-dimensional time τ = t/ D/g = 10. Our results show that a grid finer than 120 × 120 × 300 had no significant effect on the predicted bubble shape. Therefore, we use the 120 × 120 × 300 grid with a diameter of 30 cells for the following simulations.
The test cases listed in Tables 3 and 4 cover almost all typical bubble shapes, including spherical, oblate Figure 13. Effect of mesh resolution on the simulated bubble shape. ellipsoidal, oblate ellipsoidal cap, spherical cap and skirted shapes. The predicted bubble shapes using our consistent solver agree very well with the experiment (Bhaga and Weber 1981). It is worth noting that cases A5, A6, B7 and B8 have relatively high Reynolds numbers, and the motion of the bubble does not reach

B8
Eo = 116 Mo = 8.60 × 10 −4 the steady state, making the bubble easy to break up. Therefore, we just select one of the snapshots in the simulation to match the experiment. These phenomena are also observed in the reference numerical computation (Hua and Lou 2007).

Liquid Jet in Cross Flow
To demonstrate the stability of our method as a practical two-phase flow solver, we simulate the problem of a three dimensional jet in an air cross-flow. This is a suitable validation of the importance of consistent transport and momentum conservation when applied to realistic two-phase flows with a high density ratio and high Reynolds number. An experimental study on this problem has been reported in Sallam, Aalburg, and Faeth (2004), and numerical studies have also been published in Li and Soteriou (2012), Li and Soteriou (2016).
The geometric and physical parameters follow the computation conditions described in Mirjalili and Mani (2021). The domain size measures 2.0cm × 1.5cm × 3.5cm and is discretised by the 256 × 192 × 448, 384 × 288 × 672 and 512 × 384 × 896 uniform grids with different resolutions. The liquid nozzle of diameter d = 0.8 mm is placed on the bottom wall and centred at (0.2cm, 0.75cm, 0cm). The air is blown in from the left boundary in the positive x direction with a constant velocity V g , and the liquid is injected from the nozzle along the positive z direction with a constant velocity V l . The physical properties are listed in Table 5. In case 1, the density ratio is 100, while in case 2, the density ratio corresponds to the actual value for water and air. However, the momentum ratio q, Weber number and Reynolds number are the same for both cases, indicating that both cases belong to the multimode breakup regime (Sallam, Aalburg, and Faeth 2004). No-slip conditions are imposed on the bottom wall z min except for the nozzle. Outflow conditions are applied to the right boundary x max . Other walls at y min , y max and z max have slip conditions. Wu et al. (1997) proposed the equation of the liquid column trajectory as where the drag coefficient fitted by the experiment of Sallam et al. for the multimode breakup regime is     Figures 14 and 15 show the liquid surfaces at t = 2.0 × 10 −3 s obtained by the inconsistent and consistent methods for case 1. By comparing these results with the experimental results in the literature (Sallam, Aalburg, and Faeth 2004), it is obvious that the consistent method gives more robust simulation results on the interface shape and the trajectory of the jet. The results of the inconsistent method show that the jet column breaks up within a short distance after leaving the nozzle, influenced by the air cross-flow. Moreover, the droplets formed after break up are quickly blown downstream with the air. Even after refining the mesh, the liquid column breakup characteristic obtained by the inconsistent solver does not improve significantly. Such a phenomenon in the numerical simulation is in line with our expectations, if the mass and momentum transports are inconsistent, the large momentum stored in the heavy fluid will leak into the light fluid, causing the velocity of the light fluid to increase dramatically. In turn, the light fluid with a high velocity will have a strong non-physical effect on the motion of the heavy fluid. Finally, we use our consistent method to compute case 2 on the finest uniform grid. The time evolution of the jet interface till 2.0 × 10 −3 s and the velocity field on the central section of the computational domain are depicted in Figure 16. The computation remains stable, although the two-phase flow with a high density ratio becomes very turbulent. The trajectory of the liquid column is also in good agreement with other numerical simulations (Li and Soteriou 2016;Mirjalili and Mani 2021) and an experimental study (Sallam, Aalburg, and Faeth 2004). Considering the complexity of the liquid jet phenomenon, further numerical simulation with much higher grid resolutions are needed to better understand the formation and evolution of the jet structure.

Conclusions
Building on our previous work on developing the evolving pressure projection method for solving weakly compressible Navier-Stokes equations, we propose a novel conservative and consistent momentum transport method for simulating incompressible twophase flows. To achieve consistent transport of mass and momentum, this method employs a collocated grid system on which the finite volume discretisation is straightforward to simultaneously solve the VOF and momentum equations in conservative form. A simple consistent formulation can be constructed by calculating the momentum flux with the help of the volume flux in the VOF equation. Our new implementation of the evolving pressure projection method also effectively couples the pressure and velocity on the collocated grid, thereby enhancing the numerical stability of the present solver.
The proposed method's ability to preserve momentum conservation at an arbitrary high density ratio is substantiated through various benchmark tests for two-phase flows, encompassing rising bubble, transport of a heavy droplet, oscillating drop, collapsing water column and jet in a cross flow. The conservative and consistent method ensures accurate and physically realistic results by conserving momentum and accounting for the interactions between the phases. This is particularly important for simulating violent two-phase flows with high density ratios, where momentum transfer between the phases can have a significant impact on the behaviour of the complex interfacial dynamics. In addition to its simplicity and efficiency in calculating pressure evolution equations, our momentum-conserving weakly compressible Navier-Stokes solver serves as a robust and high-fidelity tool for two-phase flow simulations.