💻 computer science

Topology-Preserving Mesh Adaptation for Sharp-Interface Multiphase PFEM

This paper introduces a robust, fully Lagrangian Particle Finite Element Method (PFEM) framework that employs a dynamic mesh adaptation strategy to preserve sharp interfaces and decouple topological changes from grid resolution, enabling accurate and scalable simulations of multiphase flows with an arbitrary number of immiscible phases.

Original authors: Félix Ruyffelaere, Michel Henry, Jonathan Lambrechts, Jean-François Remacle

Published 2026-07-13

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Félix Ruyffelaere, Michel Henry, Jonathan Lambrechts, Jean-François Remacle

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are trying to film a chaotic dance party where different groups of dancers (fluids) are swirling, merging, and splitting apart. In the world of computer simulations, these dancers are usually represented by a grid of tiny squares or triangles, like a digital mosaic. The problem is, when these dancers move wildly, the mosaic gets stretched, twisted, and eventually turns into a messy, unrecognizable blob.

For decades, scientists have tried to fix this. Some methods use a "fuzzy" approach, blurring the lines between groups so the grid doesn't break, but this makes the boundaries soft and inaccurate. Others try to keep the lines sharp but often get stuck because the grid itself forces the dancers to merge or split just because the tiles are too big or too small. It's like trying to cut a piece of paper with scissors that are too big; you end up cutting things you didn't mean to.

This paper introduces a clever new way to handle these "dance parties" using a method called the Particle Finite Element Method (PFEM). Think of the PFEM as a magical camera that doesn't just take a picture; it actually rebuilds the entire stage every single second.

The Magic Trick: The "Protecting Ball"
The core of this new method is a rule the authors call the "protecting-ball criterion." Imagine every edge of a dancer's boundary is holding a transparent, empty bubble (a disk). The rule is simple: as long as that bubble stays empty of other dancers, the edge must stay there.

In the past, if the grid got messy, the computer would just redraw the lines however it felt like, often accidentally snapping two groups together or tearing them apart. This new method says, "No way!" Before the computer redraws the grid, it checks every boundary edge. If an edge is about to get squished or flipped, the computer does two things:

Splits the edge: It chops the long edge into smaller pieces, making the "protecting bubble" smaller and easier to keep empty.
Filters the crowd: It politely asks any stray dancers (nodes) that are hovering inside those bubbles to step aside.

Because of this, the computer can rebuild the grid using a standard, fast mathematical recipe (Delaunay triangulation) and still guarantee that the sharp lines between the different fluids stay exactly where they belong. It's like having a bouncer who ensures that even if the crowd moves, the VIP section (the interface) remains perfectly intact.

The Results: Simulating the Chaos
The authors tested this idea with some wild scenarios to see if it holds up.

The Vortex Test: They simulated two droplets of liquid spinning in a box. The droplets stretched into incredibly thin, hair-like threads. The new method kept these threads sharp and didn't let them snap or merge by accident. When the spin reversed, the threads retracted, and the droplets returned to their original round shapes, proving the method doesn't lose "mass" (liquid) along the way.
The 3-Phase Rain: They simulated three layers of fluids with different densities falling under gravity (a Rayleigh-Taylor instability). The fluids mixed and formed complex, finger-like structures. The simulation matched other famous methods very closely, but with a trick: it used far fewer total points (nodes) because it only added extra detail right where the mixing happened, keeping the rest of the simulation light and fast.
The Rising Bubble: They watched a bubble rise through a liquid. The bubble stretched and squished, changing shape dramatically. The new method predicted the bubble's speed and shape almost perfectly, matching high-end reference models but using significantly fewer computational resources.

The "What If" and the "Not Yet"
The paper is very clear about what it doesn't do yet.

No Magic for Triple Points: The method currently handles two fluids perfectly. However, when three different fluids meet at a single point (a "triple point"), the math gets tricky. The authors admit they haven't solved this yet, so for now, their simulations with more than two fluids (like a 16-phase experiment they ran just to show off the geometry) ignore surface tension forces.
Grid Size Still Matters (A Little): While the method separates the physics from the grid size better than before, the authors show that if you force a "cut" in a thin filament too early (by flipping an edge), the simulation changes completely. This proves that the underlying physics of when to break a thread is still a mystery that needs a specific physical model, not just a grid rule.
2D Only (For Now): All the tests in this paper were done in two dimensions (flat, like a drawing). The authors are confident the math works in 3D (real life), but they haven't run those specific 3D tests yet.

The Bottom Line
This paper doesn't claim to have solved every problem in fluid dynamics. Instead, it offers a robust, topology-preserving framework that stops the computer grid from accidentally changing the story of the fluids. By using "protecting balls" to guard the boundaries, the method keeps the interfaces sharp and the physics honest, allowing for complex simulations of mixing, splitting, and swirling fluids with much less computer power than before. It's a solid step forward, proving that with the right geometric rules, we can simulate the messy dance of fluids without losing our minds—or our data.

Technical Summary: Topology-Preserving Mesh Adaptation for Sharp-Interface Multiphase PFEM

Problem Statement
The accurate numerical simulation of immiscible multiphase flows presents significant challenges due to the multi-scale disparity between macroscopic convective transport and microscopic interfacial phenomena. Traditional Eulerian methods, such as Volume-of-Fluid (VOF), Level-Set, and Diffuse Interface models, often suffer from inherent limitations: Diffuse Interface methods introduce non-physical numerical diffusion and prevent strict enforcement of sharp pressure jumps; Level-Set methods are prone to numerical mass loss; and VOF methods require complex geometric reconstruction for accurate curvature. Conversely, while Lagrangian approaches like the Particle Finite Element Method (PFEM) avoid convective terms and handle large deformations well, they struggle with maintaining mesh integrity during extreme nodal advection. Standard PFEM remeshing often relies on constrained triangulations or $\alpha$ -shape filters, which can inadvertently alter interface connectivity through edge flips. This leads to topological changes dictated by mesh resolution rather than underlying physics, potentially introducing non-physical artifacts such as premature filament pinch-off or artificial satellite bubble generation.

Methodology
The authors propose a fully Lagrangian framework based on PFEM that integrates a robust, topology-preserving mesh adaptation strategy. The core of the method is a remeshing procedure that decouples interface geometry preservation from the global mesh size, ensuring that sharp interfaces are maintained without relying on classical constrained triangulation.

Governing Equations and Advection: The incompressible Navier-Stokes equations are solved using a stabilized $P1-P1$ Finite Element Method (FEM). To handle surface tension in multi-phase systems, an enriched pressure field ( $p^*$ ) is introduced to capture discontinuities at interfaces, though this specific formulation is currently limited to two-phase flows. The Lagrangian framework eliminates the convective acceleration term. To improve mass conservation, the authors utilize a second-order explicit Adams-Bashforth (AB2) scheme for node advection, leveraging temporal velocity history to avoid the topological constraints associated with higher-order Runge-Kutta sub-steps.
Protecting-Ball Criterion (Theoretical Foundation): The method is grounded in the duality between Voronoï diagrams and Delaunay triangulations. The authors establish a "Protecting-ball criterion" (Theorem 1): if every segment of a discretized interface is associated with an open disk containing no other vertices (an empty disk), that segment is guaranteed to appear in the unconstrained Delaunay triangulation of the point cloud.
Interface Recovery and Adaptation: To satisfy the protecting-ball criterion during remeshing, the algorithm performs two sequential operations on the advected point cloud:
- Edge Splitting: Interface edges that are too long or lack an associated empty disk are recursively split to reduce the required disk diameter.
- Volume Node Filtering: Volume nodes (those not belonging to the interface) that encroach upon the empty disks of interface segments are removed. Crucially, interface nodes themselves are never removed, preserving the discrete topology.
- Coarsening: To prevent unbounded node growth, over-refined interface regions are coarsened based on a geometric flatness tolerance, provided the removal does not alter the topology of triple points or violate the minimum size threshold ( $w_{min}$ ).
Bulk Mesh Adaptation: Following interface recovery, a refinement algorithm based on Chew's work is applied to the fluid bulk. This ensures a smooth transition in element size from the highly refined interface region to the coarser bulk, governed by a user-defined sizing field. This step optimizes element quality (aspect ratio) while respecting the protected interface disks.
Data Transfer: Because the new mesh is conforming to the previous interface topology, physical quantities (phase indicators and velocity) are transferred by querying the centroid of new elements against the previous mesh, ensuring mass conservation and phase assignment.

Key Contributions

Topology Preservation without Constrained Triangulation: The paper introduces a method that guarantees the recovery of sharp interfaces in a standard Delaunay triangulation by enforcing a diametral empty-disk condition through edge splitting and node filtering. This eliminates the computational overhead and 3D generalization issues associated with constrained Delaunay triangulations.
Decoupling Physics from Mesh Resolution: The framework separates the decision of topological transitions (merging, pinching) from the mesh size. Transitions are only permitted when the interface segments reach a critical minimum size ( $w_{min}$ ), below which the empty-disk criterion is relaxed. This allows for the integration of sub-grid physical models to govern topological changes independently of the discretization.
Robust Mass Conservation: By employing a second-order Adams-Bashforth advection scheme and a conformal remeshing strategy, the method significantly reduces mass loss compared to first-order explicit Euler schemes, as validated against benchmark cases.
Scalability to Multiphase Systems: The method is demonstrated to handle an arbitrary number of immiscible phases, successfully simulating complex scenarios including a 16-phase Rayleigh-Taylor instability.

Results and Validation
The framework was validated against standard multiphase benchmarks:

Vortex in a Box: The method successfully preserved thin fluid filaments and recovered sharp boundaries through successive remeshing steps. The use of the Adams-Bashforth scheme demonstrated superior mass retention compared to first-order schemes, with results closely matching high-order Runge-Kutta references.
Three-Phase Rayleigh-Taylor Instability: The simulation results showed strong agreement with Volume-of-Fluid (VOF) and Moving Particle Semi-implicit (MPS) benchmarks. The adaptive meshing significantly reduced the total node count compared to uniform grids by clustering nodes near interfaces, while still capturing thin filaments during the mixing phase.
Isolated Rising Bubble: The method accurately predicted bubble rise velocity and circularity, matching reference data from Level-Set and ALE solvers. It achieved these results with significantly fewer nodes than classical diffuse-interface approaches.
Mesh Dependence and 16-Phase Flow: The authors demonstrated that simulation outcomes are highly sensitive to mesh configuration. By artificially introducing an edge flip, they showed how minor mesh variations could trigger different topological outcomes (e.g., satellite bubble formation). The method's ability to handle a 16-phase system confirmed its geometric versatility and robustness.

Significance and Claims
The paper claims that this framework provides a robust solution for simulating immiscible multiphase flows by rigorously preserving geometric interfaces and physical topology during the remeshing process. The primary significance lies in the ability to maintain sharp interfaces without the numerical diffusion of diffuse-interface methods or the mass loss of Level-Set methods, all while avoiding the computational complexity of constrained triangulations.

The authors emphasize that their approach allows physical models, rather than discretization artifacts, to dictate topological changes. By preserving the interface topology down to a specific tolerance ( $w_{min}$ ), the method prevents non-physical phenomena such as premature filament pinch-off. While the current implementation is restricted to two dimensions and two-phase surface tension modeling, the authors assert that the underlying geometric philosophy and theoretical guarantees are fully valid in three dimensions, representing a logical next step for the framework. The work highlights that while mesh adaptation improves efficiency, the coupling between mesh resolution and physical outcomes remains critical, necessitating careful selection of resolution parameters or the integration of sub-grid models for topological transitions.

Technical Summary: Topology-Preserving Mesh Adaptation for Sharp-Interface Multiphase PFEM

More like this