Geometric Reinitialization for Capillary Flows: a Comparative Study with State-of-the-Art Conservative Level-Set Methods

Abstract: Simulations of immiscible flows involving surface tension (ST) require a robust high-fidelity framework. State-of-the-art multi-phase models, such as the Conservative Level-Set (CLS) approach, rely on Eulerian representations of the fluids and interface and require reinitialization methods to ensure volume conservation and accurate ST force modeling. This work focuses on the complete description of a CLS solver and proposes a novel geometric reinitialization method, based on the level-set literature. It includes a quantitative and objective comparison of this new geometric method to two reinitialization approaches: the PDE-based reinitialization proposed in the original CLS method and a simple projection-based approach. This comparison tackles three 3D application cases: the rise of a bubble, the capillary migration of a droplet, and the Rayleigh-Plateau instability development in a capillary jet. The PDE-based and geometric methods lead to high-quality, spatially-converged results in good agreement with benchmark and analytic solutions, while the projection-based reinitialization fails to capture complex 3D interfacial dynamics. The results also highlight the robustness of the novel geometric method which offers a two-parameter framework in comparison to the PDE-based method that necessitates a case-dependent selection of four parameters.

• The paper compares three reinitialization strategies for Conservative Level-Set methods within a consistent FEM framework, determining that a novel geometric method outperforms traditional PDE- based and projection-based approaches in handling capillary flows.
• The geometric reinitialization method provides better volume and interface position preservation at coarser resolutions, demonstrating up to 1.15% interface volume fluctuation over the PDE-based result of 1.93%.
• The study underscores geometric method's superior performance for conservative physics and plots method parameter tuning sensitivity for potential accuracy degradation in capillary migration and jet breakup cases.

Overview and motivation

This paper presents a complete Conservative Level-Set (CLS) two-phase flow solver within a finite element method (FEM) framework, implemented in the open-source code Lethe built on deal.II, and uses it to conduct a controlled comparison of three interface reinitialization strategies: the PDE-based reinitialization of Olsson and Kreiss (the original CLS approach), a projection-based sharpening scheme following Aliabadi et al., and a novel geometric reinitialization method adapted from level-set redistanciation techniques of Mut et al. and Ausas et al. The central claim is that reinitialization quality—not merely advection accuracy—governs the fidelity of capillary-driven simulations, and that this can only be assessed objectively when all methods are embedded in an identical numerical framework, a comparison the authors identify as missing from the literature.

The solver treats the one-fluid incompressible Navier–Stokes equations with PSPG/SUPG stabilization on equal-order $Q_1$ elements, BDF2 time integration, explicit coupling between phase-indicator transport and momentum, a Continuous Surface Force model supporting Marangoni stresses via a linear surface tension coefficient $\gamma(x)$, and $\mathcal{L}^2$ projections of the normal and curvature following Zahedi et al. Dynamic mesh adaptation based on Kelly's error estimator refines around the interface.

The three reinitialization methods

The PDE-based method solves a compressive–diffusive artificial-time PDE to steady state, yielding a hyperbolic tangent profile $\phi = 0.5 - 0.5\tanh(d/2\varepsilon)$; it requires four parameters ($\varepsilon$, $\Delta\tau$, a steady-state tolerance, and a maximum iteration count). The projection-based method applies a piecewise power-law sharpening of $\phi$ toward the iso-level $c=0.5$ with exponent $\alpha=1.5$. The proposed geometric method reconstructs the interface via marching cubes, computes exact signed distances from the reconstructed geometry using the discrete counterpart of the distance property $|d(y)| = \min_{x\in\mathcal{C}}(\|y-x\| + |d(x)|)$, applies local and global volume corrections, and converts the distance field back to a CLS profile via the same tanh transformation. Its extensions over prior work are support for quadrilateral/hexahedral elements, hanging-node constraints for adaptively refined grids, and distributed parallel execution. Notably, both the PDE-based and geometric methods produce the same target profile, differing only in how the signed distance is obtained—yet they require different values of $\gamma(x)$0 to achieve their best results, a discrepancy the authors explicitly leave unexplained.

Rising bubble benchmark

The 3D rising bubble case (Re = 35, Eo = 10, following Adelsberger et al.) is evaluated against reference solutions from DROPS, NaSt3DGPF, OpenFOAM, and FeatFlow across reinitialization frequencies $\gamma(x)$1. The PDE-based and geometric methods agree well with DROPS and NaSt3DGPF for barycenter height and rise velocity, with negligible sensitivity to frequency. The volume enclosed by the interface is conserved to within roughly 1–2% loss for these two methods, though the PDE-based results exhibit jumps in sphericity at low reinitialization frequency, indicating artificial displacement of the interface at each reinitialization event—a known weakness of the original CLS formulation that the geometric method avoids, producing smooth sphericity evolution in closer agreement with references.

The projection-based method fails decisively: the solver diverges at $\gamma(x)$2, and even at the lowest frequency the final bubble volume drops to approximately 92% of its initial value (75% at $\gamma(x)$3), with staircase-shaped contours locked to the Cartesian mesh and nonphysical sphericity oscillations. This constitutes the paper's strongest negative result: a simple projection cannot sustain 3D capillary dynamics in this framework.

Capillary migration

The droplet migration under a linear surface tension gradient admits the analytical steady-state velocity $\gamma(x)$4, making it a stringent test since motion arises purely from tangential interfacial forces. The PDE-based and geometric methods recover the analytical migration velocity at all frequencies, with velocity dips between reinitialization events at $\gamma(x)$5 underscoring how diffusion of the phase indicator directly degrades Marangoni-driven transport. A limitation emerges for the geometric method: at the highest frequency, where the reinitialization CFL drops to roughly 0.01, the final radius distribution develops sharp, front–back asymmetric oscillations, whereas the PDE-based method remains unaffected. The projection-based method fails entirely, producing strong radius oscillations and no migration.

Rayleigh-Plateau instability

Spatial convergence is assessed on four adaptive mesh levels (up to ~4.6×10⁷ cells) at fixed CFL = 0.25, monitoring relative jet volume and surface area through breakup. Both surviving methods converge toward similar solutions, but the geometric method converges faster: its relative volume oscillates around unity before breakup, while the PDE-based method exhibits a spurious volume increase in the first regime and slower convergence after droplet pinch-off. For both methods, convergence of the relative surface area degrades after necking begins (~0.02 s), indicating that post-breakup droplet shapes remain discretization-sensitive even when volume is converged. Two case-specific adjustments are disclosed: the Heaviside regularization is deactivated for the PDE-based method (using $\gamma(x)$6 directly), without which volume conservation and spatial convergence deteriorate, and parasitic capillary waves corrupt the coarse-mesh result.

Parameter sensitivity and robustness

A recurring theme is parameter burden. The PDE-based method requires four coupled parameters whose selection is case-dependent: $\gamma(x)$7 works for the rising bubble but $\gamma(x)$8 is needed for capillary migration, and higher values cause either poor artificial-time convergence or nonphysical oscillations. The geometric method uses only two parameters ($\gamma(x)$9, $\mathcal{L}^2$0), and the same settings ($\mathcal{L}^2$1, $\mathcal{L}^2$2) transfer across the rising bubble and capillary migration cases, with $\mathcal{L}^2$3 having no significant effect on solution quality beyond added cost. The paper also argues, implicitly through its metric selection, that global quantities such as barycenter position alone are insufficient discriminators; shape-sensitive metrics (sphericity, radius distribution, enclosed volume, surface area) are required to expose reinitialization artifacts.

Limitations and open questions

The authors concede several points. The reason the optimal $\mathcal{L}^2$4 differs between the PDE-based and geometric methods despite identical target profiles is unresolved. The geometric method's degradation at low reinitialization CFL (observed in capillary migration) is documented but not analyzed mechanistically. The Rayleigh-Plateau study fixes a single reinitialization frequency and modifies the Heaviside treatment specifically for the PDE-based method, so the convergence comparison is not fully symmetric between methods. The absence of theoretical error bounds for regularized discontinuous coefficients in full two-phase Navier–Stokes problems means the choice of resolvable interface thickness rests on empirical evidence rather than analysis. Whether the geometric method retains its robustness under extreme deformation regimes—such as laser powder bed fusion with large evaporation- and Marangoni-induced pressure jumps—is posed as an open application question rather than demonstrated.

Conclusion

Within a single FEM/CLS framework, this study establishes that PDE-based and geometric reinitialization both deliver accurate, spatially converged solutions for 3D capillary flows, that the geometric variant better preserves interface position and volume at coarser resolutions, and that simple projection-based sharpening is inadequate for realistic 3D interfacial dynamics. The practical contribution is a demonstrably more robust two-parameter geometric alternative to the four-parameter, case-tuned PDE-based standard, validated against benchmarks, analytical solutions, and a convergence study on jet breakup.