A fast and improved mass-conservative approach for a sharp interface-level set method on co-located grid

1 Introduction

Developing an accurate and efficient numerical methodology for simulating multiphase flows with high density ratios remains challenging, particularly when aiming to preserve a sharp interface, utilize a co-located grid, and minimize mass loss. Among the various interface-capturing techniques—such as the Volume of Fluid (VOF) (Hirt and Nichols, 1981), the Phase Field Method (PFM) (Wang and Sekerka, 1996), and interface tracking approaches like Front Tracking—the Level Set Method (LSM) remains a popular choice due to its robust handling and accurate computation of the interface geometrical properties. Nevertheless, a major limitation of LSM is its inherent mass loss, which can compromise the fidelity of long-term simulations. In this study, we propose a novel formulation of the sharp interface level set method (SI-LSM) that not only addresses mass conservation but also achieves improved computational efficiency. Furthermore, on a co-located grid, the proposed formulation avoids pressure-velocity decoupling (Rhie and Chow, 1983) and pressure-interfacial force imbalance Chirammel et al. (2023b) for multiphase simulations with high density ratios.

The level set method captures the interface implicitly as the zero contour of a signed distance function, which is evolved through a standard advection equation and a reinitialization equation to preserve the signed distance property (Osher and Sethian, 1988; Sussman et al., 1994). A consequence of these level set advection and reinitialization equations is the shifting of the zero contour of the level set function, resulting in mass loss (Solomenko et al., 2017). Various strategies to resolve this mass loss can be categorized into three approaches: (i) hybrid methods, (ii) modifications to the conventional level set equations, and (iii) compressive-diffusive term-based formulations. In the hybrid approach, the level set method is combined with other mass-conserving methods, such as VOF (Rudman, 1998) in the coupled level set/volume-of-fluid (CLSVOF) method (Sussman and Puckett, 2000) or with Lagrangian particles in the hybrid particle level set method (Enright et al., 2002). Although these hybrid approaches are effective, they introduce substantial additional computational steps and computation cost—the particle level set method requires numerous particles per cell for accuracy, and the CLSVOF method involves an additional algorithm compared to the simpler level set method. In the second category, which modifies the conventional level set equation, strategies include adding a source or sink term to the level set advection equation to offset mass imbalance (Yuan et al., 2018), solving an additional volume-reinitialization equation (Salih and Moulic, 2013), or applying a global mass correction strategy (Ge et al., 2018) by solving an additional equation for a third distance function (Zhang et al., 2010). In the final category, Olsson and Kreiss (2005) solved a compressive-diffusive equation after the level set advection equation to maintain the interface thickness and profile. This equation was later modified by Olsson et al. (2007) to avoid tangential diffusion, which otherwise advects the interface numerically. In these studies of Olsson and Kreiss (2005) and Olsson et al. (2007), a hyperbolic tangent profile was advected instead of a signed normal distance function. Subsequently, Desjardins et al. (2008) improved the numerical oscillations of interface geometric properties by recomputing the signed normal distance function from this hyperbolic tangent function and used it within the framework of ghost fluid method for simulating turbulent atomization. Overall, these methods still require solving additional equations beyond the level set advection equation to conserve mass. An alternate strategy involves locally refining the mesh near the interface when solving the level set equations (Herrmann, 2008), which introduces a complex implementation procedure and a restrictive time-step.

Significant progress in numerical methodology can be observed in recent literature on compressive-diffusive term-based formulations for mass conservation, particularly within the framework of conservative phase field methods. Original phase field methods were solved using either the Cahn–Hilliard or Allen–Cahn equations. While the Cahn–Hilliard equation suffered from mass loss, unbounded phase field values at high density ratios (Yue et al., 2007), and coarsening effects (Soligo et al., 2019), the Allen–Cahn equation is in non-conservative form and therefore does not inherently preserve mass during interface evolution. Although the Allen–Cahn formulation is computationally attractive and has demonstrated improved robustness for high-density-ratio multiphase flows, its non-conservative nature can lead to cumulative mass errors. The phase field method was later improved into a conservative phase field method (CPF) (Chiu and Lin, 2011) by solving a second-order partial differential equation derived by removing the curvature-driven term from the Allen–Cahn equation (Sun and Beckermann, 2007). Recently, Mirjalili (2019) demonstrated that a central-difference scheme is sufficient to discretize these‘s partial differential equation’s to prevent mass loss and maintain the boundedness of the interface advection function between 0 and 1 (Mirjalili et al., 2020; 2022), and proposed this method as the conservative diffuse interface (CDI) method. A key feature of these modified phase field methods is that they employ a hyperbolic tangent function as the interface advection function and solve only a single partial differential equation for interface advection without requiring any additional reinitialization equation. Furthermore, Desjardins et al. (2008) also demonstrated the capability to recompute the signed normal distance function from this hyperbolic tangent function (Satheeshchandran et al., 2026). Motivated by these findings, in this study, the approach used for the interface advection equation in CPF, along with the boundedness criterion of CDI, is extended to a sharp interface–level set method on a co-located grid (SI − LSM col) (Chirammel et al., 2023b). In contrast to the conventional LSM, where reinitialization is periodically performed to restore the signed-distance property for accurate evaluation of interface normals and curvature, the additional reinitialization equation increases the computational cost and introduce cumulative numerical mass loss due to repeated modifications of the interface location. By incorporating interface regularization directly into the interface advection equation, these methods eliminate the need for a separate reinitialization procedure while maintaining a bounded interface representation.

The novelty of the present work lies in combining a conservative compressive-diffuse term-based interface advection with a sharp-interface level-set method employing jump conditions at the interface as a boundary condition on a co-located grid. The motivation for this combination is to eliminate the computationally expensive and mass-loss-inducing reinitialization equation required in the conventional level set method by replacing it with a conservative, bounded, and mass-preserving compressive-diffusive interface advection equation, while retaining the accuracy and robustness of the sharp-interface approach. Therefore, this study has two objectives: 1) to implement the compressive-diffuse term-based interface advection along with boundedness criterion within a sharp interface framework of the level set method and also 2) to demonstrate the capability of the current jump conditions at the interface, acting as an interfacial boundary condition, to avoid pressure-interfacial force imbalance on a three-dimensional co-located grid. To achieve the first objective, a previously implemented in-house solver, SI − LSM col, is modified by incorporating the compressive-diffusive term-based interface advection equation. This modified SI − LSM col is referred to as the sharp interface conservative level set method on co-located grid (SI − CLSM col). Unlike previous SI − LSM col studies of (Chirammel et al., 2023b; Chirammel et al., 2026; Chirammel et al., 2023a; Kothare et al., 2025) which used the conventional level-set formulation with reinitialization and were restricted to two-dimensional Cartesian and axisymmetric coordinate systems, the present work introduces the conservative level-set formulation into the sharp-interface framework and further extends the jump condition-based methodology to three-dimensional Cartesian grids. This enables, the demonstration that the proposed conservative sharp-interface approach preserves the pressure–interfacial force balance on three-dimensional co-located grids while simultaneously improving mass conservation.

This paper is organized as follows: Section 2 discusses the mathematical formulation associated with SI − LSM col and SI − CLSM col, along with a brief discussion of their numerical methodology in Section Equation 3. The overall solution algorithm is presented in Section Equation 4, followed by the results in Section 5. Finally, Section 6 concludes with the major findings of this study.

2 Mathematical formulation

Following the unified mathematical framework established in our previous studies (Chirammel et al., 2023b; Chirammel et al., 2026), the present work adopts the same strategy to clearly highlight the similarities and differences between the previous sharp-interface level set method and the current sharp-interface conservative level set method. The governing equations for mass and momentum conservation, along with their solution procedure, are identical to those used in the sharp-interface methodology presented in our previous work (Chirammel et al., 2023b; Chirammel et al., 2026). The distinction lies in the interface advection equation and the associated computations involving the level set function. The governing equations presented below are valid for incompressible, immiscible two-phase flows separated by a thin, massless interface.

Contiuity equation (incompressible flow):where U ⃗ denotes the dimensionless velocity vector.

∇ ⋅ U ⃗ = 0 (1)

Unified interface advection equation:

The evolution of a function representing the interface based on the fluid flow is given by

∂ ψ ∂ τ + ∇ ⋅ U ⃗ ∇ ψ = ∇ ⋅ R ⃗ (2)

Here, ψ is an interface function unified for SI − LSM col and SI − CLSM col, defined asand R ⃗ is a regularization term (Mirjalili et al., 2020),

ψ = φ for SI - CLSM col ϕ for SI - LSM col (3)

R ⃗ = Γ ̂ ϵ ∇ φ − φ 1 − φ ∇ φ | ∇ φ | for SI - CLSM col 0 for SI - LSM col (4)

For SI − CLSM col, the parameter ϵ controls the thickness of the transition layer over which φ varies smoothly from 0 to 1, while Γ ̂ regulates the strength of the regularization term (Mirjalili et al., 2022). For SI − LSM col, ϕ is a signed normal-distance function to the interface (Sussman et al., 1994). Furthermore, φ is initialized using the normal-distance function (Mirjalili et al., 2020) as

φ = 0.5 1 + tanh ϕ ϵ (5)

Conservation of momentum:where P, χ, η, and D = 0.5 ∇ U ⃗ + ∇ U ⃗ T denote the dimensionless pressure, density, viscosity, and rate-of-deformation tensor, respectively. The Reynolds and Froude numbers are defined as

∂ U ⃗ ∂ τ + ∇ ⋅ U ⃗ U ⃗ = − ∇ P χ i + 1 χ i R e ∇ ⋅ 2 η D − 1 F r 2 j ̂ . (6)

R e = ρ 1 u c l c μ 1 , F r = u c g l c , (7)

where l c and u c are the characteristic length and velocity scales, respectively, g is the acceleration due to gravity, and ρ and μ represent the density and dynamic viscosity. The subscript i identifies the fluid phase, with i = 1 corresponding to fluid-1 and i = 2 corresponding to fluid-2.

While solving the momentum equation, the surface forces on the interface are accounted for as a jump condition at the interface, given by

P = B X ⃗ = 2 R e η N ̂ ⋅ ∇ U ⋅ N ̂ , ∇ V ⋅ N ̂ + κ W e (8)

where [ P ] is the jump in non-dimensional pressure, and [ η ] is the jump in non-dimensional dynamic viscosity. The first and second terms on the right-hand side represent the viscous stress and surface tension along the normal direction at the interface. Furthermore, κ, N ̂, W e, and σ are the non-dimensional curvature, unit normal vector, Weber number, and surface tension coefficient, respectively, defined as

N ̂ = ∇ ϕ | ∇ ϕ | ; κ = − ∇ ⋅ N ̂ ; and W e = ρ 1 u c 2 l c σ

2.1 Subsidiary equations

In addition to the conservation equations mentioned above, the two interface-capturing methods differ in their computation of the signed normal-distance function. In the SI − LSM col approach, the transported level set function ϕ is reinitialized to restore its signed normal-distance character after advection. In contrast, SI − CLSM col advects a constrained hyperbolic tangent profile in a single step and employs the boundedness-preserving procedure of Mirjalili et al. (2020) to maintain interface sharpness without solving an additional, computationally expensive reinitialization equation. Consequently, reinitialization is required only in SI − LSM col, given by

Reinitialization equation:

∂ ϕ ∂ τ s = S ϵ ϕ 0 1 − | ∇ ϕ | (9)

where τ s is the pseudo time-step and S ϵ ( ϕ 0 ) is a smoothed sign function (Sussman et al., 1994). In the present work, the pseudo time-step is chosen as τ s = Δ X / 10, where Δ X is the grid size.

However, for SI − CLSM col, the conventional partial differential equation-based reinitialization equation (Equation 9) is replaced with a simpler algebraic relation between a signed normal-distance function Φ and φ (Jain, 2022), given byNote that although both ϕ and Φ correspond to signed normal-distance functions, ϕ is computed using (Equation 9) for SI − LSM col, whereas Φ is computed using (Equation 10) for SI − CLSM col.

Φ = ϵ × ln φ + ϵ 1 − φ + ϵ (10)

2.2 Transport-physical properties

Since the above conservation equations are solved within individual fluids along with jump conditions (Equation 8) as the boundary conditions at the interface, the transport-physical properties must be defined on the computational domain. This is accomplished using a unified signed normal-distance function ψ ̃ and a unified Heaviside function H ( ψ ̃ ) for both LSMs, given by

ψ ̃ = Φ for SI - CLSM col ϕ for SI - LSM col and H ψ ̃ = H γ Φ for SI - CLSM col H Γ ϕ for SI - LSM col (11)

Since a sharp interface method is considered in this study, the corresponding Heaviside functions for SI − CLSM col and SI − LSM col are denoted as H γ Φ and H Γ ϕ, respectively, and are defined as

H γ Φ = 1 Φ > 0 0 Φ ≤ 0 and H Γ ϕ = 1 ϕ > 0 0 ϕ ≤ 0 (12)

Using the unified Heaviside function, the transport-physical properties are defined as

χ i = H ψ ̃ + χ 1 − H ψ ̃ and η i = H ψ ̃ + η 1 − H ψ ̃ (13)

3 Numerical methodology

This section describes a common numerical framework used to solve the unified governing equations for both the SI − CLSM col and SI − LSM col formulations introduced in the previous section. The continuity and momentum equations are discretized using a physical law-based finite volume method (Sharma, 2017). This formulation follows the discrete conservation approach presented by Sharma (2017). For the sharp interface method, surface tension is omitted from the prediction equation and is instead treated as a boundary condition at the interface. The discretization of the governing equations results in a set of linear algebraic equations for the flow variables, which are advanced in time using a semi-explicit pressure-projection algorithm. On a co-located grid, in the predictor step, a provisional cell-center velocity U ̃ P is computed from the momentum equation without the pressure gradient, using only the advection and diffusion terms. The normal velocity predicted at the face-centers ( U e * ) is then obtained by linear interpolation of the provisional cell-center velocities of neighboring cells. These predicted face-center velocities are used to compute a predicted volume source ( S vol,P * ). In the corrector step, the pressure field is implicitly determined by enforcing mass conservation (Equation 1), leading to a pressure Poisson equation. The final velocities at the new time level ( U P n + 1 ) and the face-center normal velocities are corrected using the computed pressure gradients, with the coefficient β = Δ τ / χ i. For the momentum equation (Equation 6), we use a central difference scheme for the diffusion terms, and a third-order Lin–Lin TVD scheme (Date, 2005) for the advection terms, and a forward Euler for temporal integration. In contrast, the interface advection and reinitialization equations are discretized using the finite difference method. Further, the interface advection equation is discretized using a third-order Runge–Kutta scheme for the temporal integration and a fifth-order WENO scheme for the spatial term, while the regularization term in SI − CLSM col is discretized using a second-order central difference scheme. The interface advection equation employs higher-order discretization to accurately capture the interface while minimizing numerical diffusion and spurious oscillations. This approach avoids explicit surface tension discretization in the momentum equation as a source term, which is essential for maintaining pressure–interfacial force balance on a co-located grid. Since this formulation has been detailed in previous literature (Chirammel et al., 2023b), it is not repeated here; however, the differences from SI − LSM col are highlighted. For example, the pressure Poisson equation is obtained by subtracting the discretized velocity prediction equation from the momentum equation, leading to the following form (see details in Chirammel et al. (2023b); Shaikh et al. (2018)):

Δ τ ∇ ⋅ ∇ P n + 1 χ i = ∇ ⋅ U ⃗ * (14)

When discretized using the FVM, this pressure Poisson equation yields a linear algebraic equation for pressure, expressed as:

a P P P n + 1 = a W P W n + 1 + a S P S n + 1 + a E P E n + 1 + a N P N n + 1 + b (15)

The coefficients and the source term ( b ) in the above equation are given as

a W = β w Δ S w δ X w , a S = β s Δ S s δ Y s , a E = β e Δ S e δ X e , a N = β n Δ S n δ Y n ,

a P = ∑ N B = E , W , N , S a N B = a E + a W + a N + a S

b = U w * − U e * Δ Y + V s * − V n * Δ X + S jump,int

The only difference in the numerical methodology between the current LSMs when solving (Equation 15) lies in the S jump,int term (shown below in Equation 16). Specifically, the β and B terms are interpolated to the interface using Φ for SI − CLSM col and ϕ for SI − LSM col, represented using a unified signed normal-distance function ψ ̃ as shown below:with

S jump,int = ∑ f = e , w , n , s S jump , f = S jump , e + S jump , w + S jump , n + S jump , s (16)

S jump , e = β ̂ e B e Δ S e Δ X ; S jump , w = β ̂ w B w Δ S e Δ X S jump , n = β ̂ n B n Δ S n Δ Y ; S jump , s = β ̂ s B s Δ S s Δ Y

β ̂ f = β 1 β 2 β 1 1 − Ψ f + β 2 Ψ f ; Ψ f = | ψ ̃ | P | ψ ̃ | P + | ψ ̃ | N B ; β = Δ τ χ i

B f = B P | ψ ̃ N B | + B N B | ψ ̃ P | | ψ ̃ P | + | ψ ̃ N B |

Once ϕ associated with SI − LSM col is replaced with φ for SI − CLSM col, the numerical methodology for the momentum equation of SI − LSM col can be straightforwardly utilized for SI − CLSM col. Finally, the cell-centered velocity field is corrected using the actual pressure field (Equation 14), given as

U ⃗ ′ Δ τ = U ⃗ n + 1 − U ⃗ * Δ τ = − ∇ P n + 1 χ i (17)

where U ⃗ ′ is the corrected velocity and the superscript n denotes the time level.

For both LSM’s on co-located grid, the non-dimensional time step is chosen based on a criterion (Shaikh et al., 2018), given aswhere Δ τ CAP, Δ τ GFN, and Δ τ CFL are the time step due to capillary, grid Fourier number, and CFL, respectively.

Δ τ = m i n Δ τ CAP , Δ τ GFN , Δ τ CFL

4 Solution algorithm

The unified solution algorithm for both SI − CLSM col and SI − LSM col is summarized below. Steps 1–3 correspond to the pre-processing stage of the simulation, while Steps 4–13 constitute the solution procedure executed at each time step.

• Establish the computational domain by defining the Cartesian dimensions ( L × H ) together with the uniform grid spacings ( Δ X , Δ Y ).
• Prescribe the initial location of the interface through the interface function ϕ, and evaluate the associated geometric quantities, namely, the unit normal vector N ̂ and the curvature κ. Note that for SI − CLSM col, a signed normal-distance function ( ϕ ) is required for initialization. In addition, initialize φ using (Equation 5) for SI − CLSM col.
• Initialize the transport-physical properties (Equation 13) of the fluids, represented by ( χ i , η i ), along with the flow variables ( U ⃗ , P ) throughout the domain.
• Impose the appropriate boundary conditions for both the flow field variables and the interface function ψ.
• Evaluate the Heaviside function H ( ψ ̃ ) using (Equation 11) (expanded in (Equation 12)) for both SI − CLSM col and SI − LSM col.
• Determine the spatial distribution of the effective transport-physical properties of the two-phase system using H ( ψ ̃ ).
• Compute the predicted cell-centered velocity field by accounting for the advection and diffusion terms.
• Evaluate the predicted face-centered normal velocity at each control-volume face.
• Iteratively solve the pressure Poisson equation described by (Equation 15).
• Obtain the velocity correction by solving (Equation 17) and combine it with the intermediate velocity obtained in Step 7 to calculate the updated velocity field.
• Update the interface position by solving the interface advection equation (Equation 2) using the corrected velocity field.
• For SI − LSM col, perform the reinitialization procedure by solving (Equation 9). For SI − CLSM col, recompute the signed normal-distance function using (Equation 10).
• Evaluate the stopping criterion. If satisfied, terminate the simulation; otherwise, return to Step 4 for the next time step.

5 Validation and performance study

In this section, the proposed SI − CLSM col is validated using several test cases and compared with SI − LSM col to demonstrate the relative improvement in computational performance. These test cases include the Young–Laplace law to verify pressure–surface tension balance, the rising bubble and falling droplet test cases to assess the combined effects of surface tension and gravity, and the droplet coalescence to evaluate the capability of handling complex interface dynamics. For all test cases, the computational time is considered as the total time taken for the simulation from τ = 0 to τ = τ f, where τ f is the final non-dimensional time of the simulation. Note that, for the present simulations, the total computational time is compared because the capillary time-step is the most restrictive stability criterion for both interface-capturing methods, resulting in identical time-steps. Consequently, the total computational time provides a meaningful measure of computational efficiency. For simulations in which the time step is primarily dictated by the fluid velocity, the total number of time steps would also be an important metric for comparison. Mass loss is evaluated using the relative mass error = mass ini − mass n / mass ini, where mass ini is the mass at the initial instant of the simulation and mass n denotes the mass at a particular time instant. To maintain the boundedness of φ between 0 and 1, ϵ and Γ ̂ are selected as 1 and | U ⃗ | max, respectively (Mirjalili et al., 2020). To minimize the transition layer thickness while ensuring boundedness and stable computations, the smallest value of ϵ satisfying the boundedness criterion is recommended. This choice improves the estimation of interfacial geometric properties, which is crucial for the accuracy of the sharp-interface approach.

5.1 Young–Laplace law

To demonstrate that SI − CLSM col establishes pressure–interfacial force balance on a co-located grid similar to SI − LSM col, the Young–Laplace law test case is considered for various density ratios. In this problem, for a static and suspended spherical droplet under equilibrium, the non-dimensional pressure difference across the interface (Francois et al., 2006) is given bywhere P 1 , Γ and P 2 , Γ are the non-dimensional pressures in fluid-1 and fluid-2 across the interface, respectively.

Δ P exact = P 1 , Γ − P 2 , Γ = 2 κ W e (18)

To simulate this test case, the computational domain shown in Figure 1a is utilized. The domain consists of an initially static and suspended spherical droplet with an interface separating fluid-1 and fluid-2. The droplet has a non-dimensional radius of 2 and is placed at the center of a 3D Cartesian domain of size 8 × 8 × 8. For the fluid properties, various density ratios of 1, 1 0 3, and 1 0 5 are considered while keeping the viscosity ratio as 1. For the simulation of the non-dimensional conservation equations, the droplet diameter d is considered as l c ( = d ) and the capillary velocity as u c ( = σ / ( ρ 1 l c ) ). The flow fields are initialized to zero, and a free-slip condition is applied along the boundaries of the computational domain. The results presented here are obtained on a grid of size 40 × 40 × 40. To evaluate the accuracy of the solvers, two error metrics are considered: the maximum magnitude of spurious currents E ( U ⃗ max ) and the maximum relative error in pressure E ( Δ P ) max, defined aswhere

E U ⃗ max = max | U ⃗ | and E Δ P max = | Δ P max − Δ P exact | Δ P exact (19)

Δ P max = P max − P min (20)

Table 1 shows that the computational accuracy of both LSMs is similar, demonstrating the capability of the present SI − CLSM col to avoid pressure–interfacial force imbalance on a co-located grid for high-density-ratio simulations. As in the previous SI − LSM col, the pressure distribution in the computational domain maintains a sharp jump across the interface for SI − CLSM col. The maximum magnitude of spurious currents, E ( U ⃗ max ), is also nearly identical for both methods over the entire range of density ratios considered. This is expected because spurious currents primarily arise from the discrete imbalance between the surface-tension force and the pressure gradient, as well as errors in the computation of interface geometric quantities such as the interface normal and curvature. Since SI − CLSM col employs the same sharp-interface formulation, surface-tension discretization, pressure solver, and geometric property evaluation as SI − LSM col, the conservative transport equation does not directly alter the numerical mechanisms responsible for spurious-current generation. Instead, the principal advantage of SI − CLSM col lies in its conservative interface transport, which eliminates the need for the separate reinitialization equation required in SI − LSM col. Consequently, numerical mass loss associated with repeated reinitialization is significantly reduced while maintaining spurious currents at a level comparable to the SI − LSM col.

TABLE 1. Error metrics obtained using SI − CLSM col and SI − LSM col for different density ratios in the Young–Laplace law test case.

All errors are evaluated after 1000 iterations.

A grid refinement study performed on grid sizes of 40 × 40 × 40 , 60 × 60 × 60 , and 80 × 80 × 80 further demonstrates that these improvements are retained with increasing mesh resolution. While both methods exhibit improved mass conservation upon mesh refinement, SI − CLSM col achieves 2–3 order of magnitude of lower mass loss than SI − LSM col, together with a modest reduction in computational time owing to the elimination of the separate reinitialization equation. Furthermore, as shown in figure Figure 2, the grid convergence study based on the maximum magnitude of spurious current indicates that both formulations exhibit convergence rates between first and second order with grid refinement.

5.2 Rising bubble

A buoyancy-induced rising bubble in a three-dimensional Cartesian coordinate system is simulated in this section. Figure 1b shows an initially spherical bubble positioned at a distance of (1, 1, 1) in a computational domain of size 2 × 4 × 2. The non-dimensional properties of the fluids are shown in Table 2. These properties for the rising bubble case are obtained by considering the initial bubble diameter as the length scale and g l c as the velocity scale. This simulation is performed with a non-dimensional time-step of 1 0 − 4 on a grid of size 80 × 160 × 80.

TABLE 2. Non-dimensional properties used for the simulation of various test cases.

Figure 3 shows good agreement between the rising bubble velocities obtained using the LSMs in this study and previously published numerical results (Turek et al., 2019; Adelsberger et al., 2014). The transient evolution of the rise velocity is within ± 5 % of the results of Turek et al. (2019). The data pertaining the computational time and the mass conservation error obtained for this case, and the following ones, are reported in Figures 4, 5, respectively. Regarding the rising bubble case on a finer grid, the computational time and mass loss for SI − CLSM col are 17.17 h r and 0.048 %, respectively, while for SI − LSM col, they are 34.90 h r and 2.14 %, respectively. This represents an improvement factor of 2.03 in computational time and 44.71 in mass loss compared to SI − LSM col. Note that the SI − LSM col has been previously validated for several 2D benchmark cases (Chirammel et al., 2023c), including the rising bubble benchmark of Hysing et al. (2009).

5.3 Falling droplet

For the falling droplet test case, Figure 1c shows the three-dimensional computational setup used for performance evaluation, while a similar two-dimensional setup is employed for validation. The test case consists of a water droplet with an initial non-dimensional radius of 0.5 positioned at (2, 4) within a rectangular computational domain of non-dimensional size 4 × 6. The non-dimensionalization follows the same approach as in the rising bubble case, where the initial droplet diameter is taken as the characteristic length scale l c and the velocity scale is defined as u c = g l c. This yields the non-dimensional properties listed in Table 2. The simulation is carried out using a non-dimensional time-step of 1 0 − 4 on a uniform grid of 80 × 120. The validation study is performed in a two-dimensional Cartesian framework, whereas the computational performance assessment is additionally conducted in a three-dimensional domain of size 4 × 6 × 4 on a grid of size 80 × 120 × 80.

Figure 6 shows the validation for the falling droplet. Good agreement is obtained when comparing the temporal evolution of the droplet centroid with the numerical results obtained using LSM (Kang et al., 2000) and CLSVOF (van der Pijl et al., 2005), demonstrating the capability of the present SI − CLSM col to handle surface-tension-dominated problems. As shown in Figures 4, 5, SI − CLSM col achieves improvements in computational time and mass loss over SI − LSM col by factors of 3.25 and 7.66 h r, respectively. On a sufficiently finer grid size to capture the centroid motion, the computational time and mass loss for SI − CLSM col are 7.66 h r and 0.194 %, respectively, while for SI − LSM col, they are 24.90 h r and 2.14 %, respectively.

5.4 Droplet coalescence

In this test case, an initially stationary droplet coalesces with a quiescent liquid pool. Figure 1d shows the computational setup for simulating this problem. This surface-tension-dominated problem leads to significant interfacial deformation, thereby providing a stringent assessment of the solver’s ability to capture dynamically evolving interfaces. A droplet of non-dimensional diameter 1 is initially placed above a liquid pool of size 4.0 × 2.5 × 4.0, requiring a computational domain of 4 × 8 × 4. The initial droplet diameter D is used as l c, and the capillary velocity σ / ( ρ 1 l c ) is taken as u c. This non-dimensionalization results in the properties listed in Table 2. The simulation is performed using a non-dimensional time-step of 1 0 − 4 on a uniform grid of 200 × 400 × 200.

Figure 7 shows the temporal evolution of the non-dimensional droplet height compared with the experimental results of Blanchette and Bigioni (2006) and the numerical results obtained using the LSM reported by Shaikh et al. (2018). It can be observed that both LSMs provide good predictions of the interface evolution, as both are within ± 5 % of the experimental data. The computational time and mass loss are improved by factors of 1.66 and 4.28, respectively, with SI − CLSM col compared to SI − LSM col. On a sufficiently finer grid size to capture the coalescence dynamics, the computational time and mass loss for SI − CLSM col are 57.28 h r and 0.007 %, respectively, while for SI − LSM col, they are 94.90 h r and 0.03 %, respectively.

6 Conclusion

A previous in-house solver—the sharp interface level set method on a co-located grid, SI − LSM col—has been extended to improve mass conservation and computational performance as the sharp interface conservative level set method on a co-located grid, SI − CLSM col. From the comparative study on several three-dimensional benchmark cases, the following conclusions can be drawn:

• For sharp interface simulations, a signed normal distance function recomputed from a hyperbolic tangent function, when properly bounded, can be utilized for simulating immiscible two-phase flows with high density ratios on a co-located grid.
• As demonstrated in the numerical methodology, the present approach is easy to implement in an existing sharp interface method because the source terms due to the jump condition remain unchanged.
• The Young–Laplace law test case demonstrates that the accuracy of both LSMs is similar for various density ratios of 1, 1 0 3, and 1 0 5.
• Compared to SI − LSM col, SI − CLSM col improves mass conservation and reduces computational time for all test cases considered—namely, the Young–Laplace law, rising bubble, falling droplet, and droplet coalescence.
• The jump conditions at the interface, acting as an interfacial boundary condition, are robust for both LSMs and on three-dimensional coordinate systems for high density ratios on a co-located grid.

Overall, a solver that is straightforward to implement without the additional complexities associated with hybrid approaches is proposed for improved mass conservation and reduced computational time on co-located grids for two-phase flows with high density ratios.

Generative AI statement

The author(s) declared that generative AI was used in the creation of this manuscript. During the preparation of this manuscript, the authors used ChatGPT to refine the language of the text. After using this tool, the authors carefully reviewed and edited the content as necessary and take full responsibility for the accuracy and integrity of the published work.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Omid Mahian, Xi’an Jiaotong University, China

Reviewed by: Martin Van Der Eijk, Deltares, Netherlands

Deepak Kumar Singh, Dr DY Patil Institute of Technology Department of Mechanical Engineering, India