Source-linked AI summary
Accurate and efficient surface reconstruction from volume fraction data on general meshes
Henning Scheufler, Johan Roenby
TL;DR
The paper addresses noisy or inaccurate interface reconstruction in geometric VOF methods, particularly on unstructured meshes. It introduces iterative reconstructed-distance-function schemes and combines them with isoAdvector advection. The schemes achieve second-order reconstruction convergence and, for CFL 0.2 or below, second-order advection convergence across the tested mesh types.
Problem
Existing isoAdvector isosurface reconstruction can produce noisy interface orientations, especially on unstructured meshes, while accurate reconstruction is needed for interface advection.
Method
The paper iteratively reconstructs a distance function and interface geometry, using RDF-based schemes with isoAdvector and previous-time-step normal information for initialization.
Results
Second-order convergence is obtained for interface position and orientation on structured and unstructured 2D and 3D meshes, while advection achieves second-order convergence for CFL 0.2 or below.
Takeaways & Limitations
The reconstruction schemes provide accurate general-mesh interface geometry with straightforward implementation and reduced computational cost relative to contemporary methods.
Takeaways & Limitations
At maximum deformation, unresolved thin interface tails can produce pinched-off droplets, particularly on triangular prism meshes.
Abstract
from arXiv · showhide
Simulations involving free surfaces and fluid interfaces are important in many areas of engineering. There is, however, still a need for improved simulation methods. Recently, a new efficient geometric VOF method called isoAdvector for general polyhedral meshes was published. We investigate the interface reconstruction step of isoAdvector, and demonstrate that especially for unstructured meshes the applied isosurface based approach can lead to noisy interface orientations. We then introduce a novel computational interface reconstruction scheme based on calculation of a reconstructed distance function (RDF). By iterating over the RDF calculation and interface reconstruction, we obtain second order convergence of both the interface normal and position within cells even with a strict $L_{\infty}$ error norm. In 2D this is verified with reconstruction of a circle on Cartesian meshes and on unstructured triangular and polygonal prism meshes. In 3D the second order convergence is verified with reconstruction of a sphere on Cartesian meshes and on unstructured tetrahedral and polyhedral meshes. The new scheme is combined with the interface advection step of the isoAdvector algorithm. Significantly reduced absolute advection errors are obtained, and for CFL number 0.2 and below we demonstrate second order convergence on all the mentioned mesh types in 2D and 3D. The implementation of the proposed interface reconstruction schemes is straightforward and the computational cost is significantly reduced compared to contemporary methods. The schemes are implemented as an extension to the Computational Fluid Dynamics (CFD) Open Source software package, OpenFOAM. The extension module and all test cases presented in this paper are released as open source.
1. Introduction
Free-surface simulations need accurate interface modeling, yet existing VOF approaches can be insufficiently accurate, especially for interface normals on unstructured meshes. The paper develops iterative RDF-based reconstruction methods targeting efficient, accurate general-mesh simulation.
- Free-surface simulations support quantitative and qualitative engineering analysis, but their numerical difficulty has sustained research activity.
- VOF and level-set methods remain widely used, but are often not sufficiently accurate.
- PLIC represents each cell interface with a plane, while Youngs is easy to implement on unstructured meshes but lacks accuracy.
- The proposed methods are iterative RDF variants focused on interface reconstruction for unstructured three-dimensional general polyhedral meshes.
- Using point-neighbour information, residual stopping, and previous-time-step normals supports efficient implementation and second-order convergence with mesh refinement.
2. The fundamental VOF equation
The VOF formulation represents a sharp two-fluid interface implicitly through cell volume fractions, but those fractions do not uniquely determine the within-cell interface configuration. Reconstruction is therefore needed to recover interface geometry from neighboring volume-fraction data.
- The system consists of two immiscible, incompressible fluids described by velocity, pressure, and a dynamically evolving sharp interface.
- The density field is expressed using a three-dimensional Heaviside function distinguishing fluids A and B.
- VOF uses cell volume fractions α_i to implicitly represent the fluid interface and discretize its evolution equation.
- The same volume fractions can arise from different fluid configurations, leaving the interface distribution at cell faces unknown.
- Assuming a smooth interface supplies a local curvature length scale that constrains reconstruction when the cell size is sufficiently small.
3. Interface reconstruction schemes
The reconstruction schemes recover a planar interface position and normal from volume-fraction data in polyhedral cells. The iso-Alpha method uses interpolated isosurfaces and geometric volume matching, while iso-RDF iteratively improves the reconstruction to address poor orientation convergence, especially on unstructured meshes.
- Interface reconstruction provides the interface centre position x_S and unit normal n̂_S in cells with 0 < α_i < 1.
- iso-Alpha method: iso-Alpha interpolates cell volume fractions to vertices, constructs an internal isosurface, and uses it to split each interface cell.
- iso-Alpha method: The isoface is formed by connecting edge cut points, then its geometry and the submerged cell volume are calculated from polygonal faces.
- iso-Alpha method: The isosurface value f* is found so the submerged subcell volume fraction matches α_i within the specified tolerance.
- iso-RDF method: iso-Alpha converges slowly on structured meshes and not at all on unstructured meshes, motivating iso-RDF.
- iso-RDF method: iso-RDF initializes from iso-Alpha, calculates a reconstructed distance function, and iterates to improve x_S and n̂_S.
- iso-RDF method: Least-squares interpolation is used in iso-RDF, including for vertex volume fractions and the initial reconstruction weights.
15 end
The RDF approach replaces the noisy derivative behavior of sharp VOF data with a smoother distance-like field for interface-normal reconstruction. Its iterative PLIC-RDF scheme combines reconstructed-distance gradients, optional temporal initialization, and curvature-aware convergence checks.
- A sharp VOF field changes from α = 0 to α = 1 over one cell, making derivatives difficult for curvature calculations.
- A reconstructed distance function is smoother than raw VOF data but need not have a signed-distance gradient length exactly equal to 1.
- The RDF in each cell is computed by a weighted average of distances from the cell and point-neighbor interfaces.
- The weighting exponent is set to A = 2 because higher values can cause divergence in under-resolved regions, although exponent choice does not affect accuracy for well-resolved interfaces.
- RDF values at vertices are obtained from surrounding cell values using a least-squares reconstruction.
- PLIC-RDF method: PLIC-RDF replaces isosurface calculations in normal estimation with gradient calculations and can initialize normals from the previous time step.
- PLIC-RDF method: Curvature-aware residual evaluation limits unnecessary iterations where the interface is highly under-resolved.
- PLIC-RDF method: The iterative algorithm calculates RDF gradients, updates interface normals, and stops when normal or curvature-related residuals meet their tolerances.
26 end
The plic-RDF method initializes interface orientations using volume-fraction gradients or previous-step normals, then iterates RDF reconstruction until normal changes and local under-resolution meet stopping criteria.
- Initial estimate: Previous-step interface normals are interpolated from point-neighbour interface cells using weights based on interface area and backward-advected distance.The interpolation uses the cell velocity and the distance between the neighbouring interface location and the backward-advected cell centre.
- Initial estimate: A fallback to the volume-fraction gradient is used when the surrounding interface-normal variation exceeds 10 degrees.This addresses noisy interpolation in highly under-resolved regions, where volume-fraction gradients are considered more reasonable on coarse meshes.
- Stopping criteria: The first plic-RDF residual measures the average difference between interface normals in successive iterations.Its default tolerance is tol = 10^-6.
- Stopping criteria: The curvature-based stopping criterion prevents excessive iterations where coarse volume-fraction data cannot provide enough orientation information.This limits wasted computational resources in under-resolved regions.
4. Numerical tests
The numerical tests evaluate reconstruction and advection implementations using OpenFOAM-based code and structured or unstructured meshes generated with established mesh tools.
- Implementation: The implementation is based on OpenFOAM-v1706, compiled with g++ 4.8.4 using the default -O3 optimization flag.Performance measurements use a dual Intel Xeon 2687W v2 system with 64 GB of DDR3-1600 MHz memory.
- Mesh generation: Triangle-prism meshes are generated with gmsh, polygonal prisms with OpenFOAM’s polyDualMesh, tetrahedra with gmsh and Netgen optimization, and polyhedra with cfMesh.All two-dimensional meshes have thickness 1.
4.1. Numerical results for interface reconstruction
Reconstruction tests compare iso-Alpha, iso-RDF, and plic-RDF on circles and spheres across structured and unstructured meshes using strict maximum-error norms. The RDF-based methods consistently achieve second-order convergence across the tested two- and three-dimensional mesh types, while iso-Alpha converges less reliably on unstructured meshes and is more computationally expensive than Youngs’ algorithm.
- Test setup: The tests reconstruct a circle in 2D and a sphere in 3D on multiple mesh types, repeating each resolution 100 times with randomly shifted interface centres.The largest error from the repetitions is reported to reduce bias from mesh-interface positioning.
- Error measures: The L∞ norm selects the maximum error in the domain, making it a strict test of interface position and normal convergence.Errors are evaluated over interface cells against exact interface positions and unit normals.
- 2D reconstruction: On triangular-prism meshes, iso-Alpha does not converge in orientation or position, whereas iso-RDF and plic-RDF show second-order convergence.The same second-order RDF behaviour is reported for the maximum position error on these unstructured meshes.
- 2D reconstruction: On square-prism meshes, iso-RDF and plic-RDF have second-order convergence for orientation and position, while iso-Alpha converges more slowly for orientation.All three methods are second-order for position on square-prism meshes.
- 2D reconstruction: On polygonal-prism meshes, iso-Alpha lacks convergence in the maximum error, while iso-RDF and plic-RDF retain second-order convergence.The RDF-based methods also show second-order convergence in the reported alternative error measure.
- 3D reconstruction: For spheres on tetrahedral, hexahedral, and polyhedral meshes, iso-RDF and plic-RDF achieve second-order convergence for both orientation and position across the full resolution range.iso-Alpha’s second-order behaviour on coarse meshes deteriorates with refinement except on hexahedral meshes.
- Comparison with other methods: In comparisons using the L1 orientation norm, iso-RDF and plic-RDF achieve second-order convergence in 2D and 3D, with absolute errors comparable to or below the cited alternatives.In 3D their absolute errors are significantly lower than LVIRA’s and similar to the Height Function method.
- Computational cost: Relative to Youngs’ algorithm, iso-RDF is 6.7-11.4 times slower and plic-RDF is 2.4-7.3 times slower.iso-Alpha is the fastest method except on polyhedral meshes, where Youngs’ algorithm is fastest; plic-RDF improves on iso-RDF through its normalized curvature residual.
4.2. Numerical results for interface advection
Replacing iso-Alpha with plic-RDF in isoAdvector reduces advection errors and restores second-order convergence at sufficiently small CFL numbers across structured and unstructured meshes. The method also reduces reconstruction cost through improved normal initialization, while highly deformed interfaces remain resolution-sensitive.
- Circle in constant flow: For triangular prism meshes, plic-RDF remains second-order for CFL = 0.1 and 0.025 across the full resolution range, while iso-Alpha drops to first order on finer meshes.At CFL = 0.5 and 1.0, plic-RDF becomes first order at high resolutions but retains lower absolute error than iso-Alpha.
- Circle in constant flow: On square and polygonal prism meshes, plic-RDF preserves second-order convergence at lower CFL values, whereas larger CFL numbers reduce convergence or increase absolute error.For square prisms, second order holds for CFL = 0.025 and 0.1; for polygonal prisms, it holds for CFL = 0.025, 0.1 and 0.5.
- Circle in constant flow: plic-RDF consistently lowers absolute advection errors and achieves second-order convergence for sufficiently small time steps.This behavior is reported across the tested mesh types in the constant-flow advection case.
- Reconstruction cost: Interpolating interface normals from the previous time step reduces the initial residual and yields second-order residual convergence for all tested prism mesh types.The improved initial guess also reduces the number of reconstruction steps, averaging about 2 to 2.5 interface reconstructions.
- Rotation test: On unstructured tetrahedral meshes, the rotation test shows good volume conservation and boundedness with near-second-order convergence at all resolutions.The method achieves slightly higher convergence rates than Jofre et al., although Jofre et al. has lower coarse-mesh shape errors.
- Disc in reversed spiral flow: At CFL = 1, isoAdvector-plicRDF has larger fine-mesh errors than several comparison methods, while halving CFL to 0.5 restores its convergence.At CFL = 1, it is nevertheless reported as 10 times faster than UFVFC-Swartz; the comparison uses the stated test conditions.
- Disc in reversed spiral flow: The reversed-vortex test is limited by unresolved thin tails: visible pinched-off droplets at maximum deformation do not occur in the exact, well-resolved solution.The limitation arises at the resolutions used in this test and elsewhere in the cited literature.
5. Conclusion
The paper introduces iterative RDF-based interface reconstruction for general meshes and combines it with isoAdvector advection. The method achieves second-order convergence across structured and unstructured meshes, while retaining comparatively low computational cost.
- The iterative RDF procedure estimates local interface position and orientation from raw volume fraction data.
- Second-order convergence is demonstrated for interface position and orientation on 2D and 3D structured and unstructured meshes.
- The reconstruction methods improve convergence on unstructured meshes and for local interface orientation compared with the earlier volume-fraction-isosurface method.
- With isoAdvector advection, second-order convergence is achieved on all tested mesh types when the CFL number is 0.2 or below.
- The implementation is released as an OpenFOAM extension with the paper’s test cases, while computational cost remains roughly an order of magnitude lower than the method of Maric et al.