An enhanced momentum conservation treatment for FDM simulation of two-phase ﬂows with large density ratio

The differences of the ﬂuid properties across a ﬂuid interface in two-phase ﬂow often bring diﬃculties into computational simulations, as the conservation of mass, momentum and energy requires careful consideration at the interfacial region. Velocity advection and unsynchronised variables lead to loss of conservation of momentum across the interface, which results in an unphysical interface deformation and spurious interfacial currents. In this study, we investigate the numerical errors and instabilities in the interfacial region, and propose a new algorithm with strong temporal coupling and momentum-based velocity reconstruction, to enhance its conservation properties. The capability of the proposed algorithm is demonstrated with two idealised cases including a one-dimensional convection case of a dense droplet and a standing wave case, and one laboratory dambreak case. Results are compared with theoretical results, experimental data or existing simulations, which demonstrate the advantages of the proposed algorithm on the conservation of mass, momentum and energy, and the mitigation of unphysical interfacial transport. Without modiﬁcation of any numerical methods or discretization schemes, the algorithm keeps its simplicity and can work with the existing methods, and it is straightforward to implement.


Introduction
Two-phase flow is a frequent phenomenon in the real world, involving the interaction of two fluids with different physical properties.It is of great importance in both industry and as a research area in its own right.The presence of twophase flow is observed over a variety of scales, e.g.molecular movement [1], lab-scale breaking waves [2], ocean and coastal waves [3], or even at the planetary scale [4].With the development of computational methods, numerical simulation has become a widely used tool to study two-phase flow [5], as the field of Computational Fluid Dynamics (CFD), has flourished.
The numerical simulation of two-phase flow requires the proper treatment of the behaviours of the separate fluids and their interface.For grid-based CFD methods like the Finite Difference Method (FDM) and Finite Volume Method (FVM), there are two main approaches: Two Fluid Method (TFM) [6] and One Fluid Method (OFM) [7].In TFM, the Baer-Nunziato model is proposed to describe the two-phase flow [8], with the known problem on conservation that the convective part in the model cannot be written in a conservative way [9].Many methods have been developed to deal with the non-conservative terms, like [10][11][12].The Baer-Nunziato model often requires 4 or 7 governing equations, indicating a larger demand for computational resources.Navier-Stokes equations (NSE) are also used in TFM, in which both fluids are described explicitly, thus two sets of NSE are needed [13].The velocity slip and properties transfer between fluids are allowed [14], however, the exchange of mass, momentum and energy across the interface often suffers from large uncertainties, and surface capturing methods is still needed [15].Additionally, TFM is more applicable for simple cases rather than realistic scenarios [16].
OFM assumes the same fluid occupies in the whole domain, but its physical properties vary temporally and spatially, i.e. the fluids in the whole domain are treated as a mixture.Only one set of NSE is needed to perform the simulation, which significantly saves simulation time.The fluid interface in OFM is described implicitly [17], and an accurate interface definition is essential for the accuracy and stability of OFM.As the interaction between the two fluids dominates the shape of the interface and affects the behaviours of the two fluids, it is essential to accurately capture the fluid interfaces.Interface tracking methods (e.g., Front-tracking method [18] and Marker-And-Cell method [19]) and interface capturing methods (e.g., the Volume of Fluid (VOF) method [20] and the Level -Set (LS) method [21], or the combination of the two, as with the CLSVOF method [22]) have been developed and employed widely.To avoid errors in interface definition, these surface capturing methods need to be coupled appropriately with the NSE, to make sure the variables and surface location are updated in a proper manner, especially in high-order temporal schemes, as the variables are updated more frequently inside one single timestep.
However, conservation of momentum and mass may still be improved across the fluids interface.One origin of errors is the advection term of the incompressible NSE.Large density ratios may lead to a large momentum ratio across the interface, and increase the risk of unphysical interfacial momentum exchange.This error may lead to an unphysical velocity field across the surface, resulting in a deformed surface profile and interfacial spurious currents [23,17], or an inaccurate interface definition.
Another source of errors is the limitation of the surface capturing methods, e.g., the LS method is known to have the problem on volume and mass loss in under-resolved regions [24] while the VOF method often generates inaccurate curvature and lower smoothness at the interfacial region, which may lead to an inaccurate velocity field [25].Numerous treatments or adjustments have been developed to deal with these foregoing problems.E.g., geometric operations [26][27][28] are implemented into VOF to improve the accuracy of the interface shape; extra smoothing [29], variable density projection [30], additional diffusion equations [31,32], artificial compression and 'viscosity' have been introduced into the LS method; special boundary conditions are applied on the fluid interface [33]; adaptive mesh refinement has also been used close to the interface to improve accuracy [34]; higher order surface capturing methods have been developed to reduce the numerical error [17,35,36]; interface tracking methods have also been coupled with interface capturing methods [37], as well as particle methods [33].
However, with these modifications, algorithms have become more and more complex and surface capturing methods tend to lose their own simplicity [38].The balance among accuracy, stability, efficiency and simplicity has become a pressing problem.Meanwhile, the standard forms of these methods or algorithms are still widely used, like the classical LS method [39] or algorithms based on the classical LS method [40,41].
To address these limitations, the present study proposes an alternative approach.Rather than modifying the existing methods or discretization schemes, we are seeking a more appropriate computational procedure that can improve the accuracy and stability with least effort while reducing the need for modifications of the numerical methods.In this paper we propose a momentum-based velocity reconstruction step, coupled with the surface capturing process in a synchronised manner, which enhances the consistency and improves the interfacial smoothness, stability and conservation.As no change is made to the underlying scheme or method, our modification can work with many different schemes and improved surface-capturing methods, while maintaining simplicity and ease of implementation into existing CFD codes.
The paper is organized as follows: the basic theories are introduced in Section 2. The origination of the error is discussed in Section 3. The algorithm modifications are proposed in Section 4. The algorithm is then tested on an idealised case and a laboratory dambreak case to demonstrate its performance in Section 5. Conclusions are presented in Section 6.

Numerical methods
The FDM-based open-access CFD code REEF3D [42][43][44] is used as the underlying tool of the present study.The governing equations and key numerical schemes are introduced in the following [45].

The incompressible isothermal NSE with constant density and viscosity
The governing equations are the conservation laws of mass and momentum for incompressible flows, written in the nonconservative form as follows: where u i is the velocity, ρ the density, p the pressure, ν the kinematic viscosity, g i the gravitational acceleration.Turbulence is ignored in the present study.
A standard two-step projection method [46] is used to solve the incompressible NSE.In the prediction step, the pressure gradient term is ignored to solve for an intermediate velocity field u int as: In the correction step, the fluid velocities are updated by considering the pressure gradient as: Applying the incompressibility condition (i.e., Eq. ( 1)) to Eq. ( 4) to gives the Pressure Poisson Equation (PPE) as: Solving the PPE yields the fluid pressure at the current time step, which is then used to update the fluid velocities through Eq. (4).

Level Set method
The classical Level-Set (LS) method [21] is used to capture the water-air interface.An LS function ϕ(x, t) is defined that represents a signed distance, which satisfies |∇ϕ| = 1.According to the ϕ values of each computational node, the fluid condition of nodes can be determined as follows: The update of the LS function ϕ is achieved by solving the following advection equation: To ensure that ϕ remains a signed distance function, after the advection, the ϕ value at each node is re-initialized as follows [47]: In simulation, a Heaviside function is defined as follows according to the ϕ value: The density and viscosity of nodes are then calculated by considering the Heaviside function value, as shown in Eq. (10): With Eqs. ( 9) and ( 10), the water-air interface is represented by an interfacial region with a distance of x ( x is local grid spacing) towards both sides of the interface, which ensures the smooth variation of fluid properties across the interface.

Spatial and temporal discretization schemes
A staggered grid system is applied in the simulation.Scalar variables (e.g., density, viscosity and LS value ϕ) are defined at cell centres, while vector variables (e.g., velocity, momentum) are defined at the centre of cell edges, as shown in Fig. 1.
We have used three schemes to perform the spatial discretization for the advection term in Eqs. ( 3) and (7), and found the 5 th order weighted essentially non-oscillatory (WENO5) produces the best result (based on the cases in the present study).The central difference scheme is applied for the spatial discretization of the viscous force term in Eq. ( 3), the pressure gradient in Eq. ( 4), and the pressure Laplacian and the velocity divergence in Eq. (5).
The time propagation of the advection term in Eq. ( 3) is achieved using the 3 rd order Runge-Kutta Total Variation Diminishing (RK3-TVD) scheme, which has 3 rd order accuracy and small storage requirements.For the viscous term in Eq. (3), we used the implicit backward time centred space scheme as it allows larger time steps and has good stability properties.

Investigations on numerical errors near fluid interface
The numerical algorithms presented in Section 2 may lead to spurious interfacial kinematics and dynamics, as presented in Section 5 and previous studies [50,51].These numerical errors are attributed to the unphysical transport of momentum across the interface, which are caused by two primary factors: non-conservation of momentum in the spatial discretization of the advection term [17] and inaccuracy in the definition of fluid interface [52].Note that a higher interfacial density ratio would amplify this error.How the two factors erode the numerical accuracy are investigated below.

Non-conservation of momentum in advection term
One source of the non-conservation of momentum is the form of discretisation.Consider the conservative advection term and its nonconservative form u ∂x + ρu ∂u ∂x , and their discretised forms as follows: In a steady state, a unique horizontal flux at each point is ensured by the conservative form rather than nonconservative form.Adding up all flux terms along a given direction, (the x-direction is presented here), the conservative form provides zero net flux, while extra source terms arise with the nonconservative form.Thus, a conservative spatial discretisation is preferred for the conservation of momentum [57].
Although a conservative form is adopted, non-conservation of momentum around the fluid interface may still develop from the advection of velocity which does not consider the density change.Taking the one-dimensional (1D) advection case in Fig. 2 as an example, consider two inviscid fluids with a density ratio of ρ 1 /ρ 2 = 10 6 located either side of the interface n , (and let = 1, see Eq. ( 9)).In the n th time step, fluid 1 has a constant velocity u 0 while fluid 2 is at rest, and the velocity changes smoothly from u 0 to zero.In the subsequent (i.e.n + 1 th ) time step, the interface moves to n+1 .Following the strategy in Section 2, the velocity at point (i+0.5) is updated by velocity advection, which can be performed by different spatial schemes (e.g.WENO5, upwind or central difference).Regarding the density ratio, the case can be seen as a pure advection case and theoretically the velocity at n + 1 th timestep at point (i+0.5)should stay u 0 , with an updated momentum of u 0 ρ n+1 i+0. 5 .Table 1 shows the relative error of the numerical results of the updated momentum at point (i+0.5)obtained by taking the production of the advected velocity and the updated density (obtained from n+1 ) (Let the CFL number be 0.1).It is noted although the numerical results of WENO5 show the smallest discrepancies (due to its higher-order accuracy), errors of velocity advection close to the fluids interface (i.e. point (i+0.5))cause non-negligible momentum deviation and extra pressure.Thus, a velocity advection scheme that preserves the momentum is crucially necessary.

Temporal mismatch of fluid variables in time propagation
Time propagation schemes often consist of multi substages, which brings the consideration of temporal synchronization in the coupling of surface-capturing and NSE-solving steps.The present work adopts the RK3-TVD scheme [53] for the advection terms, which consists of three sub-stages as below: where f n , f (1) , f (2) and f n+1 are the fluid quantities at time t, t + t, t + 0.5 t, and t + t, respectively, and L( f n ) is the time derivative of fluid quantity at the corresponding time instant.Fig. 3 illustrates the updates of the LS and velocity values in a time step of the RK3-TVD scheme, which is widely used in CFD simulations [32,39] and some variants are also presented by [54,55].Note that the weak coupling between LS and NSE at the 2 nd and 3 rd substages may lead to an inaccurate interface definition, which can result in non-conservation of momentum.The first concern is that the update of fluid quantities, (obtaining ρ n+1 and ν n+1 according to ϕ n+1 ), is completed before solving the NSE, like [42], and no coupling is found between ϕ (1) and u stage1 (also ϕ (2) and u stage2 ).The second concern is the selection of u stage 1 and u stage 2 , as using u n for u stage 1 and u stage 2 may generate a temporal mismatch.Temporal extrapolation based on u n and u n−1 [17,57] fulfils the requirement of time matching, but the coupling between LS and NSE at these two substages is still not considered, as the temporal extrapolation solely relies on the velocity field.Alternatively, u stage 1 and u stage 2 can be obtained by solving the corresponding NSE at these two substages (t + t, t + 0.5 t) [45], which requires the corresponding density and viscosity field (thus ϕ (1) and ϕ (2) ).Thus, a proper coupling between LS and NSE is needed.Additionally, it is noted the above analysis is based on using the same time propagation scheme for both LS and NSE.Applying different time propagation schemes to them makes the synchronization more difficult [56].

Influences of interfacial thickness
To capture more flow details and achieve a higher accuracy, a thinner fluid interface is often needed, which requires a smaller grid spacing x (grid refining) or a smaller interfacial thickness parameter (interface sharpening).However, these adjustments may increase the risk of non-conservation of momentum [23].In the illustrative example of Fig. 2, points (i − 1) and (i + 1) are both in the interfacial region and point (i) is on the interface.According to Eqs. ( 9) and (10), the density at the three points can be calculated as: (13) Note that a lower x or increases the local density gradient across the fluid interface.As explained in [17], large density gradient close to the interface may aggravate the spurious momentum transfer and lead to the non-conservation of momentum.
Additionally, the error in the advected velocity (e.g., point (i+0.5) in Fig. 2) leads to an unphysical velocity gradient in the denser fluid, which may also increase the instability at a lower x or .According to Fig. 2, reducing x could increase this unphysical velocity gradient and results in a larger error on the local pressure field by solving the PPE (Eq.( 5)).This error may be further amplified due to the increased local pressure gradient at a lower x through Eq. ( 4).Meanwhile, reducing increases the updated density at point (i+0.5)(i.e., ρ n+1 i+0.5 ), which leads to a higher error on the updated momentum and adversely affects the momentum conservation.In summary, as will be proved in Section 5, the interfacial instability develops faster with grid refining and interface sharpening with the strategy introduced in Section 2.
Compared with the original scheme in Fig. 3(a) which only requires one re-initialization in each timestep, the proposed manner performs three times of re-initialization, which may amplify the well-known volume-loss problem of the classical LS method [24], and the re-initialized LS value (ϕ (1) , ϕ (2) ) may also break the consistency of the RK3-TVD scheme of the 2 nd and 3 rd substages.Thus, the un-reinitialized LS values (ϕ (1) * , ϕ (2) * ) are accepted for the advection of the next substage, instead of the re-initialized LS values (ϕ (1) , ϕ (2) ) which are solely used for updating fluid properties.Only the re-initialized LS value in the last substage (ϕ n+1 ) is accepted for the advection of the next timestep, to avoid the potential error from multi re-initializations in one timestep.

Momentum-based velocity reconstruction
Regarding the issue of velocity advection (Section 3.1), we propose a velocity-reconstruction technique.Based on the conservative form of NSE (Eqs.( 14) and ( 15)), we perform the conservative momentum advection, as shown in Eq. ( 16): in which M is the momentum calculated from ρ • u.The advections of velocity and density are also tracked simultaneously for further use, as shown in Eqs. ( 17) and (18).Eq. ( 19) shows the explicit form of advections of M, u and ρ.
The advected velocity field u * ,r is then restored based on the quotient of M * and ρ * , as shown in Eq. ( 20): For the spatial matching of the density and momentum of Eq. ( 20) in a staggered grid system, the density advection (Eq.(19c)) is performed in a vectorised manner based on the interpolated face-centred density.
The velocity reconstruction process shares the same framework as [52,58,59] in which the velocity is restored in the same way, but differences appear on the convective density and momentum.The density flux often relies on the surfacecapturing method, e.g. the density flux in [52] comes from the VOF flux, and [52,[58][59][60] trace the motion of the interface from the LS values and reconstruct the density flux by geometric arguments.Compared with the straightforward Eq. (19c), the geometric argument may have extra limitations, like [60] requires a sharp interface (ϕ = 0), as the convective density is hard to define if the density changes continuously across the interface.Additionally, more effort may be needed to reconstruct a complex interface shape, leading to an increasingly complex geometric argument, which makes it challenging to perform in 3D conditions [61].It is noted the intermediate density (ρ * i ) from Eq. (19c) should be matched with the intermediate momentum (M * i ), rather than the density defined by LS method at the end of the timestep (ρ n+1 i ).The momentum flux is taken as the production of advected velocity and density flux in [52,60,61].However, we still perform the momentum advection based on the advection equation (Eq.(19b)), to ensure the consideration of velocity gradient (recall that ∂ M ∂x = u ∂ρ ∂x + ρ ∂u ∂x ).As the same advection equation, velocity field and numerical operation are applied on the update of the intermediate density and momentum, better spatial and temporal matching between them is ensured.
A simple 1D advection test is designed to demonstrate the conservation of momentum during interface moving.As shown in Fig. 5, an interface divides a 1D domain with 400 gridcells, (grid spacing x=0.005 m), and the denser fluid 1 is on the left and the less dense fluid 2 on the right, (the density ratio being 10 6 ).Fluid 1 is initialized with the velocity of ±1 m/s and the CFL number is 0.1.Open boundary conditions are applied on both sides, and gravity and viscosity are ignored.No interface deformation is supposed to appear in this case.
This case is simulated by both the proposed velocity reconstruction technique and the original velocity advection method.Compared with the analytical solution, the momentum discrepancies of both methods after 1 s (1000 loops) are shown in Fig. 6, in which it is observed the momentum conservation is improved by the velocity reconstruction technique, as the  momentum change is much more stable and is no more than 0.04% of the initial momentum of one cell of fluid 1 with WENO5 and upwind schemes.
After the velocity reconstruction, the viscous and gravitational terms in Eq. ( 3) are integrated into the restored velocity field, yielding the intermediate velocity field u int .Eq. ( 4) and ( 5) are then solved for pressure update and velocity correction.
In contrast to the method proposed by [17] in which the intermediate density field ρ * i is used in Eq. ( 4) and ( 5) or similarly [61] who uses ρ n+0.5 i , ρ * i is immediately discarded after the velocity restoration in our strategy.As the LS value at the next iteration is already known, ρ n+1 is used in Eq. ( 4) and ( 5), which is defined by the corresponding ϕ n+1 .
Potential errors may rise in the velocity restoration, and we applied 2 treatments: 1) applying the same spatial schemes for variable advections (momentum M, velocity u, density ρ and LS value ϕ); 2) applying a 'limiter' treatment around the boundary between the interfacial region and less dense fluid.The details are described in the following sections.

Control the interaction between truncation errors
The truncation errors of M and ρ may interact with each other in Eq. ( 20) and reduce the accuracy of the restored velocity.To tackle this undesirable effect, we suggest that the advection of u, M and ρ follow the same spatial scheme.Consider two simple 1D convection cases as shown in Fig. 7, let the truncation error be proportional to the accurate advection term, the updated values of u, M and ρ, denoted by asterisks, are shown in Eq. (21).The coefficients a 1 , a 2 and a 3 denote the proportionate truncation errors of u, M and ρ.
In the case of an accelerating fluid of constant density, Fig. 7(a), density gradient is not accepted after the advection process.Hence, the restored velocity u * ,r i (Eq.( 20)) should be the same as the directly-advected velocity u * i (Eq.(19a)), as shown in Eq. ( 22), requiring the equivalence of a 1 and a 2 .The constant density leads to the same term of ∂u n x u n x ∂x on both sides of Eq. ( 22), thus applying the same advection scheme on u and M meets the requirement.
Now consider the case of uniform flow of a fluid with density gradients, Fig. 7(b).All the properties are transported by u 0 , and the restored velocity u * ,r i should stay at u 0 , as shown in Eq. ( 23), requiring the equivalence of a 2 and a 3 .The constant velocity leads to the same term of ∂ρ n x ∂x on both sides of Eq. ( 23), thus applying the same advection scheme on M and ρ meets the requirement.
Starting from these two simple cases, there are two basic requirements: a 1 = a 2 (same advection scheme on u and M) and a 2 = a 3 (same advection scheme on ρ and M).Ideally, we would seek an advection scheme for which a 1 = a 2 = a 3 .
In general, this is not possible as this would impose an unphysical restraint on the gradients of the advected quantities.
Nevertheless, applying the same scheme on the advection of u, M and ρ may still help to reduce the interaction between the truncation errors and its necessity is discussed in Section 5.1.

The 'limiter' treatment around the interfacial region
While density changes smoothly in the interfacial region, unbounded density may develop due to the numerical overshoot or undershoot in the density advection, which strongly influences the velocity restoration.This topic has been discussed by [52,61].
The following advection case in Fig. 8 gives a brief illustration.Fluid 1 with high density (ρ 1 /ρ 2 = 10 6 ) is initialized by a velocity u 0 and the light fluid 2 is at rest.Viscosity and gravity are ignored.Taking = 2.1 and CFL=0.1, the boundary of the interfacial region is set to be k x away from point (i-0.5).Based on Eq. ( 19) and Eq. ( 20), the restored velocity at point (i-0.5) is shown in Eq. (24).Compared with Eq. ( 25), the ratio between the original and advected density (ρ n /ρ * ) is the key contribution to the velocity reconstruction.Should ρ * become small, zero or even negative, this is likely to cause instability without further intervention.
ρ n /ρ * is lower than unity at a positive u 0 , reducing the influence of the velocity gradient in the interfacial region (see Eq. ( 24)), and reduces the unphysical effect of the velocity of fluid 2 on denser fluid 1.In contrast, a negative u 0 may strengthen the effect of fluid 1 on the less dense fluid 2. Fig. 9 shows the results of ρ n /ρ * at point (i-0.5)calculated by 3 commonly used schemes (WENO5, central difference and upwind).As shown in Fig. 9(b), with large density gradient and low local density, the density advection of Eq. (19c) may lead to an unbounded density at a negative u 0 , as ρ * might be lower than the less dense fluid (or even negative, as observed in [57]), resulting in an unphysical value of ρ n /ρ * , which is especially noticeable with the central difference scheme.This numerical overshoot or undershoot is also observed in the advection of the volume fraction in VOF method which relies on the same advection equation as the volume fraction may be lower than 0 or higher than unity [62,63].As mentioned in [61], the density flux must meet the critical requirement of boundness.A limiter treatment is then applied on the nodes adjacent to the boundary with density close to the lighter fluid.Thus, the velocity restoration follows Eq. ( 26) instead of Eq. ( 20), as shown below: in which ρ * , f i is a filtered density value following Eq.( 27): For numerical stability, the 'threshold density' in Eq. ( 26) is set to be 0.05ρ 1 + ρ 2 , and a slightly larger threshold density is suggested when applying the central difference scheme due to the stronger numerical overshoot and undershoot.

The proposed algorithm
Combining the strong time-coupling (Section 4.1), velocity reconstruction (Section 4.2) and the limiter treatment (Section 4.4), the flow chart of the proposed algorithm which enhances the momentum conservation is shown in Fig. 10.As the proposed algorithm focusses on momentum conservation, we have termed it the 'Enhanced Momentum Conservation' (EMC) treatment.For simplicity, the FDM code with the implementation of EMC treatment is termed 'FDM-EMC'.
It is noted the proposed algorithm does not include any modification of existing schemes or methods.As explained in the Introduction, the aim of this study is to establish a proper computational procedure that improves the accuracy and stability with the least requirement for modification of numerical methods.The strong time coupling (Section 4.1) is a variation of the well-established RK3-TVD scheme, and the velocity reconstruction (Section 4.2) utilizes the existing spatial and temporal schemes to perform the advection of momentum, density and velocity instead of only velocity.The effectiveness of these two techniques is tested separately in Section 5.1.
Although not aiming to modify the surface-capturing method, the mass-loss problem of the classical LS method is still partly countered by the EMC treatment, which is shown in Section 5.As all the schemes and methods used (e.g., RK3-TVD, WENO5, classical LS and its re-initialization, projection method) remain the same, the EMC treatment retains the simplicity of these methods, and it can be implemented in a relatively straightforward manner.The EMC approach should also work with many existing variants, like the conservative LS method [48] or VOF method.

Numerical results
Two idealised numerical cases and one laboratory dambreak case are presented to investigate the performance of EMC treatment.In Section 5.1, a 2D pure advection case with large density ratio is chosen to detect the interfacial instability and study the conservation of mass, momentum and energy.In Section 5.2, an inviscid standing wave case is used to investigate the conservation in the absence of intensive interface deformation.In Section 5.3, an experimental dambreak case with a dry bed is chosen to demonstrate the advantage of EMC on realistic cases.

Table 3
Code settings for the sensitivity studies 1 to 3.

Code settings Modification applied
Advection scheme for ρ Advection scheme for M and u Density ratio R ρ *VR -Velocity Reconstruction; TC -Temporal Coupling; S0 is the original FDM code (REEF3D) without modification.

Table 4
Code settings for sensitivity study 4.

Code settings Modification applied
Convection scheme for ρ Convection scheme for M and u Density ratio R ρ

Advection of a dense inviscid droplet
The advection of a dense droplet has been investigated by various researchers [17,57,64], as a large density ratio makes it easier to detect unphysical momentum transport and spurious current across the interface.In this case, a droplet with diameter D=0.3 m is initialized with a horizontal speed u 0 (D/u 0 =0.3 s) while the surrounding air is at rest.The density ratio (R ρ ) between the droplet and air is 10 6 (10 3 if the simulation crashes too early).Viscosity, surface tension and gravity are ignored and open boundary conditions are applied at each side of the domain.The CFL number is set to be 0.1.Due to the extreme density ratio, it is supposed that the influence of the air on the droplet is negligible, and the droplet travels with a constant speed (u 0 ) without shape change.
As shown in Tables 2 to 4, 5 uniform grid systems (G1 to G5) and different code settings are used for these sensitivity studies: 1) grid convergence; 2) the individual contributions of the strong time-coupling (Section 4.1) and the velocity reconstruction (Section 4.2); 3) the necessity of applying the same spatial scheme on density and momentum (Section 4.3); 4) the model performance under different advection schemes.
Fig. 11 shows the results of the grid convergence study of code settings S0 and S1 under grids G1 to G5. Due to the non-conservation of momentum in the velocity advection (Section 3.1) and the temporal mismatch (Section 3.2), no obvious grid convergence is observed in the results of S0, which also loses the smooth disk shape of the droplet.The stability and interface smoothness are downgraded under a thinner interface thickness (G3) or a finer grid (G4, G5), which agrees with the analysis in Section 3.3.However, lower x or are allowed in the simulation with EMC treatment, as S1 achieves grid convergence and retains the smooth disk shape of the droplet, while its location also agrees with the reference one.It is noted the interface deformation is visible at a resolution close to G4 in [17], while S1 allows a resolution of 1 time finer (G5).
Fig. 12 shows the evolution of the volume, horizontal momentum and kinematic energy of the droplet, presenting a quantitative comparison between the results of S0 and S1.Suffering from the instability explained in Section 3 and the wellknown mass loss problem of LS method [49], the time series of mass, volume and momentum show obvious oscillations and a decreasing trend in the results for S0 in Fig. 12.Although there is no modification to the classical LS method, the implementation of EMC shows a significant improvement in its conservation properties, especially under the finer grids (G4, G5), as the loss of all 3 variables is significantly reduced, (compare the solid and dashed line with the same colour in Fig. 12(a) to (c)).As an approximate means of filtering out the influence of mass loss on energy loss and study the energy change from the interfacial transport, we compare the ratio between energy and volume in Fig. 12(d), from which   it is concluded the energy loss is mainly from the volume loss, as the ratio is always close to 100% with both S0 and S1.However, the superiority of S1 is still confirmed by its more stable and converged time series and the negligible energy change per unit volume.
To identify the individual contributions of the strong time-coupling (Section 4.1) and velocity reconstruction (Section 4.2), the results of S0 (original FDM code) S1 (EMC fully implemented), S2 (only strong time-coupling implemented) and S3 (only velocity reconstruction implemented) under grid G4 are compared in Fig. 13(a) and Fig. 14.As shown in Fig. 13(a), the smooth disk shape only retains in the result of S1 while slight interface twisting is already observed at an early stage (T=1) in S0 (most obvious), S2 and S3.At T=2.33, the disk shape is distorted in S3 and S0 (slightly more intense) while S2 roughly keeps the disk shape despite the deformation at the rear side of the droplet.Runs S0 and S3 crash at T=2.50 and T=3.03, respectively.At T=3, significant interface shape deformation is eventually observed in S2.Fig. 14 compares the time series of volume, momentum, energy and energy per unit volume of the four code settings, in which S0 is the most unstable.All S1, S2 and S3 show a decreasing trend on volume, momentum and energy, while S1 shows the best performance (reduction is less than 0.25% at T=3.3) and S2 outperforms S3.It is noted the reason of the rising trend in Fig. 14(d) is that the energy loss is slower than the volume loss.Compared with S0, both S2 and S3 show improvement on interface smoothness, numerical stability and property conservation, and the performance of S3 is significantly better, as the velocity reconstruction directly focuses on the momentum conservation and the aim of strong time coupling is the synchronization of variables.However, satisfactory disk shape is only observed in S1 in which the time coupling and velocity reconstruction are both implemented.Thus, it is concluded that the main contribution of EMC treatment is the velocity reconstruction, but the strong time-coupling is also essential.
To prove the necessity to use the same advection scheme for both density and momentum, the results of S4 and S5 are compared in Fig. 13(b).Although both WENO5 and the WENO5 for Hamilton-Jacobi equation (WENO5-HJ, which is more suitable for Hamilton-Jacobi equations [65] and frequently used for LS advection, but less conservative than normal WENO5) can achieve 5 th order accuracy, S4 and S5 crash immediately under a density ratio of 10 6 .As shown in Fig. 13(b), in the  simulation with density ratio of 10 3 , both of them failed to keep the disk shape of the droplet even at T=1, indicating the virtue of using the same spatial scheme for the advection of ρ, M and u.The same conclusion is obtained from different codes, for example [64], in which the distorted interface and larger velocity errors are detected when dissimilar convection schemes for mass and momentum are applied.
To test the performance of the EMC treatment under different spatial schemes, simulations under grid G2 are performed with different convection schemes shown in Table 4, and the interface shapes are compared in Fig. 15.With lower order schemes, S0 quickly loses the disk shape, and instability can still be observed with a high order scheme like WENO5-HJ.In comparison, with all schemes tested, the EMC treatment significantly improves the stability and smoothness of the interface although the performance is much better with higher order schemes.

Standing wave case
A standing wave case is designed to test the effectiveness of EMC in a case with more gentle surface variation.As shown in Fig. 16(a), the computational area is 5.5 m long with still water level 0.5 m.A sinusoidal wave with wave height 0.01 m and 2 m wavelength (period 1.182 s) is generated by a relaxation area of 2 m length.The wave propagates towards right and hits a perfectly reflecting boundary, the superposition of the incident and reflected waves forms a standing wave in an area of 3.5 m length.The reflected wave propagates out of the active area and is absorbed by the relaxation area [66,67].The densities of water and air are set to 1000 kg/m 3 and 1 kg/m 3 respectively, and g is set to 9.81 m/s 2 and the fluids are considered inviscid.A wall boundary condition is applied to the top and bottom boundaries and a reflection boundary condition is applied at the right wall.WENO5 and RK3-TVD schemes are used for the advection of variables and the CFL number is 0.1.The interfacial thickness parameter of LS ( ) is set to 2.1.Three uniform grid systems are applied, namely a coarse grid ( x=0.02 m), a medium grid ( x=0.01 m) and a fine grid ( x=0.005 m).The simulation is performed for 120 s with both the original FDM-based code and FDM-EMC.The simulated wave profiles of the standing wave between 57.3 s and 57.9 s are shown in Fig. 16(b).
Fig. 17 and 18 compares the vertical momentum of water and air in the last 20 seconds.FDM-EMC shows a good grid convergence, and the momentum of water oscillates around zero with a stable amplitude, either water or air (Fig. 17(a) and 18(a)).However, FDM leads to a downgraded grid convergence, and a less stable momentum oscillation, especially with coarser grids, indicating the unphysical momentum change across the water-air interface (Fig. 17 The time series of the total energy of water in the last 20 seconds are compared in Fig. 19 and the average volume of water of the last 10 wave periods is shown in Table 5.Although there is no intensive water-air interaction and visible interface deformation under this low wave height (0.02 m in 0.5 m water depth), FDM-EMC still provides a more stable     oscillation and a better convergence for the total energy.Both codes produce a minimum volume change (no more than 0.33%) in this case, but FDM-EMC still slightly outperforms the original FDM because of the reduced volume loss.

Dambreak case with a dry bed
The third benchmark case is a lab-scale dambreak case with a dry bed with settings and parameters shown in Fig. 20 and more information can be found in [68].Three uniform grid systems are applied to the simulation, namely coarse  ( x = 0.01m), medium ( x = 0.0075m) and fine grid ( x = 0.005m), and the interfacial region parameter is set to be 2.1 [45].A wall boundary condition is applied on the left, right and bottom sides while an open boundary condition is applied to the top boundary.
The grid convergence studies of the FDM and FDM-EMC (S0 and S1 in Table 3) are compared in Fig. 21.It is seen FDM-EMC better keeps the smoothness of the water-air interface and produces a sharper overturning tongue than FDM.A better grid convergence is also found in the results of FDM-EMC as the water-air interfaces of different grid resolutions overlap together while larger differences are found in the results of FDM.
Fig. 22 compares the surface profile simulated by FDM-EMC and several existing simulations [70,69,71] with the experimental snapshots [68].As shown in the top row of Fig. 22, the simulation results all agree well with the experiment before the collision between wavefront and left wall.The intense water-air interaction leads to the twisted water-air interface since t=573.3ms in VOF simulation [69] while other codes still agree well with the experiment, as seen in the middle row of Fig. 22.For the fully developed breaker (t=1023.3ms), FDM-EMC generates a much sharper tongue and lower splash than VOF [71] and SPH [70].At 1166.6 ms, FDM-EMC and compressive interface tracking method [71] generates a highly similar entrapped air pocket, while the secondary jet from FDM-EMC is closer to the experiment.
The pressure distributions, water-air interface shapes and velocity vector maps reproduced by FDM-EMC and FDM are compared in Fig. 23.A twisted air-water interface is already observed in FDM from t=0.30 s as negative pressure areas and vortices develop at the front border which propagate in a retrograde manner along the interface.Numerical instability is visible on the interface shape at t=1.0 s which is characterised by unphysical negative pressure and water 'branching' around the breaking tongue.The instability at t=1.16 s causes the simulation to terminate, as the local vortex and negative pressure (as low as -2000 Pa) dominates the flow field around the secondary jet.However, FDM-EMC effectively suppresses the negative pressure and vortex around the water-air interface, as the water movement is less influenced by the air due to its improved momentum conservation.
The collision between the overturning breaker and underlying water body proves the necessity of the limiter, as shown in Fig. 24.At t=1.06 s, high gradients of density and velocity develop between the two separate waterbodies and the  surrounding air.According to Eq. ( 24), the unphysical density ratio leads to an incorrect velocity field and pressure field, culminating in the termination of the simulation without the limiter before the water-air interfaces merge.

CPU performance
As described in Section 4, the EMC treatment uses the strong time-coupling and momentum-based velocity reconstruction to enhance the conservation.Compared with the original code, additional advection of momentum and density are performed, leading to the reduction of numerical efficiency.The CPU time per loop of FDM and FDM-EMC are compared in Table 6.It is seen the CPU time rises about 50% at a lower grid resolution and 100% at a higher resolution.Future work would be conducted to improve the computational efficiency.

Conclusion
Simulation of two-phase flow suffers from spurious interfacial currents and instability.The associated lack of conservation of momentum arises chiefly from inaccuracies in velocity advection and synchronization of the variables.Grid refinement and interface sharpening may also amplify the interfacial instability.To deal with these problems, the EMC treatment is proposed, which consists of the strong time-coupling between LS and NSE, the momentum-based velocity reconstruction and the limiter treatment around the interfacial region.The EMC treatment has been implemented into an existing FDM code and significantly improves its performance.The simulation of the convection of a dense inviscid droplet demonstrates the improvement on the accuracy, stability and the conservation of mass, momentum and energy, and the capability of dealing with two phase flows with an extreme density ratio.Through the comparison study, it is concluded that the momentumbased velocity reconstruction contributes most to the EMC treatment, but strong time-coupling is also essential.A standing wave case confirms that the EMC treatment can provide improvements even under milder conditions without large interface deformation.In the dambreak case, the shape of the breaking wave is well reproduced, and the interfacial smoothness is also improved considerably.The distribution of pressure and velocity around the water-air interface has been analysed in detail, confirming the advantages of the EMC approach.
The EMC treatment requires neither modifications to the underlying numerical scheme or method, nor special treatments or additional boundary conditions applied to the interface.It therefore provides an additional way to improve the conservation of mass, momentum and energy as well as the smoothness of the fluid interface.Furthermore, the EMC treatment can work with existing algorithms and surface tracking methods, and it is straightforward to implement because of its simplicity.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 2 .
Fig. 2. A 1D case to illustrate the errors of momentum and density.

Fig. 3 .
Fig. 3. Updates of the LS (a) and velocity (b) values in a time step of the RK3-TVD scheme.

Fig. 7 .
Fig. 7. Two idealized cases for the demonstration of applying the same spatial scheme.

Fig. 8 .
Fig. 8.The initial settings of the advection case.

Fig. 10 .
Fig.10.The computational procedure of the proposed algorithm in one substage of RK3-TVD scheme.

Fig. 11 .
Fig. 11.Grid convergence study for S0 and S1 (time t is normalized by T=t*u0/D).Solutions are symmetric about the 'equator'.Results from each grid are shown for a hemisphere.See text for further explanation.

Fig. 12 .
Fig. 12.Comparison of the results of G1 to G5 on the evolution of (a) the total volume; (b) the total energy; (c) the total momentum of the droplet towards the x direction; (d) the energy per unit volume.(Dashed line -S0; Solid line -S1).

Fig. 13 .
Fig. 13.(a) comparison of the contributions of the modifications (b) results of EMC treatment when applying different spatial scheme for ρ, M and u.

Fig. 14 .
Fig. 14.Comparison of the results of S0 to S3 on the evolution of (a) the total volume; (b) the total energy; (c) the total momentum of the droplet towards the x direction; (d) the energy per unit volume.

Fig. 16 .
Fig. 16.(a) the settings of the standing wave case (not in scale); (b) the profiles of the simulated incident wave and standing wave.

Fig. 17 .
Fig. 17.The vertical momentum of the water in the active area (a) FDM -EMC (b) original FDM code.(blue -coarse grid, red -medium grid, black -fine grid).

Fig. 18 .
Fig. 18.The vertical momentum of the air in the active area (a) FDM -EMC (b) original FDM code.(blue -coarse grid, red -medium grid, black -fine grid).

Fig. 19 .
Fig. 19.The total energy of the water (a) FDM with EMC (b) original FDM code.(blue -coarse grid, red -medium grid, black -fine grid; energy is normalized by the initial potential energy of still water).

Fig. 20 .
Fig. 20.Settings of the dambreak case with a dry bed.

Fig. 23 .
Fig. 23.Comparison of the pressure distribution, interface shape and velocity vector around the leading front (white line is the water-air interface).

XFig. 24 .
Fig. 24.The interface shapes and pressure contours of FDM-EMC with and without limiter at t=1.06s.

Table 1
Relative error of updated momentum at point (i+0.5)by velocity advection.

Table 2
Grid settings in case 1.

Table 5
The comparison of volume loss (unit: %).

Table 6
The comparison of CPU time of the code with and without EMC treatment (All simulations are based on 8 threads).