Source-linked AI summary
A Computational Method for Sharp Interface Advection
Johan Roenby, Henrik Bredmose, Hrvoje Jasak
TL;DR
Interface advection requires estimating fluid transport across mesh faces while preserving a sharp surface. isoAdvector combines isosurface-based reconstruction with evolving face-intersection geometry for arbitrary polyhedral meshes, and its tests report satisfactory conservation, boundedness, sharpness, and efficiency, with convergence orders from 1.7 to 3.2. The method’s reconstruction sacrifices global isoface continuity when different isovalues are used across neighboring surface cells, while optional clipping can sacrifice strict volume conservation.
Problem
The paper addresses how to accurately and efficiently estimate fluid-A volume transport across mesh faces from cell volume fractions and velocity data during a time step.
Method
isoAdvector reconstructs interfaces with volume-matching isofaces and integrates the evolving submerged face area obtained from the moving face-interface intersection line.
Results
The tests on 2D and 3D structured and unstructured meshes report satisfactory shape preservation, volume conservation, boundedness, interface sharpness, and efficiency, with convergence orders from 1.7 to 3.2.
Takeaways & Limitations
isoAdvector provides a geometric VOF interface-advection algorithm designed for general structured and unstructured meshes.
Takeaways & Limitations
Using different isovalues in neighboring surface cells makes the union of isofaces discontinuous, while optional clipping ensures strict boundedness at the cost of strict volume conservation.
Abstract
from arXiv · showhide
We devise a numerical method for passive advection of a surface, such as the interface between two incompressible fluids, across a computational mesh. The method is called isoAdvector, and is developed for general meshes consisting of arbitrary polyhedral cells. The algorithm is based on the volume of fluid (VOF) idea of calculating the volume of one of the fluids transported across the mesh faces during a time step. The novelty of the isoAdvector concept consists in two parts: First, we exploit an isosurface concept for modelling the interface inside cells in a geometric surface reconstruction step. Second, from the reconstructed surface, we model the motion of the face-interface intersection line for a general polygonal face to obtain the time evolution within a time step of the submerged face area. Integrating this submerged area over the time step leads to an accurate estimate for the total volume of fluid transported across the face. The method was tested on simple 2D and 3D interface advection problems both on structured and unstructured meshes. The results are very satisfactory both in terms of volume conservation, boundedness, surface sharpness, and efficiency. The isoAdvector method was implemented as an OpenFOAM(R) extension and is published as open source.
1. Introduction
The paper frames interface advection as estimating fluid transport across mesh faces while preserving a sharp interface. It develops a VOF formulation for this task and targets accurate, efficient operation on arbitrary meshes.
- VOF formulation: VOF methods represent the interface implicitly through cell volume fractions and redistribute one fluid between adjacent cells across mesh faces.The paper distinguishes geometric schemes, which reconstruct the interface explicitly, from algebraic schemes, which do not.
- Design goal: The isoAdvector ambition is a VOF method for arbitrary meshes that retains geometric accuracy while minimizing geometric operations and calculation time.The stated target includes automatically generated unstructured meshes with complex geometries.
- Problem setting: Interface advection advances a surface moving in a prescribed velocity field, often representing the boundary between two fluids.The paper focuses on passive advection with the velocity field assumed known in advance.
- VOF formulation: The paper derives the update from the control-volume continuity equation for a discontinuous density field and treats face transport as the fundamental quantity to estimate.The formulation defines an indicator field for fluid A and integrates transport over the time step.
- VOF formulation: The central numerical question is how to use volume fractions and velocity data to estimate fluid-A transport across each face during a time step.This transported volume determines the updated cell volume fractions and therefore the implicitly represented interface.
2. The isoAdvector concept
isoAdvector reconstructs an interface inside each surface cell with an isoface, then models its motion and the resulting submerged face area over a time step. Face transport is integrated from this evolving geometry and shared consistently between neighboring cells.
- Interface reconstruction: The reconstruction represents the interface with an isoface formed from edge intersections, splitting a polyhedral cell into fluid-A and fluid-B subcells.The isoface is chosen to reproduce the cell’s fluid volume fraction rather than simply using a fixed visualisation isovalue.
- Advection step: For constant velocity in space and time, the face-transport formulation is exact; its approximation error becomes negligible when cells and time steps resolve velocity variations.This conclusion follows from the stated approximation of interface motion within each interval.
- Advection step: The method estimates the time evolution of the face area submerged in fluid A and integrates that area to obtain transported fluid volume.The submerged-area evolution can be discontinuous and non-differentiable because it depends on interface orientation, motion, and face shape.
- Advection step: The advection step tracks the face-interface intersection line as the reconstructed isoface moves from the upwind cell across a polygonal face.The isoface motion is estimated from velocity data in surrounding cells during the time step.
- Algorithm overview: Each face transport is calculated once and reused for both cells sharing an internal face, guaranteeing local and global conservation of both fluids.The algorithm then updates volume fractions, redistributes fluxes when values leave [0,1], and optionally clips them for strict boundedness.
3. Implementation details
The implementation controls time-step size through a target interface Courant number and separates initialization and surface-cell detection from the detailed face-transport calculation.
- Time-step control: A target interface Courant number sets the time-step size so that the chosen limit is not exceeded in surface cells.The interface Courant number concerns only interface-normal velocity in surface cells.
- Algorithm setup: The implementation begins by initializing face transport with upwind values and identifying surface cells using a user-specified volume-fraction tolerance.The detailed geometric transport calculation starts after these two preliminary steps.
3.1 Calculating the initial isoface in a surface cell
The initial isoface is constructed by interpolating volume fractions to mesh vertices, cutting cell edges, and selecting an isovalue that reproduces the cell’s fluid volume fraction.
- Vertex interpolation and edge cutting: Volume fractions are interpolated from cells to mesh points, after which each cell edge is cut where its endpoint values straddle the selected isovalue.The implementation uses inverse point-to-cell-centre interpolation, though other interpolation choices are possible.
- Isoface construction: The edge-cutting points are connected across faces and then into an isoface that divides the cell into polyhedral fluid-A and fluid-B regions.The resulting fluid-A region is used to calculate its volume relative to the cell volume.
- Volume matching: The selected isovalue f* satisfies the reconstructed fluid-A volume fraction ˜α(f*) = αi, with ˜α varying continuously and monotonically from 0 to 1.Between neighboring vertex values, the implementation models this variation with a cubic polynomial and solves for f* using Newton’s method.
- Stencil: Cell-to-vertex interpolation gives the isoface an effective stencil containing the surface cell and all cells sharing one of its vertices.This extends the local information used by reconstruction beyond the cell’s own volume fraction.
3.2 Estimating the isoface motion during a time step
The method computes isoface geometry and interpolates a normal propagation speed to approximate interface motion during each time step, while noting extensions for acceleration and rotation.
- The isoface centre xS is computed as an area-weighted average of triangle centres formed from the polygon’s average point.
- Velocity data are interpolated to xS by decomposing the cell into tetrahedra and linearly interpolating within the containing tetrahedron.
- The scalar normal speed US = US · n̂S determines isoface propagation, with positive US filling fluid A and negative US filling fluid B.
- The implementation treats US as constant throughout the time step.
- Future work could estimate isoface acceleration and rotation using prior velocities and surrounding-cell velocity gradients.
3.3 Evolution of the face-interface intersection line
The face-interface intersection line is advanced by estimating when a moving isoface reaches each face vertex, then interpolating the line between successive vertex-hit times.
- For each downwind face, vertex-hit times are estimated from the isoface centre, normal, and scalar normal speed.
- The vertex-hit estimate is tk ≈ t + (Xk − xS) · n̂S/US for each face vertex k.
- At an intermediate time τ, a linear edge-cutting interpolation reconstructs the face-interface intersection line using vertex times instead of volume fractions.
- Between successive vertex-hit times tk and tk+1, the line endpoints are interpolated so the line sweeps a quadrilateral.
- This procedure approximates the intersection-line evolution over the complete time step.
3.4 Time integral of submerged face area
The submerged face-area integral is evaluated by splitting the time step at vertex-hit events, integrating each quadrilateral-sweep interval, and combining the resulting contributions to update volume fractions.
- A sorted time list begins at t, ends at t + ∆t, and includes every vertex-hit time lying inside the time step.
- The total submerged-area integral is split across the resulting consecutive subintervals.
- On each subinterval, the intersection line sweeps a quadrilateral and defines the intermediate submerged area Aj(τ).
- Analytically calculated polynomial coefficients Pk and Qk determine each subinterval’s time-integral contribution.
- Summing all subinterval contributions yields the estimate for transported volume ∆Vj(t, ∆t).
- The updated volume fractions αi at t + ∆t are calculated by inserting the face transport volumes ∆Vj into the update equation.
3.5 Bounding procedure
The bounding procedure redistributes small volume-fraction overshoots and undershoots through downwind face fluxes, while optional clipping enforces strict bounds at the cost of volume conservation.
- The transport estimate can produce updated volume fractions outside the physical range 0 ≤ αi(t + ∆t) ≤ 1, especially near completely emptied or filled cells.
- For overfilled cells, surplus fluid A is passed to downwind cells because upwind neighbours are typically already filled with fluid A.
- Redistribution begins when an overfilled cell becomes full, after which pure fluid A flows through its downwind faces.
- Surplus V+ is distributed among downwind faces using flux weights φj, with each face limited to its transport capacity φj∆t.
- If capacities prevent complete redistribution, the process repeats over remaining faces until the surplus is passed onward or no capacity remains.
- Lower-bound correction transforms the problem using βi ≡ 1 − αi and the corresponding fluid-B transport volumes.
- Optional clipping enforces strict boundedness but makes the method not strictly volume conserving, so total fluid-A volume must be monitored.
4. Results
The tests evaluate shape preservation, volume conservation, boundedness, sharpness, convergence, and efficiency across structured and unstructured 2D and 3D advection cases. isoAdvector generally outperforms the reference schemes, although severe interface stretching can introduce pinching-related distortion.
- Shape preservation: isoAdvector preserved disk shape best on square meshes, with E1 error at least 7 times smaller than the best reference scheme in every run.MULES and HRIC tended toward mesh-aligned shapes, while CICSAM approached a circular solution only as the time step decreased.
- Shape preservation: ∼2.4 was isoAdvector’s convergence order under mesh refinement for the square-mesh Co = 0.5 runs.Its error improved only slightly from Co = 0.5 to 0.2 and became slightly worse from Co = 0.2 to 0.1 on a fixed mesh.
- Conservation and boundedness: Machine-precision volume conservation and bounded volume fractions were achieved by isoAdvector across the tested square-mesh cases.MULES and HRIC also remained bounded, while CICSAM had severe bounding problems; isoAdvector preserved volume even on the coarsest mesh.
- Sharpness and efficiency: isoAdvector produced interface thickness close to the exact volume-fraction representation and was slightly faster than CICSAM and HRIC while 2–4 times faster than MULES.CICSAM was the sharpest reference scheme, whereas HRIC showed substantial interface smearing.
- Unstructured meshes: On unstructured 2D meshes, isoAdvector performed significantly better than the reference schemes, with calculation times similar to HRIC and CICSAM and lower than MULES.Reference-scheme behavior varied strongly by mesh: HRIC remained diffusive, while CICSAM developed threads and unboundedness on triangular meshes.
- Spiralling disc: ∼10 times slower was MULES than isoAdvector in the spiralling-disc comparison, while its E1 error was ∼50% larger on square meshes and 5 times larger on polygon meshes.At Co = 0.5, MULES results on all three meshes became completely disintegrated.
- Three-dimensional advection: 2.6–3.2 was the calculated convergence-order range for the 3D sphere advection tests on random tetrahedral meshes.MULES had significantly worse shape preservation and was 20 or 5 times slower than isoAdvector in the Co = 0.1 and Co = 0.5 comparisons, respectively.
5. Conclusion
isoAdvector is a geometric VOF algorithm for interface advection on general structured and unstructured meshes, with novel reconstruction and advection steps. Tests found satisfactory accuracy and efficiency, while the explicit method has Courant-number constraints.
- isoAdvector derives numerical interface advection from the control-volume integrated continuity equation for a discontinuous density field.
- Isosurface calculations estimate fluid distributions inside cells and avoid gradient calculations that may cause problems on unstructured meshes.
- Time-step subdivision enables analytical face-flux calculation under steady interface motion, without assuming a particular polygonal face shape.
- Mesh-refinement tests produced convergence orders between 1.7 and 3.2, alongside satisfactory shape preservation, volume conservation, boundedness, sharpness, and efficiency.
- The explicit advection step is principally limited to Courant numbers in [0, 1], with optimum accuracy around Co ∼0.5 and only small degradation near 1.
- The code is published as an open-source OpenFOAM extension, and the method may also apply to passive tracer fields beyond fluid interfaces.