Iterative Volume-of-Fluid interface positioning in general polyhedrons with Consecutive Cubic Spline interpolation

A straightforward and computationally efficient Consecutive Cubic Spline (CCS) iterative algorithm is proposed for positioning the planar interface of the unstructured geometrical Volume-of-Fluid method in arbitrarily-shaped cells. The CCS algorithm is a two-point root-finding algorithm specifically designed for the VOF interface positioning problem, where the volume fraction function has diminishing derivatives at the ends of the search interval. As a two-point iterative algorithm, CCS re-uses function values and derivatives from previous iterations and does not rely on interval bracketing. The CCS algorithm only requires only two iterations on average to position the interface with a tolerance of $10^{-12}$, even with numerically very challenging volume fraction values, e.g. near $10^{-9}$ or $1-10^{-9}$. The proposed CCS algorithm is very straightforward to implement because its input is already calculated by every geometrical VOF method. It builds upon and significantly improves the predictive Newton method and is independent of the cell's geometrical model and related intersection algorithm. Geometrical parametrizations of truncated volumes used by other contemporary methods are completely avoided. The computational efficiency is comparable in terms of the number of iterations to the fastest methods reported so far. References are provided in the results section to the open-source implementation of the CCS algorithm and the performance measurement data.

An essential task of any multiphase simulation method is the identification and tracking of immiscible fluid phases, schematically shown in fig.1a.Each phase with distinctive physical properties can be associated with a subset Ω ± (t) of the domain Ω = Ω + (t) ∪ Ω − (t), separated by the fluid interface Σ(t) with outward-oriented normal vectors n Σ .An indicator function is introduced to distinguish Ω ± (t), as It is impossible to calculate χ (t, x) analytically when the interface motion is driven by a complex two-phase flow velocity.Therefore, approximations such as those shown in fig.1b are used, specifically by the geometrical Volume-of-Fluid method, to approximately calculate Ω ± (t).The domain Ω is discretized as a union of non-overlapping polyhedrons V c , such that Ω = ∪ c V c , c ∈ C. For each polyhedron V c , a volume fraction is defined as or, in other words, a fill level of the polyhedron V c with 0 ≤ α c ≤ 1.The Piecewise Linear Interface Calculation (PLIC) [7,8,9] is still prevalently used by the geometrical VOF method to approximate the indicator function χ (t, x) in each V c using where H c (p c (t), n c (t)) is the positive halfspace of the PLIC plane given by the plane position p c and plane normal n c at time t.The interface reconstruction algorithm computes the halfspace H c (p c (t), n c (t)), by approximating the interface unit normal vector n c and computing the point p c , such that the intersection between the positive halfspace and the polyhedron V c satisfies The computation of p c in each V c is the interface positioning part of the reconstruction algorithm, given {n c } c∈C , {α c } c∈C , under the condition 0 < α c < 1, and a polyhedral domain discretization Ω = ∪ c V c , c ∈ C. Since the goal of the interface positioning is to find p c (t) at a fixed t, t can be disregarded from now on.An explicit (closed-form) solution of eq. ( 4) for p c does not exist for an arbitrary volume V c .For example, the positioning of the PLIC plane in an endo-dodecahedron1 from fig. 2a is shown in fig.2b with α c = 0.5 and the normal vector n c that is collinear with the Z-axis.The resulting intersection between the positive halfspace of the PLIC plane H c (n c , p c ) and the endo-dodecahedron is the shaded volume in fig.2b.To simplify the solution process for p c in eq. ( 4), the position of p c on the n c axis can be parameterized with a scalar s, as shown in fig.2b.Minimal and maximal point x min,max of the volume V c , projected onto n c , as shown in fig.2b, are defined as where x ref is a reference point for the s-axis.Any point can be chosen as x ref and usually the origin of the coordinate system is used.However, using the origin as a reference point increases the error of floating-point calculations used in eq. ( 5), and using a point that belongs to V c (e.g. the first point, or the centroid of V c ) reduces these errors [10, section 2].The set P in eq. ( 5) is the set all points x p of the volume V c .Using eq. ( 5), the halfspace position p c is parameterized with where s ∈ [s min , s max ].Note that the parametrization determines the values of the parameter s: for example, with the x min , x max , parametrization by eq. ( 7), s ∈ [0, s max ].This, in turn, leads to the parametrization of the intersection volume in eq. ( 4) as The parameterization of the intersection volume reformulates eq. ( 4) as where α c is the given volume fraction for V c .The root s * of α(s) is sought to compute p c (s * ), the position of the PLIC plane (p c (s * ), n c ) in V c such that eq. ( 4) is satisfied, as shown in fig.2b.The function α(s) is visualized in fig. 3 for the endo-dodecahedron shown in fig.2a, using n c = (0, 0, 1) and interface positioning sometimes separates V c into disjoint volumes, similar to the result shown in fig.3b for s ≈ 2. Geometrical operations used for the intersection must therefore support disjoint sets.Non-convex cells with non-planar faces are decomposed into tetrahedrons using the centroid of the cell and centroids  of the non-planar faces.The magnitude of the truncated volume that lies inside the halfspace H c (p c (s), n c ) is, therefore, equal to the sum of volume magnitudes of the truncated tetrahedrons.The tetrahedral decomposition is the most straightforward approach to the volume truncation of non-convex cells with non-planar faces [11], but it is more computationally expensive, compared to sophisticated approaches such as [6].The choice of the truncation is arbitrary for the proposed CCS method: an alternative volume truncation reduces the total CPU time, however, it has no influence on the total number of iterations, whose very straightforward reduction the main contribution of the proposed method.
Following boundary conditions are given for α(s): Boundary conditions given by eq. ( 13) complicate the root finding in eq. ( 9), because the derivatives of α(s) are diminishing at s = s min , s max .Furthermore, although highly improbable, it is possible that the normal vector n c becomes collinear with the normal of a planar face of the cell.In this case either α (s min ) = 0 or α (s max ) = 0.This special case of the boundary condition is handled by the proposed CCS method and addressed in detail below.
Positioning algorithms can be categorized as iterative or bracketing algorithms.Iterative algorithms rely solely on root-finding to solve eq. ( 9) up to a prescribed tolerance.Bracketing algorithms intersect the volume V c incrementally, until a bracketed volume is found that contains the interface.In the bracketed volume, modern positioning algorithms [4,5,6] exactly position the interface using an exact function for the volume fraction within the bracketed interval.
Kothe et al. [12], Rider and Kothe [9] have proposed an iterative algorithm that relied on Brent's root finding method, combining inverse quadratic interpolation with bisection to avoid divergence near interval boundaries.Bracketing intersects the volume V c with the halfspace given by H(p p , n c ), p ∈ P , where {p p } p∈P are cell corner points.This brackets a volume V c as resulting in the bracketed volume B c,i,j that contains the root s * of α given by eq. ( 9).Within the bracketed volume, Brent's algorithm requires additional iterations to find the root of eq. ( 9) up to a prescribed tolerance.This procedure was widely adopted for interface positioning [13,14,15,16].Scardovelli and Zaleski [17] have proposed an analytical interface positioning method for rectangular meshes that (used in [18]) and Yang and James [19] have proposed an analytical positioning method for tetrahedral and triangular meshes.López and Hernández [20] have extended the work of Scardovelli and Zaleski [17] and Yang and James [19].They have used and indexed face-set2 as the boundary representation of a polyhedron Ghali [21,Chapter 28] and an analytical expression for a volume of a convex polyhedron [22] to position the interface exactly within the bracketed interval.López and Hernández [20] introduce the Central Sequential Bracketing (CSB) procedure that utilizes sorted signed distances associated with the corner points i p of the volume V c , calculated using the interface normal n c .The goal of the algorithm is to locate indices k min , k max in this list, such that  14 I p + 1) and O(log 2 I p ) if so-called Binary Bracketing (BB) is used to set k c = (k min + k max )/2, where I p is the total number of points of the volume V c .They utilize BB for cells with I p > 8 because then the logarithmic complexity outperforms the linear complexity of the CSB bracketing algorithm.López et al. [5] have published the source code that implements their methods from [20].
Ahn and Shashkov [23] have proposed a stabilized bisection-secant method as an iterative positioning algorithm, used by the author in [24,25].Their algorithm benefits from the super-linear (golden ratio) convergence of the secant method outside of the regions with diminishing derivatives, and relies on the bisection method to ensure convergence otherwise.They rely on the tetrahedral decomposition of the volume for the calculation of the truncated volume, however the root finding algorithm is independent of that choice.
Diot et al. [3], Diot and François [4] have proposed an analytical positioning method that does not rely on the analytical expression for the volume of the convex polygon (polyhedron).Instead, the polygon (polyhedron) is decomposed into sub-volumes whose magnitudes are exactly calculated using mixed vector and scalar products.The proposed method is very accurate, however it involves relatively complex subdivisions of the V c volume.
López et al. [26] present a detailed review of interface positioning algorithms and propose an enhancement of their CSB and BB bracketing algorithms from [20] on convex polyhedral cells.Their new Interpolation Bracketing (IB) algorithm uses signed distances to linearly interpolate H(p c (s), n c ), in effect interpolating linearly the volume fraction function shown in fig.3b.The authors couple their IB algorithm with a new explicit function for interface positioning in the bracketed interval.The authors state that the explicit positioning function costs as much as 1.7 times the volume truncation operation, which must be taken into account when its overall computational complexity is considered.The new explicit function is combined with the IB algorithm into the new Coupled Interpolation-Bracketed Analytical Volume Enforcement (CIBRAVE) method.
Recently, Chen and Zhang [2] have noticed that the fundamental theorem of calculus can be used to increase the convergence of the iterative positioning algorithm.In fig.4, the intersection of the plane (p(s), n c ) and the volume V c is the polygon with the area A(s) that lies on the plane (p(s), n c ): the so-called cap polygon.The area of this polygon (cap area) is computed geometrically, with machine-epsilon accuracy, at any s, and it can be used to integrate parametrized volume V (s) as It is obvious that A(s) = V (s) is the exact derivative of the volume V (s) from eq. ( 9), whose root determines the position of the interface p(s * ).A direct consequence of this is Chen and Zhang [2] use the exact derivative α (s) in the Newton's root finding method to solve eq. ( 9).Knowing the exact derivative ensures the quadratic convergence of the Newton's method, provided that the condition α (s) = 0 is fulfilled.However, it does not solve issue of divergence of the Newton's method near interval boundaries s min , s max .If an iteration of the Newton's method computes a root outside of (s min , s max ), Chen and Zhang [2, Algorithm 1] apply the Hermite cubic spline interpolation at s, across (s min , s max ).Unfortunately, interpolating across the whole search interval (s min , s max ) introduces large interpolation errors and substantially slows down convergence for volume fraction values α c ≈ 0, α ≈ 1.To understand why, consider the case when α c ≈ 1: the derivative α (s) ≈ 0, Newton's method shoots out of the search interval s min , s max , and then Hermite interpolation over s min , s max is performed.If α ≈ 1, the root lies near s min .However, the Hermite interpolation introduces the conditions on α(s) at s max , that lies on the other end of the search interval.Therefore, the accuracy is lost in this case, and vice versa for α c ≈ 0. The method of Chen and Zhang [2] is named Newton Cubic Spline (NCS) method to facilitate comparison with the proposed method.Similar to Brent's method, which is a multipoint method3 , the algorithm proposed here reuses function values α(s) and derivatives α (s) from previous iterations to effectively double the order of accuracy of the interpolated volume fraction.This is the basis of the proposed Consecutive Cubic Spline (CCS) interpolation.An additional stabilization for the ends of the interval [s min , s max ] is developed, that effectively handles diminishing derivatives at the ends of the search interval.
Chen and Zhang [2] increase the computational efficiency of the Moment-of-Fluid Method [27], by expressing the centroid of the cap polygon A(s) as a function of spherical interface orientation angles θ, φ.This is relevant for increasing the efficiency of interface positioning in the context of algorithms that improve the orientation of the VOF interface.The CCS algorithm proposed in this manuscript achieves significantly higher computational efficiency compared to NCS, especially for challenging values α c ≈ 1, α c ≈ 0, even without utilizing the information about the interface orientation.Of course, CCS can be coupled to MoF or any other algorithm that improves interface orientation in the same way as NCS.Since the convergence order of the CCS algorithm, demonstrated in the following section, is significantly higher than NCS without the use of interface orientation, adding this information would only increase the computational efficiency of the interface reconstruction algorithm.In future work, CCS will be coupled with the simplified Swartz reconstruction algorithm [25].

Chen and Zhang
known to be second-order convergent if α (s n ) = 0 within the search interval, which is not the case for interface positioning, where the function derivatives diminish on the boundaries of the search interval.
The proposed Consecutive Cubic Spline (CCS) algorithm interpolates the volume fraction α(s) with a Hermite polynomial between two consecutive iteration steps when they contain the root, by re-using volume fraction values and derivatives from previous iterations.
Using function value and function derivative from two subsequent iterations makes the CCS a two-point root finding algorithm [1, chap. 2].If the interpolation interval does contain the root, the cubic interpolation exactly recovers the volume fraction α(s), and its root is therefore the root of α(s).The volume fraction α(s) is approximated with the polynomial of order K < 4, P k (s) = k=0...K a k s k ≈ α(s), using the conditions that determine k, if s * ∈ [s a , s b ] ⊂ (s min , s max ), where s * is the root of α(s) given by eq. ( 9).
To understand why P k (s) recovers the volume fraction exactly, consider the following.The volume V c is bounded by planar polygons or, if the boundary of V c consists of non-planar faces, by sets of triangles resulting from the triangulation of non-planar faces.In both cases, the cap area A(x = x min + sn c ) in fig. 4 is at most quadratic in (x, y, z) in three dimensions.Because the parametrization x = x min +sn c given by eq. ( 7) is linear, parametrized A(s) is then also at most quadratic in s, in three dimensions.By eq. ( 15), α(s) is then at most cubic in three dimensions.We also know that α(s) is piecewise-polynomial between the points x p of the volume V p : it is C 1 continuous at the boundaries of intervals bracketed by the points x p of V c , because its derivative, the cap (intersection) area A(s), is C 0 continuous in s, in cells that are used in unstructured meshes for the discretization of Partial Differential Equations.Therefore, α(s ] with ((α(s 4 ), α (s 4 ), (α(s 3 ), α (s 3 )).
Figure 5: A schematic representation of the CCS method in the worst case scenario: when the root s * happens not to lie between two successive iterations.Because of the high accuracy of the spline P k (s), consecutive iterations are very close to each other and immediately "land" within a bracketed interval.Within this interval P k (s) exactly recovers α(s), and thus the interface position.
inside the interval bracketed by points x p of the volume V c , and it is C 1 continuous on the boundaries of bracketed intervals.The α(s), defined by eq. ( 9), is calculated geometrically using the truncated volume V c .Similarly, α (s) is geometrically calculated as A(s), the cap area in fig. 4. Therefore, the conditions given by eq. ( 18) are satisfied exactly up to machine epsilon.Consequentially, P k (s) exactly interpolates α(s) within the bracketed interval and an exact root of P k (s) is equivalent to the root of α(s) within the bracketing interval.Note that the calculation of α (s) comes at no additional computational cost, because A(s) is a byproduct of the volume truncation used to compute V (s) for α(s).
This very straightforward way of interpolating α(s) that exactly recovers α(s) within a bracketed interval is the main contribution of the CCS algorithm, compared to other methods that require complex geometrical parameterizations [3,4,5,6] to achieve the same.Furthermore, directly compared to NCS [2], CSS requires overall significantly less iterations, especially at the ends of the search interval, where the volume fraction derivative diminishes.The two-point root finding approach is novel compared to only using α (s) in a Newton iteration by the NCS algorithm of Chen and Zhang [2].The number of iterations of the CCS algorithm are not impacted by the way α(s), α (s) are computed.Different geometrical models and algorithms can be used to compute α(s), α (s).As long as α(s), α (s) are computed near machine-epsilon accuracy (e.g. using algorithms from [6], or from [11]), the number of iterations performed by CCS will not change.Because the quadratic nature of the intersection area results in eq. ( 19), extending the two-point CCS approach to include a polynomial interpolation with order K > 3, by storing data from more than two consecutive iterations brings no benefit.In fact, increasing the polynomial order beyond 3 destabilizes the root finding process because higher-order polynomials often have more than a single root s * in the search interval (s a , s b ).
The first part of the proposed CCS iterative positioning method is shown in algorithm 1 and the CCS interpolation part is shown in algorithm 2. The initial guess of the CCS method is calculated using the Hermite cubic interpolation over [s min , s max ], with boundary conditions given by eq. ( 18) at s min , s max , as shown schematically in fig.5a.To use the CCS interpolation between consecutive iterations in next iterations, the interval [s a , s b ] is updated in algorithm 1 with data from the current iteration s n (s 2 in fig.5a).
In the first iteration, there is only a single root candidate available, shown as s 2 schematically in fig.5b.Hermite cubic spline is therefore interpolated in the next step over the initial interval [s min , s max ] ([s 0 , s 1 ] in fig.5b), together with the first root candidate s n and its value and derivative pair (α(s n ), α (s n )).
For all further iterations (n > 0), as shown in algorithm 2, the CCS interpolation with P k (s) is performed between [s a , s b ] using values and derivatives from previous root candidates.If [s a , s b ] is a subset of a bracketing interval, and the root of P k (s) is inside [s a , s b ], the position of the interface (root of α(s)) is exactly equal to the root of P k (s).For the results presented in this manuscript, Brent's method [28] is used to find the root of the polynomial with machine epsilon tolerance.Finding the root of P k (s) using the Brent's method costs only a small fraction of overall computational cost and the choice of the polynomial root-finding algorithm has no impact on the number of iterations of the CCS method as long as the polynomial root is calculated within a machine-epsilon tolerance.p(s n ) = x min + s n n c 13: end if An important detail regarding P k (s) is the accurate calculation of of the polynomial coefficients for P k (s).An error in the calculation of P k (s) coefficients means P k (s) cannot exactly interpolate α(s) in [s a , s b ], resulting, in turn, with an erroneous interface position, even if the root of P k (s) is calculated exactly.Parameterizing s such that s ∈ [0, 1], enables a numerically stable way to exactly calculate the coefficients of P k (s) as An unnecessary intersection is removed by line 37 in algorithm 2, in the case when the Algorithm 2 Consecutive Cubic Spline (CCS) iterative positioning: Part II

27:
if n == 0 then A single root candidate is available in the first iteration.28: Interpolate P k (s) over [s 0 , s n , s 1 ].

45: end if 46:
end if 47: end for root value of P k (s) falls below the positioning tolerance, because when s * ∈ [s a , s b ], P k (s) recovers α(s) exactly.
In the case when s * / ∈ [s a , s b ], a Newton iteration is used, as shown schematically in fig.5c.A separate special case of s * / ∈ [s a , s b ], happening over multiple iterations, is caused by the collinearity of the halfspace normal n c with the normal vector of a planar face of V c .This leads to α (s b ) = 0 (α (s a ) = 0) where s b = s max (s a = s min ), which is contrary to the boundary conditions 18.This special case is handled in a very straightforward way by CCS, by lines 40 − 42 in algorithm 2. In the case where CCS needs more than 2 iterations and one of the boundary derivatives is zero, this is a possible special case with a falsely assumed zero derivative at one of the interval boundaries.In this case a simple bisection step is performed, because it will force a truncation and thus the calculation of the actual α (s n+1 ), replacing the assumed zero derivative at s min or s max with the calculated derivative at s n+1 .

Results
The CCS positioning algorithm is tested in cells with different shapes, using challenging volume fraction values and interface orientations, to ensure robust con-vergence on unstructured meshes.The cells used for testing are shown in fig.6, namely: the tetrahedron (TET), the unit cube (CUBE), the dodecahedron (DOD), the endo-dodecahedron (ENDO) and the non-convex dodecahedron with non-planar faces (NPDO).The endo-dodecahedron represents a model of a non-convex polyhedron that has planar faces.The non-planar dodecahedron is ubiquitous in unstructured polyhedral meshes that are generated by agglomerating tetrahedral cells into polyhedral cells.
The CCS positioning is compared to NCS [2] using following test parameters where (φ, θ) are the spherical angles that parametrize the interface normal as Values of α c in eq. ( 27) are based on the reconstruction tolerance R = 1e−9, proposed in [23] and used by the author in [24,25].It is crucial to include values boundary volume fraction values (α c ≈ 0, α c ≈ 1) in the test data, because positioning algorithms that rely on nonzero derivatives fail at the boundaries of the α c ∈ [0, 1] interval.Boundary volume fractions are common in multiphase flow simulations, where the so-called wisps (dimensionally unsplit VOF) and splitting-errors (dimensionally split VOF) introduce small variations in α c in cells that are otherwise either full (α c = 1) or empty (α c = 0).Since the derivatives of α c (s) vanish for these volume fraction values (cf.fig.2b), a positioning algorithm that does not take into account diminishing derivatives will experience either slower convergence, or in the worst case, divergence.The positioning condition is based on the positioning tolerance R < 1e − 12, also used by López et al. [5]: this forces the CCS and NCS iterative algorithms to perform one more iteration, which further reduces P below machine epsilon for the vast majority of tests.The test results presented here are generated with N θ = 40 and N φ = 20, because there is no correlation observed between the interface orientation and the number of iterations or CPU time.Furthermore, as the tests results confirm, there is no need to use many α c values between [1e − 03, 1 − 1e − 03], as the average iterations do not vary significantly within this interval.That is why N α = 50 was used for the tests presented here.First let us consider the overall results of the CCS algorithm, compared with the NCS algorithm [2].The increase in efficiency in terms of the average number of iterations is presented in fig.7a, and the decrease in the average CPU time is shown for different cell shapes in fig.7b.An iteration is defined in this context as the calculation of the next interface position s n+1 .Timing was performed on an architecture given by table 1, and the measurements were performed using the chrono C++ standard template library.The same tests were performed on a High-Performance Computing cluster without a significant difference in reported results.A Singularity image [29] that contains the source code and the computing environment is publicly available [30] and can be used to easily reproduce the results across different platforms.Results presented in this section using the computing architecture from table 1 are also publicly available, together with the source code and binary executables [31] 4 .The results presented in fig.7a show that CCS algorithm requires approximately 2.9−3.7 times less iterations than the NCS algorithm for different polyhedron shapes.In terms of the CPU time shown in fig.7b, the CCS algorithm is approximately 1.7−3.0times faster than the NCS algorithm for different polyhedrons.The stabilized secant / bisection of Ahn and Shashkov [23] and Brent's method were not used in this comparison, because Chen and Zhang [2]   Chen and Zhang [2] use the NCS algorithm in combination with the Moment-of-Fluid (MOF) method and extend the interface positioning problem with the interface orientation angles θ, φ that define the interface orientation vector in eq. ( 28).Using the implicit formulation for the interface plane n c,x x + n c,y y + n c,z z + d = 0, the interface position d is defined by Chen and Zhang [2] as a multivariate function d := d(α c , θ, φ).Chen and Zhang [2] then extend the multivariate function d into Taylor series which results in a prediction of the interface position in the new iteration based on ∂ θ d and ∂ φ d.Consequently, when the information about the changing interface orientation is available from an algorithm that tries to improve the orientation (e.g.MOF), it is possible to find the position of the interface faster.However, it is crucial to note that the orientation information does not impact the partial derivative ∂ αc d and the NCS algorithm is used in Chen and Zhang [2] independently of the interface orientation prediction.The CCS positioning algorithm improves significantly the estimation of ∂ αc d by the CCS polynomial interpolation of α(s).Chen and Zhang [2, Figure 12] report an average of 4 iterations when the orientation information is given by N θ = 100(N φ = 200), and approximately 3 iterations when the orientation information is given by N θ = 1000(N φ = 2000) for a test case with α The CCS algorithm positions the interface with only 2 average iterations without relying on the interface orientation information (cf.fig.8), for different cell shapes,  and with a more challenging set of α c values given by eq. ( 27).This is confirmed by Chen and Zhang [2, Figure 13], with more challenging values α c ∈ [1e − 09, 1e − 09], where the average iterations are reported as follows: NCS 13.53 average iterations, predicted NCS (N θ = 100) 6.57 average iterations, and predicted NCS (N θ = 1000) 2.89 average iterations.Therefore, it takes the orientation prediction in [2] 1000 subdivisions of the interval θ ∈ [0, π] to provide enough information to the Taylor series expansion in order to reach the average number of iterations given by the CCS that does not require this information.Because the MOF method requires significantly less than 1000 of iterations to improve the interface orientation, Chen and Zhang [2] report an overall increase of 60 − 66% in efficiency in terms of the average number of iterations.It is relevant to note that the implementation of NCS used here for comparison is numerically unstable: a test configuration for the nonplanar dodechedron in fig.8 causes the NCS not to converge even after 100 iterations.
The CIBRAVE method by López et al. [26,5] uses Interpolation Bracketing based on linear interpolation of the bracketing volume, and an explicit (analytical) calculation of the interface position within the bracketed interval.López et al. [26,Fig. 14] report a log 2 (log 2 (I p )) complexity for the average number of volume truncation operations for the CIBRAVE method, where I p is the number of vertices of the polyhedron.This results in 2.6 bracketing truncations on average for a dodecahedron.However, unlike the iterative algorithms, bracketing algorithms additionally require the evaluation of a geometrically parameterized explicit function to position the interface, and this comes at a cost (cf.[26,Fig 15]).López et al. [26] state that the cost of the coefficient calculation for explicit positioning function is "around 1.7 times the CPU time needed to make a truncation operation" [26].López et al. [26] use the test parameters from Diot and François [4], namely α ∈ [1e − 03, 1 − 1e − 03], which is not as challenging as volume fraction values used here for the CCS method and by Chen and Zhang [2, Figure 13] for the NCS algorithm, so the results are difficult to compare directly.Even though López et al. [26] also use different cell shapes, the average number of iterations of CCS is 2, which is comparable to CIBRAVE, while being significantly easier to implement.
The distribution of the CPU time between sub-algorithms of the CCS method is shown in fig.9.In fig. 9 the CPU time distribution is reported for every test polyhedron and a comparison between the NCS algorithm (blue color) and the CCS algorithm (orange color) is given.The CPU time is distributed between the numerical root finding for the Consecutive Cubic Spline (ROOT), the Consecutive Cubic Spline polynomial interpolation (POLY) and the geometrical volume truncation operation (GEOM).As expected, most of the computational cost is caused by the volume truncation operation, which is why the reduction in the number of iterations in fig.7a correlates well with the speedup in terms of CPU time in fig.7b.The GEOM CPU time would need to be reduced by approximately an order of magnitude for some polyhedrons in fig. 9 to reach the absolute level of ROOT and POLY times and thus impact the relative increase in computational efficiency.For the tetrahedron, this is not relevant, as the average positioning CPU time is only 4.5 microseconds.López et al. [6] extend the volume truncation algorithm used by CIBRAVE [26,5] to non-convex polyhedrons with the aim to maintain the same computational efficiency.López et al. [6] compare the computational efficiency of the CIBRAVE method with and without their novel truncation using the Brent's method as the basis for comparison.An order of magnitude faster positioning is achieved with CIBRAVE and the new truncation, and four times faster when with CIBRAVE and tetrahedral decomposition.The CCS algorithm would benefit in the absolute CPU time by introducing a faster truncation algorithm (e.g. the one proposed by López et al. [5,6]), however this does not change its relative speedup, resulting from the reduction of the number of truncations.The truncation of a cell in an unstructured mesh is closely related to the data structure used to implement the unstructured mesh.To reduce the CPU time required for the interface positioning further, by introducing a new truncation the average number of iterations even when R = 1e − 09, which has not yet been reported for other contemporary positioning algorithms, to the best of the author's knowledge.

Conclusions
A straightforward iterative algorithm is developed that significantly improves the computational efficiency of the VOF interface positioning problem on arbitrary unstructured meshes.The proposed Consecutive Cubic Spline (CCS) algorithm outperforms Brent's method, the stabilized secant-bisection method of Ahn and Shashkov [23], and the Newton Cubic Spline method by Chen and Zhang [2].The CCS algorithm is comparable with the CIBRAVE method of López et al. [26,5,6] in terms of the average number of volume truncations, with a challenging volume fractions testing sequence and without relying on a relatively complex geometric parameterization of the truncated volume.Its relative simplicity and the usage of geometrical data that are already available in the geometrical VOF method simplifies the adoption of the proposed CCS algorithm in existing numerical codes.

Figure 1 :
Figure 1: Multiphase domain and its discretization with the Volume-of-Fluid method.

Figure 3 :
Figure 3: Geometry and the volume fraction α c (s) of the endo-dodecahedron.
k is the truncated volume given by the plane normal n c and the vertex i p associated with the k-th signed distance in the list.The algorithm starts with the central index k c = IN T [(I p + 1)/2] in the signed distance list, and computes the truncated volume V T passing through this point.If |V T | > α c |V c |, k max = k c , otherwise k min = k c .The next iteration continues with the reduced list of indices between [k min , k max ].The authors state that the algorithm complexity in terms of the CPU time is O(

Figure 4 :
Figure 4: Cap area A(s) used as the exact volume derivative V (s).

Figure 8 :
Figure 8: Iteration distributions of the NCS and CCS algorithms for test polyhedrons.

Figure 9 :
Figure 9: CPU time distribution between the sub-algorithms of the NCS and CCS algorithm.
a )| < or |α (s b )| < then Non-zero derivative at s a or s b .
already show that those algorithms are significantly outperformed by the NCS algorithm.The reduction of iterations achieved by the CCS algorithm is reflected in fig.8in the iteration distribution for each test polyhedron, when compared to the NCS

Table 1 :
Computer architecture used for testing.algorithm.Only 2 iterations are required on average to position the interface within test polyhedrons.Results shown in figs.7a, 7b and 8 confirm that the CCS algorithm significantly improves the computational efficiency of the interface positioning algorithm on polyhedrons, including star-shaped non-convex polyhedra with non-planar faces.