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 of 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 two-dimensional and three-dimensional interface advection problems on both structured and unstructured meshes. The results are very satisfactory in terms of volume conservation, boundedness, surface sharpness and efficiency. The isoAdvector method was implemented as an OpenFOAM® extension and is published as open source.
In this paper, we address the numerical challenge of advancing a surface moving in a prescribed velocity field. We will refer to this as the interface advection problem, because the surface often constitutes an interface e.g. between two fluids. Simple as the problem may appear, there is a large variety of problems in science, engineering and industry where its solutions are far from trivial. Our motivation for addressing this problem is rooted in our usage of computational fluid dynamics (CFD) as a practical engineering tool for calculating wave loads on coastal and marine structures. Whether it is an offshore wind turbine foundation, or an oil and gas platform, accurate estimation of the peak loads from violent breaking waves is paramount for the correct dimensioning of the structure. In our view, CFD has a large unexploited potential to improve wave load estimates, and to reduce both cost and risks in the design phase of coastal and offshore structures.
Owing to the omnipresence of interfacial flows, the list of areas that may benefit from improved solution methods to the interface advection problem is almost endless. Some examples are bubble column reactors, oil–gas mixtures in pipelines, inkjet printing, automotive aquaplaning, ship manoeuvring, tank sloshing, dam breaks, metal casting processes and hydraulic jumps.
During the past 40–50 years, both Lagrangian and Eulerian strategies have been employed to develop a wide range of methods for advecting a sharp interface . Today, most CFD codes for practical engineering calculations use variants of the volume-of-fluid (VOF) method for the interface advection step in their interfacial flow solvers. This includes current versions of ANSYS Fluent®, STAR-CCM+®, Gerris , OpenFOAM® [3,4] and many others. In the VOF methodology, the interface is implicitly represented via the volume fractions of one of the fluids in computational cells. The advection is done by redistributing the content of this fluid between adjacent cells by moving it across the mesh faces. Since the first VOF methods appeared in the literature , a large variety of VOF schemes have been developed. They may be divided into two categories: geometric methods involving an explicit reconstruction of the interface from the volume fraction data, and algebraic methods making no such attempt. Algebraic VOF schemes are typically much simpler to implement, more efficient and are not restricted to structured meshes. They are, however, founded on much more heuristic considerations and are not as accurate as the geometric VOF schemes . Geometric VOF schemes, on the other hand, involve complex geometric operations making their implementation cumbersome and their execution slow. Geometric VOF methods for unstructured meshes is an active area of research [7–13].
Our ambition in the development of the isoAdvector algorithm is to develop a VOF-based interface advection method that works on arbitrary meshes, retains the accuracy of the geometric schemes by explicitly approximating the interface, and yet keeps the geometric operations at a minimum in order to obtain acceptable calculation times. An efficient VOF scheme yielding accurate results even on automatically generated unstructured meshes of complex geometries has a huge potential for speeding up the simulation process and making CFD an integrated part of design processes involving interfacial flows.
In the remainder of this section, we give an introduction to the interface advection problem and its formulation in the VOF framework. In §2, we present the new ideas and concepts of the isoAdvector method, and give an overview of the steps involved in the numerical procedure. The implementation details and considerations involved in each step are described at length in §3. In §4, we demonstrate the performance of the new method with a series of simple test cases. Finally, in §5 we summarize our findings.
1.1. VOF formulation of the interface advection problem
We consider a computational domain in which a surface is embedded. The surface may consist of any number of closed surfaces and may also extend to the boundaries of the domain. We will think of the surface, , as the interface between two incompressible, immiscible fluids denoted by A and B, and occupying the two closed regions, and , satisfying and .
The fluid particles are assumed to be passively advected in a continuous, solenoidal velocity field, u(x,t), which is defined in the whole domain, . In practical engineering applications involving incompressible two-phase flows, the time evolution of the velocity field is governed by the Navier–Stokes equations for u coupled with a Poisson equation for the pressure, p. This system of equations must be solved simultaneously with the interface advection problem. In this work, we focus entirely on the interface advection problem, thus assuming u(x,t) to be known in advance for all points, , and all times, t.
We will now represent the surface in terms of a density field, ρ(x,t), which takes one constant value, ρA, everywhere in and another constant value, ρB, everywhere in . The density field thus has a discontinuity at the interface .1 The evolution of the surface is then determined by the integral form of the continuity equation, 1.1where is an arbitrary stationary volume, is its boundary and dS is the differential area vector pointing out of the volume. In words, this mass conservation equation says that the instantaneous rate of change of the total mass enclosed in the volume is given by the instantaneous flux of mass through its boundary.
In the pure advection problem with a predetermined velocity field, the specific values of the fluid densities, ρA and ρB, are immaterial, that is, the solution does not depend on them. To remove these insignificant parameters from the problem, we define the indicator field, 1.2such that H=1 for all , and H=0 for all .
We now discretize the computational domain, , by conceptually dividing it into a large number of control volumes, or cells, , for i=1,…,NC. If two cells i and j are adjacent, their shared boundary, , is called an internal face. If cell i touches the domain boundary, the shared surface, , will consist of one or more boundary faces. All faces are labelled with integers, j=1,…,NF, and the surface of face j is denoted by . Thus, the boundary of the cell i may be represented by a list, Bi, of all the labels of faces belonging to its boundary .
With these mesh definitions in place, we can now substitute (1.2) into (1.1) with cell i as the volume of integration, 1.3Because face j has its own orientation determining the direction of dS, we have introduced the auxiliary factor sij=+1 or −1, such that sij dS points out of cell i for face j.
The natural next step is to define the volume fraction of fluid A in cell i, 1.4where V i is the volume of cell i. Substituting (1.4) into (1.3), and formally integrating (1.3) from time t to time t+Δt, we obtain the following equation for the updated volume fraction of cell i, 1.5We stress that this equation is exact with no numerical approximations introduced yet. It is the fundamental equation from which one must derive any consistent interface advection method. The time integral on the right-hand side is the total volume of fluid A transported across face j during the time interval from time t to t+Δt. It is the fundamental quantity that we must estimate in order to advance αi, and hence implicitly the surface , in time. We will call this quantity ΔV j(t,Δt): 1.6The fundamental equation (1.5) can then be formulated as 1.7
Before we move on to present the basic ideas of the isoAdvector method, we will need to consider how the velocity field is represented. In the finite volume treatment of the fluid equations of motion, the natural representation of the velocity field is in terms of cell-averaged values, 1.8As the convective terms in the governing fluid equations give the transport of mass, momentum, etc. across cell faces, other important velocity field representations are the volumetric fluxes across mesh faces, 1.9
The question we will try to answer in the following can now be formulated as follows:
How do we most accurately and efficiently exploit the available information at time t, i.e. the volume fractions, αi, and the velocity data, ui and ϕj, to estimate the fluid A volume transport, ΔVj(t,Δt), across a face during the time interval [t,t+Δt]?
2. The isoAdvector concept
We will now present the general ideas behind the isoAdvector method, starting with the interface representation using isosurfaces, then introducing the concept of a face–interface intersection line moving across a face, and finally giving an overview of the steps involved in the numerical procedure. For the sake of clarity, we focus on ideas in this section, and postpone the detailed description of the implementation to §3. For full implementation details, the reader is referred to the source code provided with this article .
2.1. The interface reconstruction step
The integral in (1.6) is highly dependent on the local distribution of fluid A and B inside cell i and inside its neighbour cells from which it receives fluid during the time step. However, the volume fractions hold no information about the distribution of the two fluids inside the cells. We must therefore come up with a subgrid model for this ‘intracellular’ distribution from the given volume fraction data. If the volume fraction data are ‘sharp’, only cells very close to the interface will have volume fractions significantly different from 0 and 1. Then, if cell i is on the interface, its neighbours in one direction will mainly contain fluid A, while its neighbour cells in the opposite direction will mainly contain fluid B. In words, we want our subgrid model to capture this local distribution information, and place the fluid A content of cell i close to the neighbours containing fluid A (which is equivalent to its fluid B content being placed near the neighbours containing fluid B). The implicit assumption made in this model is that the interface is sufficiently well resolved by the mesh such that the local radius of curvature is larger than the cell size. Whenever this is satisfied, an isosurface calculation will provide a good estimate of the required local fluid distribution information.
The idea of using isosurface calculations in the interface reconstruction step is inspired by our work with postprocessing of interfacial CFD data. Here, it is customary to visualize the fluid interface (e.g. using ParaView® ) by showing the 0.5-isosurface based on the volume fraction data. Such numerically calculated isosurfaces are topologically consistent continuous surfaces and are straightforward to calculate on arbitrary polyhedral meshes. The numerical representation of an isosurface in a polyhedral cell is a list of the points where the isosurface cuts the cell edges. See red points in figure 1a for an illustration. This list of points represents a face which cuts the cell into two polyhedral subcells, with one completely immersed in fluid A and the other completely immersed in fluid B. We will call such a face an isoface. See the green patch in figure 1a for an example. We note that if an isoface has more than three vertices, it will generally not be exactly planar.
When calculating an isosurface from the volume fraction data, we have the freedom of choosing an isovalue between 0 and 1. Which isovalue should we choose? For surface visualization from volume fraction data, we usually plot the 0.5-isosurface. This, however, is not a good choice for the surface reconstruction step in an interface advection algorithm, because the isoface in cell i with isovalue 0.5 does not in general cut it into two subcells of the volumetric proportions dictated by the volume fraction, αi. It may, for instance, occur that the cell has αi=0.8, and yet is not even cut by the 0.5-isosurface. In this case, a surface reconstruction model based on the 0.5-isosurface would say nothing about how the 80% fluid A and 20% fluid B is distributed inside this cell. There will, however, exist an isoface with another isovalue, which will cut the cell into subcells of the correct volumetric proportions. An important component of our proposed scheme is an efficient method for finding this isovalue for a given surface cell (see §3 for details). We note that with the use of different isovalues in different surface cells, the union of isofaces is no longer a continuous surface, as it would be if the same isovalue was used in adjacent cells. This is the price we must pay to ensure that isofaces cut surface cells into subcells having the correct volumetric proportions.
2.2. The advection step
Most Navier–Stokes solvers for interface flows use a segregated solution approach, in which the coupled system of equations governing the flow are solved in sequence within a time step. This means that at the point where the interface is to be advected from time t to time t+Δt, we only have information about the velocity field up to time t. But as seen in the integrand in (1.5), calculation of the updated αi requires information about the velocity field on the interval [t,t+Δt]. We must therefore estimate the evolution of the velocity field during the time step. The simplest such estimate is to regard the velocity field as constant (in time) during the whole time step. With this assumption, we can write u(x,τ)≈u(x,t) in (1.6). Another assumption we will make in (1.6) is that u on the face dotted with the differential face normal vector, dS, can be approximated in terms of the volumetric face flux, ϕj (defined in (1.9)), as follows: 2.1where dS≡d|S|, and the face normal is given by 2.2Substituting this into (1.6), we obtain 2.3The remaining surface integral in (2.3) is then simply the instantaneous area of face j submerged in fluid A, which we will call Aj(τ): 2.4Using this definition, we may now write (2.3) as 2.5An important point is that in the special case, where the velocity field is constant in both space and time, (2.5) is exact. Thus, if the cells become sufficiently small compared with the gradients of the velocity field, and the time steps become sufficiently small compared with the temporal variations in the velocity field, the error committed in the above approximation becomes negligible.
As is seen from (2.5), the challenge in constructing a VOF scheme is to estimate the time evolution within a time step of the submerged (in fluid A) area of a face, and then integrate this area in time. The time scale on which Aj(τ) changes is not dictated by the time scales of the flow, but by a complicated combination of the relative orientations of the face and interface, the direction of motion of the interface and the shape of the specific polygonal face. As an example, consider a planar interface approaching a planar face to which it is parallel. In this case, Aj(τ) will be a discontinuous function of τ. This in turn makes ΔV j(t,Δt) non-differentiable with respect to Δt. The discontinuous and non-differentiable nature of these quantities is what makes the interface advection problem so difficult to attack with the traditional weaponry of numerical analysis, which relies on the existence of a Taylor expansion of the sought solution.
In the isoAdvector advection step, when we calculate Aj(τ) for face j, our starting point is the isoface in the cell upwind of face j at time t, because this is the cell from which the face receives fluid during the time step. The motion of this isoface within the time step [t,t+Δt] may be approximated by using the velocity data in the surrounding cells. Figure 1b shows an example of how the isoface may appear at three times during the time step. For details on our approximation of the isoface motion, the reader is referred to §3.2.
Knowing the isoface position and orientation inside cell i at any time within [t,t+Δt], we also know for its downwind face j the face–interface intersection line (see blue lines in figure 1a) at any time during the time interval. With this information, the time integral in (2.5) can be calculated analytically, to finally obtain our estimate of the total volume of fluid A transported across face j during the time interval [t,t+Δt].
We stress that the fluid A transport across a face is only calculated once for each face, and that, for internal faces, this value is used to update the volume fractions of both of the two cells sharing the face. This guarantees local and global conservation of each of the two fluids A and B.
2.3. Algorithm overview
We here give an overview of the steps taken in the isoAdvector algorithm to advance the volume fractions from time t to time t+Δt:
Step 1 For each face j, initialize ΔV j with the upwind cell volume fraction, ΔV j=αupwind(j)ϕjΔt.
Step 2 Find all surface cells, i.e. cells with ϵ<αi(t)<1−ϵ, where ϵ is a user-specified tolerance (we typically use 10−8).
Step 3 For each surface cell i, do the following:
3.1 Find its isoface, i.e. the isosurface inside the cell with isovalue such that it cuts the cell into the correct volumetric fractions, αi(t) and 1−αi(t) (details in §3.1).
3.2 Use the velocity field data to estimate the isoface motion during the time interval [t,t+Δt] (details in §(b)).
3.3 For each downwind face j of surface cell i, use the isoface and its motion to calculate the face–interface intersection line during the time interval [t,t+Δt] (details in §3.1).
3.4 For each downwind face j of surface cell i, use the motion of its face–interface intersection line to calculate ΔV j(t,Δt) from the time integral in (2.5) (details in §3.4).
Step 4 For each cell calculate αi(t+Δt) by inserting the ΔV js of its faces in (1.7).
Step 5 For cells with αi(t+Δt)<0 or αi(t+Δt)>1 adjust the ΔV js of its faces using a redistribution procedure and recalculate αi(t+Δt) by inserting corrected ΔV js in (1.7). This step also includes an optional subsequent clipping of any cell values αi<0 or αi>1 to ensure strict boundedness before proceeding to the next time step (details in §(e)).
3. Implementation details
In this section, we provide the implementation details of the procedure outlined in §2.3. We first note that the time step size may vary between time steps. The user can specify a target interface Courant number, Co, based on which time step size is set at the beginning of each time step, to ensure that Co is not exceeded in any surface cell. The interface Courant number only concerns the velocity of the interface normal to itself in surface cells.
Step 1 in §2.3, where we initialize ΔV j with upwind values, and Step 2, where we find all surface cells with ϵ<αi<1−ϵ, need no further explanation. We will therefore jump to step 3, which contains the actual calculation of the volume transport across faces.
3.1. Calculating the initial isoface in a surface cell
The first step in calculating the isosurface is to interpolate the volume fractions to the grid points of the mesh. The value in a grid point will in general be a linear combination of the volume fractions in the cells sharing this point. Here, we have chosen to use the inverse of the distances between a grid point and the surrounding cell centres as the weights in this interpolation step. Other options could be to use the cell volumes or solid angles as interpolation weights. A systematic study of the influence of this choice is left for future work.
Let us temporarily denote the N vertices of cell i by X1,…,XN and the corresponding interpolated volume fractions by f1,…,fN. The cell edges are straight lines between pairs of points in the vertex list. To construct the f-isoface for cell i, we go through all cell edges and cut them by linear interpolation of the edge vertex values: if the edge (Xk,Xl) has values fk<f and f<fl, the edge is cut at the point 3.1Once all such edge cutting points have been found for cell i, they can be connected across faces to form the face–interface intersection lines, which can again be connected to form the isoface inside the cell (figure 1a). The isoface splits cell i into a polyhedral cell, , entirely in fluid A, and another cell, , entirely in fluid B. We can calculate the volume of relative to the cell volume, 3.2This will vary monotonically and continuously from 0 to 1, as the isovalue f varies from the maximum vertex value, max( f1,…,fN), to the minimum vertex value, min( f1,…,fN). As argued in §2.1, the correct isovalue to use is the one recovering the cell volume fraction. That is, we should find f* such that . In the current implementation, f* is found by the following procedure. First, we geometrically calculate for the vertex values, f1,…,fN. Of these, we then find the closest value below and above f*, say, fk and fl, such that f*∈[fk,fl]. On this interval, we know that varies monotonically like a cubic polynomial. Thus, evaluating geometrically at two points between fk and fl, we have four equations for the four polynomial coefficients. The resulting linear 4×4 Vandermonde matrix system we solve using LU decomposition. With a polynomial expression for at hand, we can use Newton’s root finding method to efficiently find f* such that , where ϵ is a user-specified tolerance, typically set to ϵ=10−8. In rare cases, the LU solution does not give useful coefficients because the 4×4 matrix is ill-conditioned, so the method does not converge. In these cases, we use Newton’s root finding method with direct geometric evaluation of instead of the much cheaper polynomial evaluation.
We note that, due to the cell-centre-to-vertex interpolation of the volume fractions, the effective stencil contributing to the isoface inside a surface cell consists of the cell itself with all its point neighbours, that is, all surrounding cells with which it shares a vertex.
3.2. Estimating the isoface motion during a time step
We first calculate the geometric face centre, xS, and the unit normal vector, , of the isoface (figure 1a). The procedure for doing this is the same as for any other mesh face in OpenFOAM®. The average point between the N vertex points of the N-gonal face is calculated. Then the face is decomposed into N triangles, all sharing this average point as their common apex. The face centre, xS, is then calculated as the area-weighted average of the geometric centres of these N triangles. Likewise, the face normal vector, nS, is calculated as the area-weighted average of the N triangle area vectors.
The next step is to interpolate the velocity data, ui(t), to xS. This is done by first decomposing the cell into tetrahedra, all sharing the cell centre as their common apex. Then we find the tetrahedron containing xS and interpolate the velocity field into its vertices. Finally, we interpolate linearly from the tetrahedral vertices to obtain the velocity vector US at xS. We note that, for stationary meshes, the weightings in this interpolation procedure only need to be calculated and stored once at the beginning of a simulation.
The next step is to dot US with the isoface normal, , to obtain the isoface motion normal to itself, . We will make the convention that is directed from fluid A into fluid B. Thus, positive US means that the cell is filling up with fluid A, while negative US means that it is filling up with fluid B. In the current implementation, we regard the US of an isoface as constant during the whole time step. Possible improvements could be (i) using velocity data from previous time steps to estimate the isoface acceleration during the time step and (ii) calculating the velocity gradient from surrounding cell velocity data to approximate the isoface rotation around its two tangential axes during the time step. Work along these lines is left for future development.
3.3. Evolution of the face–interface intersection line
We now use xS, and US to approximate the time evolution of the face–interface intersection line of a face j, which is downwind of surface cell i. This we do by estimating the times at which the isoface, travelling with velocity US normal to itself, will reach the vertex points of face j (figure 2a). Let us temporarily denote the N vertex points of face j by X1,…,XN, and the times at which the isoface passes these points by t1,…,tN. Then we can estimate these times as 3.3To obtain the face–interface intersection line at a given time τ∈[t,t+Δt], we can now apply a linear interpolation-based edge-cutting procedure equivalent to the one used to find the initial (i.e. at time t) isoface from the volume fractions. Only now the function values at the vertices are the times from (3.3), rather than the interpolated volume fractions.
More specifically, let us temporarily denote by AB the line segment at time tk, and by CD the line segment at time tk+1, such that ABCD is the grey quadrilateral shown in figure 2a,b. Then, at an intermediate time, τ∈[tk,tk+1], we will assume the two endpoints of the face–interface intersection line segment to be 3.4as illustrated in figure 2b. This concludes our approximation of the face–interface intersection line evolution during a time step.
3.4. Time integral of submerged face area
To calculate the time integral of the submerged area, Aj(τ) in (2.5), we first generate a sorted list of times, , starting with , and ending with , and with a sorted list of all the tks from (3.3) satisfying t<tk<t+Δt in between. Then, the time integral in (2.5) can be split up as follows: 3.5On each of these subintervals, the face–interface intersection line sweeps a quadrilateral as the one shown in figure 2b. Using the definition in (3.4), the submerged area at the intermediate time is 3.6Here, Pk and Qk are polynomial coefficients that can be calculated analytically from and . The sign of US in cell i accounts for the direction of propagation of the isoface, i.e. whether the cell and face are gaining or losing fluid A during the time interval. Once these coefficients are obtained, the contribution to the time integral in (3.5) from the sub-time interval is simply 3.7Adding up all these sub-interval contributions, as devised by (3.5), and substituting the result into (2.5), we finally reach the sought estimate for ΔV j(t,Δt).
As stated in §2.3, the above procedure should be repeated for all downwind faces of a surface cell. On all other faces, we set ΔV j equal to the volume fraction in their upwind cell multiplied by ϕjΔt. The updated αis at time t+Δt can now be calculated by inserting the ΔV js into (1.7).
3.5. Bounding procedure
The procedure described above gives an accurate estimate of the fluid transport across faces in many simple cases. It does not, however, guarantee strict boundedness. That is, there is nothing preventing the algorithm from producing updated volume fractions outside the physically meaningful range 0≤αi(t+Δt)≤1. Experience shows that slight unboundedness may be produced in cells just behind (i.e. upwind of) the interface. The explanation for this is that while the method’s estimate of the ΔV js is typically very good, it is not exact. Thus, in cases where a cell is completely emptied or filled during the time step, the small error committed causes the algorithm to miss 0 or 1 by a small amount. If the produced over- and undershoots are sufficiently small, one might be tempted to simply introduce a step in the algorithm that chops αi(t+Δt) at 0 and 1 before proceeding to the next time step. However, because this corresponds to removing and adding fluid in cells, this method destroys strict volume conservation and is not true to the VOF idea of only allowing redistribution of fluid among cells. While such a step may be practically necessary in order to ensure strict boundedness, it should be used with caution as it may potentially cause severe lack of volume conservation, in particular for long-duration simulations. Can we instead introduce a bounding procedure, which is not adding or removing fluid from the domain, but only redistributes it in order to achieve boundedness? In the following, we will first explain our upper bounding procedure for redistributing the surplus of fluid A in cells with αi(t+Δt)>1. Then, we show how the exact same procedure can be used for lower bounding.
3.5.1. Upper bounding
Cells with αi(t+Δt)>1 are typically just upwind of the interface, in regions where the interface is moving into fluid B (i.e. US>0). Therefore, the cells just upwind of an overfilled cell i are filled with fluid A, and are therefore not good candidates for taking over the surplus of fluid A in cell i in a redistribution step. On the other hand, the cells just downwind of cell i are only partially filled with fluid A, and are therefore able to receive cell i’s small surplus of fluid A. But if cell i has more than one downwind cell, how should its surplus of fluid A be distributed among these? We argue as follows: the overshooting of cell i starts at the time t*∈[t,t+Δt], where the cell becomes filled, i.e. αi(t*)=1. From this time on, all its faces must be completely filled with fluid A. Therefore, pure fluid A will flow through its downwind faces from time t* and onwards. It is therefore natural to pass cell i’s surplus of fluid A through its downwind faces using the face fluxes, ϕj, as the weighting factors. So if the fluid A surplus in cell i is V +, and the cell has N downwind faces with fluxes ϕ1,…,ϕN, then the fraction of V + passed on through the j’th of these faces should be . However, we will not permit more fluid A to be passed through face j than ϕj(t)Δt. Therefore, we will clip the extra flux through face j to . If a face reaches its maximum fluid A transport capacity, so the surplus flux is clipped in this way, the result is that not all the surplus V + in cell i is passed on to downwind cells in this first redistribution step. In that case, the step is repeated to pass on the remaining surplus of fluid A through the remaining downwind faces, still using the ϕjs as weightings, and clipping if the maximum capacity of a face is reached. The step is repeated until either all surplus fluid A in cell i is passed on to the downwind neighbours, or there are no more downwind cells that can take up more fluid A.
3.5.2. Lower bounding
The procedure for lower bounding (i.e. correcting cells with αi(t+Δt)<0) follows simply by changing our perspective from that of fluid A to that of fluid B. We introduce the volume fraction of fluid B, βi≡1−αi, and the volume of fluid B transported across faces during Δt, . Now αi<0 is equivalent to βi>1, and we can apply the upper bounding procedure outlined above to correct the s. With the s corrected, we calculate and insert in (1.7) to obtain the updated volume fraction αi(t+Δt).
It is our experience that the redistribution process outlined above succeeds in bounding most cells. However, occasionally all downwind faces of an overfilled cell will reach their maximum flux capacity before the cell is fully bounded. This only happens on rare occasions, and when it does it only has a minor effect on the overall quality of the solution. Nevertheless, some applications may require strict boundedness at all times, and so we have to introduce an optional clipping of the volume fractions after the bounding procedure described above and before proceeding to the next time step. When this clipping is switched on the method is not strictly volume-conserving, and one should therefore monitor the evolution of the total volume of fluid A, to ensure that it only varies within acceptable limits.
In the following, we present the results of simple test cases with isoAdvector. The numerically advected volume fractions should reproduce the solution to an interface advection problem as accurately as possible with the given mesh and time step size. A simple check of this is to advect a confined volume of fluid A across the computational mesh in a uniform velocity field, and observe to what extent the method preserves the shape of the volume as it should.
The other type of test we will perform exploits the time reversibility of the advection problem. If we advect a confined volume of fluid A in a spatially and temporally varying velocity field for a period of time, the interface will be distorted. If we then reverse the flow, and run it backwards for the same amount of time, the volume should return to its initial position and shape.
The following error measures will be used to quantify the solution quality:
— Shape preservation. Our quantitative measure of shape preservation will be 4.1where the sums are over all cells and is the volume fraction representation of the exact solution at time t.
— Volume conservation. The relative change in the total volume of fluid A in the domain relative to the initial fluid A volume, 4.2should be zero in simulations, where no fluid A enters or leaves the domain.
— Boundedness. For the volume fractions to be physically meaningful, we should have 0≤αi≤1 for i=1,…,NC. Our measures of unboundedness will be mini(αi) and maxi(αi), where the minimum and maximum are taken over all cells at the end of a simulation.
— Sharpness. For a sharp interface, the width of the region where αi changes from 0 to 1 should be similar to the cell size. As the quantitative sharpness measure, we use the volume between the α=0.01 and 0.99 isosurfaces of the volume fraction data divided by the corresponding volume for the volume fraction representation of the exact solution. We will call this quantity δWrel.
— Efficiency. Here, we give the simulation times, Tcalc, in seconds. All simulations were executed on a single core of an Intel Xeon 3.10 GHz CPU (E5-2687W) on a Dell Precision T7600 Workstation.
For benchmarking the isoAdvector algorithm, we compare its performance with three algebraic VOF schemes, which were all developed for arbitrarily unstructured meshes:
— MULES. The interface compression scheme implemented in the OpenFOAM® interfacial flow solver, interFoam. For a good description of this scheme, see Deshpande et al. . The scheme does not have a name, but because it uses the Multidimensional Universal Limiter with Explicit Solution (MULES) to keep the VOF data bounded, we will refer to it simply as MULES. All MULES calculations presented in the following were executed using the interFoam solver in OpenFOAM-2.2.0 with the velocity–pressure coupling calculation switched off. As our CFD work is mainly based on OpenFOAM®, and our primary aim is to improve its interFoam solver, the main emphasis will be on benchmarking against MULES in the subsequent test cases.
— HRIC. The High Resolution Interface Capturing scheme , which is, for instance, used in the commercial CFD software STAR-CCM+®.
— CICSAM. The Compressive Interface Capturing Scheme for Arbitrary Meshes , which is, for instance, one among several available options in ANSYS Fluent®. In all subsequent CICSAM simulations, we use the recommended blending factor value 0.5. It should be noted that Fluent also has a ‘Geo-Reconstruct’ option, which is a geometric VOF method. We have been unable to find a detailed description of this method in the literature.
For the HRIC and CICSAM calculations, we use our own implementations of the schemes in OpenFOAM®. The schemes are available together with the isoAdvector code in the repository , where all set-up files for the following test cases may also be found. To also benchmark isoAdvector against geometric VOF schemes, we include for the test case in §4.4 a comparison with error measures listed in the literature for that case.
4.1. Disc in steady uniform two-dimensional flow
We start by considering a very simple two-dimensional case on a mesh consisting of square cells: a circular region of fluid A of radius R=0.25 moving in a constant and uniform velocity field, u=(1,0.5). The initial volume fractions are obtained from the R-isosurface of the function , where (x0,y0)=(0.5,0.5) is the initial position of the disc centre. Figure 3 shows the volume fraction representations of the exact initial and final interface in grey scale, with white and black cells meaning empty and filled with fluid A, respectively. The α=0.5 contour is shown in blue, and the 0.01 and 0.99 contours are shown in green to indicate the minimal interface width on the given mesh resolution. In the top left corner of figure 3, we also show a zoom on the initial configuration with the exact circle shown in red.
4.1.1. Square meshes
In figure 4, we show in four columns (left to right) the final volume fraction solutions obtained with isoAdvector, MULES, HRIC and CICSAM with five combinations of mesh and time resolution. In rows 1–3, we investigate the effect of refining the mesh resolution with fixed Courant number, Co=0.5. Then, in row 3–5, we use the finest mesh and reduce Co from 0.5 to 0.2 and 0.1. Error measures and calculation times are displayed in table 1. From figure 4 and table 1, the following observations can be made.
Shape preservation. The visual impression from figure 4 is that isoAdvector is superior at preserving the shape of the disc on all shown mesh-Courant number combinations. MULES has a tendency to align the interface at 45 degrees with the mesh faces. Therefore, the MULES solution converges to a tilted square shape as cell and time step sizes are refined (2nd column in figure 4). The HRIC scheme shows a tendency to align the interface with the mesh faces, as also reported in Nielsen . This causes the initially circular interface to converge to a square (3rd column in figure 4). For all the Co=0.5 runs (4th column, rows 1–3 in figure 4), CICSAM does not perform very well in terms of shape preservation. However, it is the only one of the reference schemes which converges to something resembling a circular interface solution as the time step is decreased (lower right corner in figure 4). Table 1a quantifies these observations, showing that the isoAdvector E1 error is at least a factor of 7 smaller than the best of the other schemes for all runs. The table also reveals that the isoAdvector solution only improves slightly, when going from Co=0.5 to Co=0.2, and becomes slightly worse from Co=0.2 to 0.1. Increasing errors with decreasing time step size on a fixed mesh was also reported in Ubbink & Issa . From the three Co=0.5 errors in table 1a, we calculate isoAdvector’s order of convergence with mesh refinement to be ∼2.4.
Volume conservation. From table 1b, we see that isoAdvector is the only scheme with volume preservation down to machine precision even on the coarsest mesh. On the finest mesh MULES also performs very well, followed by HRIC. CICSAM is the worst performing scheme in this comparison.
Boundedness. From tables 1c and 1d, we see that isoAdvector keeps the volume fraction data bounded to within machine precision. Also MULES and HRIC produce bounded volume fractions, whereas CICSAM has severe bounding problems even on the finest mesh.
Sharpness. Table 1e shows our sharpness measure, δWrel. For all simulations, the isoAdvector thickness is very close to the best one can expect, i.e. the thickness of the volume fraction representation of the exact solution on the given mesh. The MULES interface width is only 30–50% larger than the width of the exact solution. HRIC performs rather badly in terms of interface sharpness with a smearing of the interface which is clearly visible in figure 4 (column 3). CICSAM keeps the interface sharp for all runs and is the best performing of the reference schemes in this respect.
Efficiency. From table 1f, we see that, for this simple test case, isoAdvector is slightly faster than the fastest reference schemes, CICSAM and HRIC, for most simulations, and two to four times faster than MULES. It is remarkable that the isoAdvector scheme can obtain a significantly improved accuracy with this significantly lower usage of computer resources.
4.1.2. Unstructured meshes
In figures 5 and 6, we show a sequence of simulations similar to those in figure 4, but now on triangular and polygonal meshes, respectively. Again the columns show ( from left to right) the solutions obtained with isoAdvector, MULES, HRIC and CICSAM. From row 1 to 2 we refined the mesh, keeping the Courant number at =0.5. From row 2 to 3 we retain the mesh, but go from Co=0.5 to 0.1. As the meshes have no preferred direction, we use velocity u=(1,0) for these simulations. The disc radius is still R=0.25, and the solutions are shown at time t=4. Inspection of figures 5 and 6 and the quantitative measures (here only E1 is shown in tables 2 and 3) reveals that most of the observations listed above for the square mesh also hold for the triangle and polygon meshes. There are, however, a number of differences concerning the performance of the reference schemes.
First, the tendency of MULES to align the interface at 45 degrees with the mesh faces is no longer visible due to the random face orientations, which presumably causes this systematic error to cancel out. However, on the triangle mesh, MULES still does not seem to converge to a circular interface due to the development of ‘wings’ on the sides (relative to the flow direction) of the fluid A region. On the polygon mesh, MULES does significantly better in terms of shape preservation, though with a tendency to squeeze the interface along the direction of motion.
Second, HRIC is much better at preserving the interface shape on both the triangle and polygon mesh than on the square mesh. It is, however, still very diffusive. This is a good example of a case where the α=0.5 contour (blue) alone would give an impression of good performance, but where the α=0.01 and 0.99 contours (green) reveal the excessive smearing of the interface.
Third, CICSAM performs very poorly on the triangle mesh with threads of fluid B piercing into the disc volume from behind. On the polygon mesh, these threads are not present and the solution quality is similar to the square mesh solution. On both the triangle and the polygon mesh, CICSAM has the same problems with unboundedness that we saw on the square mesh.
We conclude that also on unstructured meshes in two-dimensional the performance of isoAdvector is significantly better than the reference scheme with calculation times that are similar to HRIC and CICSAM, and significantly lower than MULES.
4.2. Spiralling disc
After two-dimensional uniform flow tests, our next step is to test the solver performance in a spatially varying flow. We adopt the set-up shown in figure 7, which has become a standard case for testing the ability of an interface advection schemes to deal with severe interface stretching [17,19,20–25]. The domain is the unit square with a disc of radius R=0.15 initially placed at (x,y)=(0.5,0.75). The velocity field is given by 4.3where the period of the flow is set to T=16. This flow stretches the disc into a long filament until, at time t=4, the flow is completely attenuated by the temporally varying cosine prefactor. Then, the flow reverses and the volume of fluid flows back into its original shape at time t=8. At this time, our shape preservation error measure, E1, can be calculated by comparing the computed final state with the initial state. As we know the flow in advance, we use the fixed intermediate velocity, u(x,y,t+0.5Δt), on the whole time interval [t,t+Δt].
In figure 8, we show the square, triangle and polygon meshes in three different resolutions on which the isoAdvector method was tested. The results are shown in figure 9 using the same arrangement of the meshes. All simulations are run with Co=0.5. In each panel, the exact initial and final interface shape is shown with a red circle overlaid with the α=0.5 contour (blue) of the final (i.e. at time t=8) volume fraction data. The spiral-shaped volume of fluid at time t=4, where it is maximally stretched, is also shown in each panel.
All simulations show some degree of pinching at t=4. This occurs when the filament thickness reaches the cell size as is to be expected. The phenomenon is therefore most pronounced on the coarsest meshes. We note that although the exact mathematical solution does not pinch, the 0.5-contour of its volume fraction representation will indeed pinch if the mesh is coarse enough. As such, pinching does not have to be an error. However, as droplets pinch off, and the local interface curvature becomes comparable with the cell size, the isofaces are not able to represent the significant interface curvature inside a cell. The isoface-based approximation of the advection then becomes faulty, leading to errors in the estimate of the droplet motion similar to those reported in Cerne et al. . The irreversibility of the introduced errors causes a distortion of the final disc in its upper region, which is made up of the previously pinched-off fluid.
The mesh sizes, error measures and calculation times are shown in table 4. From the E1 values in table 4, the orders of convergence with mesh refinement are calculated to be 1.9, 1.7 and 1.9 for the square, triangle and polygon meshes, respectively.
For a comparison, we show in figure 10 and table 5 the results obtained with MULES on the intermediate resolution meshes of the three types, using Co=0.1. For the square mesh, the MULES E1 error is ∼50% larger than the corresponding isoAdvector error. For the triangle mesh, the final interface is completely disintegrated. On the polygon mesh, MULES also gives acceptable results, although the E1 error is five times larger than the isoAdvector error on the same mesh. In terms of calculation times, MULES is ∼10 times slower than isoAdvector. This is in part because MULES is run with smaller time steps. However, we also ran the simulations with Co=0.5, in which case the MULES results on all three meshes were completely disintegrated like the triangle mesh solution in figure 10.
4.3. Sphere in steady uniform three-dimensional flow
In this test, we go back to a uniform flow, but now in three dimensions. The velocity is U=(0,0,1), and the initial interface is a sphere of radius R=0.25 centred at (0.5,0.5,0.5). The simulations are run on three meshes consisting of 49 868, 343 441 and 1 753 352 random tetrahedral cells covering the domain, [0,1]×[0,1]×[0,5]. The meshes and the 0.5-isosurface of the initial volume fraction data are shown in figure 11. The simulations are run with Co=0.5 until t=4, where the sphere has moved to (0.5,0.5,4.5). The results are shown in figure 12 and table 6. In figure 12a, we show the exact final sphere (red) and the 0.5-isosurface of its volume fraction representation on the three mesh resolutions. In figure 12b, we show the exact sphere (red) together with the 0.5-isosurface of the final volume fraction data obtained with isoAdvector. As seen from table 6, the E1 error on the coarsest mesh is fairly large. From figure 12 (lower left panel), we see that this lack of overlap is mainly due to an overestimation of the propagation speed rather than a lack of ability to retain the spherical interface shape. On the finer meshes E1 is reduced significantly, although the tendency to be slightly ahead of the exact solution is still visible in figure 12. The linear cell size is reduced by a factor of 1.9 from the coarse to intermediate mesh, and by a factor of 1.7 from the intermediate to fine. Based on these ratios and the E1s in table 6, the convergence order is calculated to lie in the range 2.6–3.2.
For comparison, we show in figure 13 and table 7 the results obtained with MULES on the finest mesh running with Co=0.1 and 0.5. In both cases the shape preservation is significantly worse than the isoAdvector results. It is also notable that the MULES simulations with Co=0.1 and Co=0.5 are, respectively, 20 and 5 times slower than the corresponding isoAdvector simulation with Co=0.5.
4.4. Sphere in non-uniform three-dimensional flow
Our final test case is also in three dimensions, but now with a non-uniform velocity field. We adopt a configuration from LeVeque , which has become a standard test case [19,28–30] for testing interface advection methods and their ability to deal with highly distorted interfaces in three dimensions. The domain is the unit box, and the initial interface is a sphere of radius R=0.15 centred at (0.35,0.35,0.35). This surface is advected in the velocity field, 4.4where the period is set to T=6. This flow stretches the sphere into a thin sheet creating two bending and spiralling ‘tongues’. The maximum deformation is reached at t=1.5, where the temporal cosine prefactor completely quenches the flow. From here on the flow reverses, and the interface returns to its initial shape and position at time t=3. In figure 14, the isoAdvector results are shown at time t=1.5 in the top row, and at time t=3 in the bottom row, on three cubic meshes with and . In the lower panels, the exact final spherical shape is also shown in red. From ODE calculations with the velocity field (4.4), we have measured the sheet thickness at t=1.5 to be ∼0.0063. This, and the fact that an edge can at most be cut once by the prescribed isosurface routine, explains why there are holes in the 0.5-isosurface of the volume fraction data on the two coarsest meshes with dx≈0.016 and dx≈0.0078, and no holes in the finest simulation with dx≈0.0039. The error measures and calculation times for the three simulations are shown in table 8. Based on the E1s in this table, the order of convergence is calculated to be 2.2.
Error measures for this test case for three different geometric VOF approaches are listed in table 7 of Liovic . The error is given as the L1-norm, which is related to the E1 measure by L1=E1V A, where V A=4/3πR3≈0.0141 is the volume of fluid A in the domain. Thus, the isoAdvector E1 values in table 8 for the (64,128,256) meshes translate into L1=(3.0×10−3,6.4×10−4,1.2×10−4). Comparison with table 7 in Liovic  shows that all listed errors are similar in magnitude, with the exception that, on the finest mesh, the error of the method denoted by ‘CVTNA+PCFSC unsplit’ in Liovic  is approximately three times smaller than the other approaches in that paper and almost twice as small as our error. We note that the schemes described in Liovic  are developed for orthogonal hexahedral meshes. Our error measures are also similar in magnitude to those presented in table IV of Xie et al. , which shows the L1-norm for their method (developed for unstructured meshes consisting of hexahedra), and those described in Hernandéz et al.  (hexahedral meshes) and Jofre et al.  (general meshes).
As a final test, we have repeated this simulation on a mesh consisting of random tetrahedra. To get sufficient resolution to avoid holes in the 0.5-isosurface of the solution, we used a mesh with 10 131 041 cells covering [0.15,0.90]×[0.15,0.80]×[0.15,0.80]. A cut through this mesh and the 0.5-isosurface of the volume fraction representation of the initial sphere are shown in figure 15a. In figure 15b, we show the isoAdvector solution at time t=1.5, where the stretching is maximal. This panel also contains a solution obtained by integrating the velocity field with a Runge–Kutta ODE solver for 160.000 points evenly distributed on the initial sphere (green dots). The visual impression from figure 15 is that there is a good match between the ODE and the isoAdvector solutions. At time t=3, the shape was preserved with an error of E1=0.014. Owing to the clipping procedure which was activated in the bounding step, δV rel was −0.63%. Bounding errors were in the order of 10−5. The simulation took 67 h on a single core and 23 h on eight cores. As the parallelization of the code is based on domain decomposition, the scalability depends on the extent to which surface cells are evenly distributed among the processors.
We have developed a new algorithm, isoAdvector, for numerical interface advection across general structured and unstructured computational meshes. The method is derived from ‘first principles’, i.e. from the control volume integrated continuity equation for a discontinuous density field. The isoAdvector scheme belongs to the class of geometric VOF methods, but with novel ideas implemented in both the interface reconstruction step and in the interface advection step.
The novelty in the reconstruction step is the usage of efficient isosurface calculations to estimate the distribution of fluids inside computational cells. This is a very robust method even on unstructured meshes. It avoids the gradient calculations traditionally used in the geometric VOF reconstruction step, which may cause problems, because the numerically estimated gradient is a cell volume averaged Dirac δ-function.
In the advection step, the approach taken in existing geometric VOF methods for general meshes is to calculate so-called flux polyhedra and their intersection with the submerged part of the mesh cells [10,11,13]. This is an expensive and complicated procedure, which we avoided in our work. Instead, we focus on the time evolution of the submerged part of a face during a time step. We approximate the face–interface intersection as a line sweeping the face during the time step, and divide the time step into sub-time intervals defined by the times at which the line passes the face vertices. On such time intervals, we can then analytically calculate the volume of fluid A passing a face under the assumption that the line is moving steadily across the face during the interval. In the development of this procedure, no assumptions are made on the shape of a face, and therefore also the advection step is by design applicable on general meshes.
We have given a proof-of-concept by testing the isoAdvector method on various simple flow–interface combinations both in two-dimensional and three-dimensional structured and unstructured meshes. The results are very satisfactory both in terms of shape preservation, volume conservation, boundedness, interface sharpness and efficiency. The order of convergence with mesh refinement varies between 1.7 and 3.2 for the test cases presented here. Also, in spite of the geometric nature of some of the steps involved, the implementation of the new algorithm is relatively straightforward.
The isoAdvector advection step relies on local geometric considerations based on available information from a surface cell and its nearest neighbours. Thus, while there seems to be no strict limit of Co<1, the underlying geometric considerations will become questionable if the time step is so large that the interface moves across many cells during one time step. We therefore do not expect the method to perform well with interface Courant numbers significantly higher than 1. For Co in the range [0,1], it is our experience that the solution quality is fairly stable, albeit with an optimum around 0.5. This is to be contrasted with the explicit MULES scheme in OpenFOAM®’s interFoam solver, which in our experience is limited to Co≤0.1 if accuracy is important.2
We note that because the governing equation we solve is the passive advection equation for a scalar field in a solenoidal velocity field, the isoAdvector method may also find applications within other branches of CFD, where the advected surface is not necessarily marking the interface between two distinct fluids. There are many situations where one needs to follow a passive tracer field, e.g. representing the concentration of some substance, which is immiscible with the surrounding fluid. Another possible application could be in an immersed boundary method, where the isoAdvector scheme could provide accurate estimates of the fluid–solid interface within computational cells.
Based on the interFoam solver in OpenFOAM®, we are currently working on a consistent coupling of isoAdvector with a pressure–velocity solver. The performance of the resulting new interfacial flow solver will be presented in a future paper. Finally, we note that due to its applicability on arbitrary meshes, the isoAdvector code can be coupled with an adaptive mesh refinement routine with only minor modifications. Such a coupling will also be investigated in future work.
The isoAdvector code is published  as an open source extension to OpenFOAM®. The code has been parallelized based on the domain decomposition approach of OpenFOAM®, and we are currently modifying the code to work also on moving meshes. It is our hope that the isoAdvector concept and code will be used, tested and further developed by the CFD community, and eventually result in improved simulation quality in the broad field of applications involving interfacial flows.
The isoAdvector code is publicly available in the repository www.github.com/isoadvector together with the set-up files necessary to reproduce the results of all test cases included in this paper. The code is an extension to the open source C++ code library, OpenFOAM, which can be downloaded from www.openfoam.com. The unstructured meshes used in this paper are available in the OpenFOAM polyMesh format from the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.66840 .
J.R. conceived the isoAdvector method; implemented it in OpenFOAM®; set up, conducted and analysed the numerical experiments; and drafted the manuscript. H.B. participated in the conception of the isoAdvector scheme and in the manuscript drafting. H.J. contributed to the implementation of the isoAdvector scheme in OpenFOAM®, implemented the HRIC and CICSAM schemes in OpenFOAM®, and revised the manuscript critically for important intellectual content. All authors gave their final approval for publication.
We have no competing interests.
This work was sponsored by a Sapere Aude: DFF—Research Talent grant from The Danish Council for Independent Research — Technology and Production Sciences to J.R. (Grant–ID: DFF – 1337-00118). The grant also covers all activities of H.B. and H.J. in connection with the project. J.R. also enjoys partial funding through the GTS grant to DHI from the Danish Agency for Science, Technology and Innovation. We would like to express our sincere gratitude for this support.
J.R. is grateful to the following people for fruitful discussions and for help on improving the isoAdvector code: Tomislav Marić, Daniel Deising and Holger Marschall from the Mathematical Modeling and Analysis Group at the Center of Smart Interfaces, Technische Universität Darmstadt; Vuko Vukčević and Tessa Uroić from the Faculty of Mechanical Engineering and Naval Architecture, University of Zagreb, Bjarne Jensen from the Department for Ports and Offshore Technology, Danish Hydraulic Institute; and Henrik Rusche from Wikki Ltd.
↵1 On the surface one could set for ρ to be defined everywhere. However, as has zero volume, the value of ρ on is immaterial.
↵2 It should be mentioned that, from OpenFOAM® v. 2.3.0, a new semi-implicit MULES scheme is introduced to solve some of the issues with boundedness, stability and accuracy for large Courant numbers. The literature documenting the ideas behind and the performance of this new method is, however, still very sparse.
- Received December 19, 2015.
- Accepted October 28, 2016.
- © 2016 The Authors.
Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.