Three-dimensional cellwise conservative unsplit geometric VOF schemes

This work presents two unsplit geometric VOF schemes that extend the two-dimensional cellwise conservative unsplit (CCU) scheme (Comminal et al. (2015) [49]) to three dimensions. The novelty of the 3D-CCU schemes lies in the representation of the streaksurfaces of donating regions by polyhedral surfaces whose vertices are calculated with the 4th order Runge-Kutta scheme. Moreover, the advected liquid volumes are computed using a truncation algorithm (López et al. (2019) [62]) suited for arbitrary non-convex and self-intersecting polyhedra, which removes the need for tetrahedral decomposition. The 3D-CCU advection schemes were coupled to three interface reconstruction methods (Youngs’ method, the Mixed Youngs-Centered scheme, and the Least-Square Fit algorithm). The resulting VOF methods were tested in classical benchmark advection tests, including translation, rigid-body rotation, shear and deformation ﬂows. The proposed 3D-CCU schemes conserve the liquid volume and maintain the physical boundedness of liquid volume fractions to the machine precision. The 3D-CCU schemes perform favorably compared to other unsplit geometric VOF schemes when coupled to Youngs’ interface reconstruction method. Moreover, the 3D-CCU schemes coupled to the Least-Square Fit algorithm are more accurate than most other VOF schemes that use a second-order accurate interface reconstruction, except those where a 3D extension of the Mosso- Swartz interface reconstruction is employed. The comparison of the different VOF schemes highlights the importance of coupling accurate interface reconstruction methods with accurate unsplit advection schemes.


Introduction
The Volume-of-Fluid (VOF) method is a widely used interface capturing technique to solve free-surface and multiphase flows of incompressible immiscible fluids. The VOF method was originally developed by Hirt and Nichols [1] to supersede the earlier Marker-and-Cell (MAC) method [2]. Unlike interface tracking methods (for example the MAC [3] and fronttracking [4,5] methods), the main idea of the VOF method consists in tracking the volume of liquid phases inside mesh cells instead of the interface location. Thereby, the VOF method is inherently connected to the volume conservation of incompressible fluids. Moreover, the Eulerian representation of liquid phases has the advantage of naturally and automatically handling topological changes of the interface, i.e. interface breakup and merging, which makes the VOF method suitable for complex two-phase flows. Other Eulerian interface capturing methods include the level-set method [6,7], the phase-field method [8], and coupled methods such as the coupled level-set/VOF [9][10][11] and phase-field/VOF [12] methods. Hybrid Eulerian-Lagrangian methods have also been developed, e.g. the mixed marker VOF algorithm [13], the hybrid particle levelset [14] or the level-set/front-tracking [15] methods, to cite a few. The advantages and limitations of the main interface capturing approaches have been discusses in [16].
Overall, state-of-the-art implementations of the VOF method are considered more accurate and versatile than the other techniques, as they show good to excellent conservation of mass, momentum and kinetics energy [16]. As such, the developments in the VOF method have enabled numerical studies of complex two-phase flow phenomena, such as the atomization of a liquid jet [17][18][19][20][21], the dynamics of raising bubbles [22,23], droplet deformation [24][25][26] and cavitation [27,28]. Those examples are particularly challenging to simulate as they involve flow features on multiple scales, interface fragmentation, strong surface tension effects, and large density and viscosity ratios of the fluid phases. The VOF method has also been used in a countless number of other applications spanning a diversity of scientific and engineering fields, since it had been implemented in open source codes as well as commercial Computational Fluid Dynamics software (for example FLOW-3D ® [29]).
The VOF method can be described as follows. The volume of fluid of a liquid phase at time t inside a cell C is defined as the volume integral where x is the spatial position and f (x, t) is the indicator function representing the distribution of the liquid phase as follows f (x, t) = 1 if the liquid phase is present at (x, t), The indicator function f is continuously defined inside the liquid phases but is discontinuous at the interface. The volume fraction of liquid phase inside the cell C is defined as where V C =´C dV is the volume of the cell. Assuming incompressibility of the liquid volume, the indicator function is a conserved passive scalar quantity that is transported with the flow. Therefore, it obeys the advection equation where u is the velocity vector field. The advection Eq. (4) is coupled to the equation of motion of the flow (e.g. the Navier-Stokes equations) governing the velocity vector field. An evolution equation of the liquid volume fraction can be derived by integrating Eq. (4) over a fixed cell C , during the finite time interval [t, t + t]. A semi-discrete evolution equation of the volume fraction is obtained upon application of the divergence theorem, yielding where ∂C is the boundary of C and n ∂C is its local unit normal vector pointing outward the cell. The integral term in Eq.
where S k and n k are the surface and the outward-pointing unit normal vector of the k face, respectively. Different flavors of VOF methods exist, depending on how the right-hand side of Eq. (5) is calculated. The VOF methods can be classified into two broad families: the algebraic and geometric schemes.
The algebraic VOF schemes [30][31][32][33][34][35][36] directly interpolate the indicator function f at the cell's faces from the values of the liquid volume fraction inside the cells, using high order differencing schemes. However, as the exact indicator function is discontinuous, algebraic schemes suffer from numerical diffusion, which smears the representation of sharp interfaces. Algebraic VOF schemes are generally not discretely conservative. Nevertheless, algebraic VOF schemes still remain a popular choice because they are relatively easy to implement into an existing finite-volume solver [16].
In contrast, the geometric VOF schemes approximate the indicator function f through a geometric reconstruction of the interface. The standard reconstruction techniques involve Piecewise Linear Interface Calculations (PLIC) where the interface location inside each cell is explicitly defined by a plane equation, which maintains the interface sharpness. In addition, geometric VOF schemes give geometric interpretations of the integral terms in Eq. (5). The flux of liquid volume through a face corresponds to the volume of fluid inside the donating region of that face [37]. In other words, the donating region contains the set of all Lagrangian particles that will pass through the face during the time interval [t, t + t] and contribute to the volumetric flux (we refer to [37] for a rigorous mathematical definition of donating regions). Note that the donating region has a signed volume with respect to the face's orientation defined by the normal vector n k in Eq. (6). Geometric VOF schemes use constructive geometrical procedures to advect the volume of fluid inside the cells. However, the complexity of combinatorial computational geometry is often viewed as an obstacle to the implementation of geometric VOF schemes in 3D [16].
Directionally-split schemes [38][39][40] reduce the complexity of the geometrical operations by decomposing the multidimensional advection into a succession of one-dimensional advection steps. However, the directional splitting introduces additional approximations that produce numerical diffusion of the liquid volume fraction in rotational flows. In addition, directionally-split schemes are only suited for Cartesian grids and cannot be extended to unstructured meshes, unlike unsplit schemes.
The development of accurate, robust and efficient unsplit 3D geometric VOF schemes is an active field of research. A thorough critical review of contemporary unsplit geometric VOF schemes, as well as their historical developments was recently conducted by Marić et al. [41]. Unsplit geometric VOF schemes generally adopt one of the two approaches: either the Eulerian fluxwise advection, or the semi-Lagrangian cellwise advection.
In the Eulerian fluxwise approach, the donating regions are computed by tracing the timelines and the streaktubes of the faces; see Fig. 1 (timeline in blue, streaktube in red). Thereafter, the liquid volume flux is calculated as the intersection of the donating regions with the reconstructed PLIC (deep blue area in Fig. 1). The liquid volume flux may either be positive ( Fig. 1(a)) or negative ( Fig. 1(b)) depending on its orientation with respect to the outward-pointing normal vector of the face. The donating region may also intersect with the cell's face, as shown in Fig. 1(c). In that case, one part of the donating region has a positive contribution to the liquid volume flux and the other part has a negative contribution.
In contrast, the semi-Lagrangian cellwise approach tracks the displacement of the entire volume of the grid cells. Lagrangian tracking defines the images of the grid cells through the flow map [42]. The volume of fluid is subsequently updated by remapping the liquid volume of the cell's image onto underlying grid cells (or vice versa). Despite the difference of point of view, Zhang [42] proved that both approaches are mathematically equivalent, if the calculation of the donating regions in the fluxwise approach is consistent with the Lagrangian flow map used in the cellwise approach.
Despite the presumed superiority of the unsplit geometric VOF schemes in term of accuracy, only few of them achieve discrete volume conservation and physical boundedness of liquid volume fractions (0 ≤ α C ≤ 1). Rider and Kothe [43] laid down a list of requirements that flux-based unsplit geometric VOF schemes must fulfill in order to be conservative. Specifically, the computed donating regions of adjacent faces must not overlap or form gaps. Moreover, the incompressibility of the flow requires a null sum of the signed volume of all the donating regions of a cell's faces. López et al. (2D EMFPA scheme) [44] showed that the calculation of donating regions is more accurate when streaktubes are computed using cell's vertex velocities instead of face-centered velocities. The 2D EMFPA scheme demonstrated a substantially higher accuracy than earlier unsplit schemes [43,[45][46][47]. Cervones et al. [48] (GPCA scheme) and Comminal et al. [49] (CCU scheme) built upon the 2D EMFPA scheme and showed that the shared edge of adjacent faces must also have the same timeline in their respective donating regions; see Fig. 1(d). This condition enforces temporal consistency of the scheme and it is a requirement for the physical boundedness of liquid volume fractions [50]. Otherwise, a cell could donate some liquid volume before receiving it, which may result in unphysical volume fractions (α C outside [0, 1]), as illustrated in Fig. 1(e). Moreover, the geometric VOF algorithm must adequately handle non-convex and self-interesting donating regions, which frequently occur in vortical flows. Both the two-dimensional GPCA and CCU schemes fulfill all the aforementioned requirements and achieve discrete volume conservation and physical-boundedness. The first 3D unsplit geometric scheme using cell's vertex velocities was proposed by Liovic et al. [51] (PCFSC scheme). It was realized that linear approximations of streaklines arising from the cell's vertices define non-planar streaksurfaces of 3D donating regions. The PCFSC scheme modifies the streaklines of cell's vertices in a way that maintains planarity of streaksurfaces. Hernández et al. [52] proposed another scheme (FMFPA-3D) where adjacent polyhedral donating regions have matched faces. Both PCFSC and FMFPA-3D schemes adopt the Eulerian fluxwise approach and represent donating regions by convex 8-vertex polyhedra, which facilitates geometrical operations (i.e. polyhedral truncation and volume calculation). However, the departure from the actual strealklines of the flow map alters temporal consistency and still produces small overlaps of neighboring donating regions, which degrades volume conservation and boundedness.
Le Chenadec and Pitsch [53] proposed a 3D semi-Lagrangian cellwise VOF scheme coupled with the level-set method, where images of Cartesian cells are represented by 8-vertex polyhedra decomposed into 6 tetrahedra. The tetrahedral decomposition effectively triangulates non-planar surfaces, handles non-convex polyhedra, and reduces the complexity of geometrical operations. Moreover, the cellwise approach automatically satisfies temporal consistency of underlying donating regions. The geometric scheme of Le Chenadec and Pitsch [53] does not comprise a volume correction procedure, but volume conservation errors were reduced by combining forward and backward temporal integrations forming a trapezoidal scheme. A 3D conservative and bounded cellwise unsplit geometric scheme was designed by Owkes and Desjardins [54]. Volume conservation was achieved through a volume correction of the cell's image, represented by a 14-vertex polyhedron (with 24 triangular faces), decomposed into 20 tetrahedra. This tetrahedral decomposition properly handles non-convex cell's images and self-intersecting donating regions, however, it only applies to Cartesian grid cells. Jofre et al. [55] generalized the VOF method of Owkes and Desjardins to unstructured meshes and derived a temporally-consistent fluxwise advection scheme. Marić et al. [56,57] (voFoam and UFVFC schemes) proposed several enhancements to the schemes of Owkes and Desjardins [54] and Jofre et al. [55], suited to unstructured meshes. In particular, Marić et al. [57] implemented a "flux-aware" triangulation of non-convex ruled surfaces, which removes potential overlapping issues of the tetrahedral decomposition.
Ivey and Moin [58,59] proposed several conservative fluxwise unsplit geometric schemes (NIFPA schemes) on unstructured grids, where donating regions are computed sequentially, in a way that avoids overlapping of adjacent donating regions. The NIFPA schemes also decompose polyhedral donating regions into tetrahedral volumes. The NIFPA schemes achieve physical boundedness of the liquid volume fraction, despite not fulfilling temporal consistency of donating regions. Moreover, the sequence of calculation of the donating regions was randomized to avoid any directional bias of advection, which renders the NIFPA schemes non-deterministic.
Rønby et al. [60] proposed the isoAdvector conservative unsplit VOF scheme for unstructured meshes. Unlike other geometric schemes, the isoAdvector scheme reconstructs the interface with an iso-surface of liquid volume fraction, instead of a PLIC. The reconstructed iso-surfaces are not required to be planar. Moreover, the isoAdvector scheme computes the integral of the liquid volume flux in Eq. (6) by tracking the evolution of the immersed area fraction of the faces into the liquid volume, when advecting the interface. In some aspects, the isoAdvector scheme can be viewed as a fluxwise advection scheme with a forward temporal integration of the advected volume instead of the backward integration of donating regions.
To date, most unsplit geometric VOF schemes use a tetrahedral decomposition of polyhedral volumes, in order to reduce the complexity of geometric operations (i.e. polyhedral truncation and volume calculation). However, this approach is not the most efficient, and it does not extend to arbitrary polyhedra. It is generally desirable to have flexible routines that can operate on arbitrary polyhedra, using robust and efficient algorithms, with a minimal number of special cases to handle. López et al. [61][62][63] developed a versatile toolbox of analytical and geometrical polyhedral operations for arbitrary non-convex polyhedra. The present work applies the algorithms of the geometric toolbox of López et al. [61][62][63] to construct two 3D semi-Lagrangian cellwise conservative unsplit (3D-CCU) geometric VOF schemes, based on a backward and forward time integration. The proposed 3D-CCU schemes adopt the cellwise approach, which is temporally consistent. Both schemes use a 4th order Runge-Kutta integration to advect the position of the cell's vertices forming polyhedral donating regions. In addition, the polyhedral donating regions are corrected with pyramidal volume corrections to maintain volume conservation. Thus, the proposed 3D-CCU schemes achieve volume conservation and physical boundedness to the machine precision. Moreover, using the robust geometric truncation algorithm of López et al. [62] for arbitrary non-convex polyhedra avoids the need of tetrahedral decomposition, which greatly simplifies the implementation of the geometric VOF schemes. Polygonal faces do not need to be triangulated either, and the issue raised in [57] about degenerated decomposition of non-convex polyhedra is avoided. The 3D-CCU advection schemes have been tested with first-and second-order accurate interface reconstruction schemes, including Youngs' method [64], the Mixed Youngs-Centered (MYC) scheme [38], and a modified version of the Least-Squares Fit (LSF) algorithm [38,65]. The backward 3D-CCU scheme is an extension of the two-dimensional CCU scheme [49].

Forward and backward semi-Lagrangian tracking
This section describes the two 3D-CCU schemes, based on forward and backward semi-Lagrangian tracking. Forward and backward semi-Lagrangian advection schemes were previously discussed in [41,51]. Below is a brief summary of the two schemes. In

Enforcement of volume conservation and boundedness
Without loss of generality, polyhedral grid cells are defined by the position of their vertices and the face-vertex connectivity of the mesh. The face-vertex connectivity is invariant through the Lagrangian tracking. Thus, the calculation of cells' images is primarily concerned with computing the images of the cells' vertices. However, when a cell has faces with more than 3 vertices (e.g. in hexahedral mesh), the images of those vertices are generally not co-planar. Therefore, the images of the cell's faces have ruled surfaces. A polyhedral representation of the cell's image is obtained by triangulating each face's image using an additional face-centered vertex. As an example, the image of a Cartesian grid cell (6 faces, 8 vertices) is represented by a 14-vertex polyhedron of 24 triangular faces; see Fig. 3.
As mentioned in the introduction, the volume conservation of incompressible fluids requires that the cells' images have the same volume as their original grid cells. The proposed 3D-CCU schemes enforce this requirement by adjusting the position of the additional face-centered vertices, similarly as the pyramidal correction described in [54,55,57]. This procedure is referred to as the volume flux correction.
In virtue of equivalence between the cellwise and fluxwise approaches, the volume of the cell's image is conserved if and only if the discrete volume flux through the cell's faces is divergence-free [42]. Assuming that the discrete velocity vectors fulfill the continuity equation of incompressible fluids, the following relation between face-centered velocities holds where A k and u k are the area and the discrete face-centered velocity of the k face, respectively. The product A k u k ·n k is the net volume flux through the k face. Its integration over time must correspond to the signed volume V DR,k of the donating region. Therefore, the discrete volume fluxes of a cell are made divergence-free by enforcing that the sign volume of each where the sign in the upper bound of the integral depends whether it is a forward or backward temporal integration. The constraint (8) is used to calculate the adjustment in the position of the additional face-centered vertex in the image of each cell's face.
Volume conservation of geometric VOF schemes also requires that the interface reconstruction inside the cells enforces the correct volume of fluid. Note that the forward advection scheme requires the interface reconstruction inside the cell's image at time t + t, whereas the backward scheme needs the interface reconstruction inside the grid cells at time t − t.
In principle, the forward scheme could also advect phase-specific volumes without tracking the entire cell, as described in [41]. However, the conservativeness and physical-boundedness of the advection schemes is easier to achieve by applying the volume flux correction to the entire cell. In addition, the cellwise advection approach ensures temporal consistency of the CCU schemes and significantly reduces the possibility of overlaps between cell's image.

Implementation details
This section provides the implementation details of the two 3D-CCU VOF schemes, for a Cartesian grid mesh. The proposed advection schemes could be extended to unstructured mesh with convex cells of arbitrary shapes. However, the development of VOF-PLIC algorithms suited to arbitrary polyhedral cells is outside the scope the present work.

Lagrangian marker tracing
The image P 1 of a cell's vertex P 0 is calculated by tracing the position of a Lagrangian marker seeded in position P 0 at time t and transported with the flow during t. The temporal integration of the marker's position uses the classic 4th order Runge-Kutta scheme where U 0,P 1  velocities. We assume that discrete velocities are known at times t and t + t, in the forward scheme, and at times t and t − t, in the backward scheme. However, if the VOF scheme was coupled to a flow solver where liquid volume fractions are updated before the velocity field, the velocities would be unknown at t + t in the forward scheme, and t in the backward scheme. In that case, the intermediate velocity estimates U 0,P 1 , U 1,P 1 , U 2,P 1 , U 3,P 1 would all be evaluated at time t in the forward scheme, and t − t in the backward scheme.
For higher accuracy, the streakline from P 0 to P 1 is approximated by a 2-segment polygonal line P 0 -P 1/2 -P 1 , where the immediate points P 1/2 are calculated by tracking an intermediate marker seeded in P 0 at time t + ξ t 2 , during the

Streaktube calculation
The streaklines of the cell's vertex are used to calculate the streaktube of the cell's faces. Specifically, a streaksurface is computed for each edge of a cell's face. This streaksurface is bounded by the streaklines of the 2 vertices forming that edge. The positions of two additional points, E 1/3 and E 2/3 , are also calculated by tracing the positions of Lagrangian markers seeded at the middle of the edge at the times t + ξ t respectively. The streaksurface of the edge is patched by 8 triangular facets, according to the pattern shown in Fig. 4(b). Then, the streaktube of the cell's face is obtained by combining the streaksurfaces of all it edges; see Fig. 4(c). Different patterns were tested to triangulate the streaksurfaces, including some with more triangular facets and additional intermediates points along the streaklines of the cell's vertices and the edges' middle points. However, increasing the number of triangular patches beyond 8 facets was deemed unnecessary, as it did not improve the accuracy of the advection schemes in the benchmark tests presented in Section 4.

Volume flux correction
In order to compute the donating regions, the streaktube needs to be capped by the timeline of the face at the end of the time integration. The timeline of the face being a ruled surface, it also needs to be triangulated. The triangulation uses the barycenter F * 1 of the images P 1 , Q 1 , R 1 , S 1 of the face's vertices; see Fig. 4(d).
At this point, the donating region is represented by a polyhedron. For instance, the donating region of a rectangular face is a polyhedron made of 37 faces (36 triangular facets plus the rectangular cell's face). The signed volume V DR,k of the polyhedron is unambiguously defined and can be calculated using a general formula for arbitrary non-convex polyhedra; see Section 3.5. However, the signed volume of the polyhedral donating region is generally different from the net volume flux The volume of the polyhedral donating region is corrected accordingly by moving the position of the vertex F * 1 . Specifically, the volume of the polyhedral donating region is adjusted by moving the position of F * 1 along the direction of the average unit normal vector n F of the facets {P 1 Q 1 F * where λ is a scalar. Depending on the sign of λ, the pyramidal volume flux correction P 1 Q 1 R 1 S 1 F 1 F * 1 adds or removes the volume V cor,k onto the capping faces of the polyhedral donating regions; see Fig. 4(d). Knowing the required volume correction V cor,k , the unit normal vector n F , the position of F * 1 as well as the vertices {P 1 , Q 1 , R 1 , S 1 } of the face's image, the correct value of λ can easily be calculated analytically, as described in [54,55,57].
Once all the polyhedral donating regions have been computed and corrected, a polyhedral image of the grid cell is obtained by combining the image of all its faces. Therefore, the image of a rectilinear grid cell is a 14-vertex polyhedron comprising 24 triangular faces, as shown in Fig. 3(c).

Interface reconstruction in grid cells
Depending on the advection scheme, the interface reconstruction is either required inside the grid cell (backward scheme) or inside the cell's image (forward scheme). In practice however, it is easier to reconstruct the interface inside a Cartesian cell than in an arbitrary 14-vertex polyhedron representing the cell's image. Thus, the PLIC is always reconstructed inside the grid cells, but for the forward scheme, it is additionally mapped into the cell's image. Three alternative interface reconstruction algorithms have been used to determine the PLIC normal vectors n 0 in interfacial cells. Once the unit vector normal is known, the correct positioning of the PLIC inside Cartesian cells is calculated exactly, using the analytical expressions derived in [66].

Youngs' method
Youngs' method [64] approximates the normal vector as the gradient of the liquid volume fraction where · 2 refers to the Euclidean vector norm. The gradient of the liquid volume fraction is first calculated at the vertices of the Cartesian grid, using finite differences. The PLIC normal vectors inside interfacial cells are obtained by averaging the liquid volume fraction gradients from the 8 vertices of the cells, similarly as in [38,67]. Youngs' method provides a cheap estimate of the PLIC normal vectors, with a first order accuracy [38,47].

Mixed Youngs-Centered method
A slightly more accurate estimate of the PLIC normal vectors is obtained using the MYC method [38]. For each interfacial cell, the MYC method computes four candidates of the normal vector. One candidate is obtained using Youngs' method, and three other candidates are evaluated by the centered column differences in three directions, inside the 3 × 3 × 3 stencil of neighboring cells. The MYC algorithm selects the unit normal vector from the centered column differences that has the largest component (in absolute value) along its columns' summation direction, unless that maximum component is larger than the component obtained with Youngs' method. In that case, the normal vector from Youngs' method is chosen. The rationale behind these choices is explained in [38]. Both Youngs' method and the centered column differences schemes are formally first order accurate. Nevertheless, in practice the MYC method can often reach apparent convergence rates with respect to mesh refinements between first and second order accuracy, thanks to the careful selection of the best normal vector candidate. A known shortcoming of the MYC algorithm (as well as Youngs' method) is that it does not reconstruct exactly all 3D planes, which can increase the discontinuity gap between the PLIC of adjacent cells.

Least-Square Fit method
A second-order accurate interface reconstruction is obtained with the LSF method [38,65]. The LSF algorithm requires an initial approximation of the PLIC normal vectors (which is obtained with the MYC method in this work). Then, the LSF method produces an improved estimate of the normal vectors using information from the other PLIC approximations in the 3 × 3 × 3 stencil of neighboring cells. The LSF method tries to reduce discontinuities of the PLIC across neighboring cells. 8 For each interfacial cell k in the 3 × 3 × 3 stencil, the geometric center X k of the PLIC polygon is computed. The LSF estimate of the normal vector is determined by searching for the best plane that minimizes the least-square distance from the X k points. Note that the fitting plane is constrained to pass through the geometric center X 0 of the initial PLIC polygon in the central cell of the 3 × 3 × 3 stencil. Thus, the LSF method minimizes the function where N k is the number of neighboring interfacial cells, and m 0 is the improved estimate of the normal vector inside the central cell. Realizing that H(m 0 ) is a quadratic function of the three components of m 0 , it is minimized by solving the linear system of equations given by However, as (16) is a dependent system (the 3D plane orientation is determined by two independent variables), one of the components of m 0 has to be set to an arbitrary non-zero value, such that (16) reduces to an independent system of two equations. This is done by setting the component that has the largest absolute value in the initial approximation of the normal vector to unity. The two other components of m 0 are determined by solving the system (16). Finally, the LSF unit normal vector is obtained as n LSF 0 = m 0 / m 0 2 . The LSF procedure is applied only once, as prescribed in [38].
Note that the LSF method used in this work contains some modifications compared to the original method presented in [38]. Firstly, all X k points have the same weight in the function H . Secondly, the function H only includes the points X k from neighboring interfacial cells whose initial normal vectors form an angle lower than 45 • with the initial normal vector approximation of the central cell. Adjacent interfacial cells with normal vectors forming a larger angle than 45 • indicate that the interface is underresolved in that region, in which case the LSF method is unlikely to improve the normal approximation from the MYC method. This discrimination strategy was first proposed in [68], in the context of Swartz interface reconstruction method. Thirdly, any interfacial cell that does not have at least one full and one empty cell among its adjacent neighbors is deemed unsuited for the LSF, and its normal vector is computed with the MYC method. This proposed exclusion criterion was found to reduce advection errors due to the PLIC in underresolved regions with a radius of curvature about the grid size or lower.

Interface remapping into cells' images
In the forward CCU advection scheme, the PLIC of the grid cells is remapped into the cells' images. The intersection of the PLIC with a Cartesian grid cell produces a 3D polygon that has between 3 to 6 vertices. The PLIC polygon of the grid cell is mapped into the cell's image by isomorphism; see Fig. 5. Let I 0 be the intersection of the PLIC with an edge P 0 -Q 0 of the grid cell. The images of P 0 and Q 0 by Lagrangian tracking are noted P 1 and Q 1 , respectively. Then, the isomorphic mapping of I 0 into the cell's image gives the point where x = I 0 P 0 2 / Q 0 P 0 2 . The mapping of the PLIC inside the cell's image is generally non-planar, if the PLIC polygon has more than 3 vertices. In that case, the average normal vector n 1 of the non-planar polygon is computed using Newell's method [69]. Once an estimate of the normal vector is known, the PLIC of the cell's image needs to be repositioned such that it enforces the correct volume of fluid (that is the same as in original grid cell). The vector position x of all the points belonging to the PLIC of the cell's image is parametrized by the following equation where d is the signed distance from the origin to the plane.
The volume of fluid below the PLIC can be evaluated by truncating the polyhedron of the cell's image by the PLIC plane and computing the volume of the truncated polyhedron. In 3D, the function V PLIC,n 1 (d) that returns the volume of fluid below the PLIC as a function of d, for a set normal vector n 1 , is a continuous piecewise third order polynomial. Thus, each piece of the function V PLIC,n 1 (d) is defined by four coefficients c i , i = {0, 1, 2, 3}, where The piecewise polynomial function changes coefficients every time a vertex of the polyhedron changes side with respect to the plane, when d is increased. Moreover, V PLIC,n 1 (d) is strictly monotonous. The correct value of d enforcing the required volume of fluid υ C is computed by solving the cubic equation However, we first have to determine the values of the four coefficients c i . The first step consists in finding the two vertices of the polyhedron that bracket the correct interval containing the solution of Eq. (20). Those two bracketing vertices can be obtained after computing the volume of the polyhedron truncated by each plane passing through one of the vertices. The lower and upper bracketing vertices correspond to the vertices that return the larger and lowest truncated volume below and above υ C , respectively. However, computing the polyhedral truncations is relatively expensive. Therefore, it is desired to use a sequential procedure that only requires a minimum number of polyhedral truncations. For that, the vertices of the polyhedron are sorted with respect to their signed distances to an initial guess of the PLIC position, which is a very cheap computational task compared to the polyhedral truncation. The plane with normal vector n 1 , passing through the center of mass of the vertices of the isomorphic image of the PLIC polygon is taken as initial guess. Then, the lower and upper brackets are iteratively searched through the sorted list of vertices. The proposed bracketing procedure requires a small number of polyhedral truncations, if the initial guess is sufficiently close to the correct brackets. Other searching strategies of the brackets are discussed in [70].
Once the brackets are found, two additional polyhedral truncations with planes positioned at 1/3 and 2/3 of the bracketed interval are computed. Those 2 polyhedral truncations, plus the polyhedral truncations of the upper and lower bracketing vertices provide 4 points of the third order polynomial V PLIC,n 1 (d). Thereby, the 4 coefficients c i in Eq. (19) are determined by curve fitting, similarly as in [60]. Once the four coefficients c i are known, the roots of the cubic Eq. (20) are calculated analytically, using Cardano's formula. The cubic Eq. (20) either has 1 or 3 real roots, but only one of the roots lays inside the bracketed interval, as V PLIC,n 1 (d) is strictly monotonous. Note that the proposed PLIC remapping does not require the cell's image to be convex.

Polyhedral operations
The 3D-CCU schemes heavily rely on two polyhedral operations: (i) a truncation algorithm returning the cut of a polyhedron by an infinite plane, and (ii) a function that computes the volume of a given polyhedron. Those two polyhedral operations are used to compute the volumes of fluid inside a cell and a cells' image, by cutting its polyhedron with the PLIC plane. Moreover, since the grid cells are convex, the intersection of a polyhedral cell's image with the underlying grid cells can be computed by successive truncations of the cell's image with all the planes of the grid cell's faces.
In addition, those two subroutines require a flexible data structure capable of holding non-convex polyhedra with arbitrary numbers of faces and vertices (where faces also have an arbitrary number of vertices). Note that the polyhedral truncation may return a polyhedron with holes and disconnected regions, which must be correctly represented into the data structure. In vortical flows, the donating region polyhedron may also be self-intersecting with the cell's face. In that case, the volume on one side of the face has a positive contribution to the net volume flux, whereas the volume on the opposite side has a negative contribution.
The polyhedral truncation subroutines and the polyhedral data structure of [62] are adopted. The polyhedral data structure essentially uses two arrays. First, a vertex table holds the 3D coordinates of every vertex of the polyhedron. Secondly, a face table contains an ordered list of vertices, for each face of the polyhedron. The orientation of the faces is implicitly defined by the ordered vertex list, according to an orientation convention (e.g. vertices traveled clockwise/counterclockwise when looking from the inside/outside of the face). A hole inside a face is represented by an additional polygonal face with the reversed orientation. The overlapping areas of the two faces having opposite orientations cancel each other and produce a logical hole into the polygonal face. The orientation of the faces also determines the sign of the polyhedral volume. Polyhedral regions bounded by faces with an inverted orientation (faces turned inside out) have negative volumes. Specifically, the sign volume V POLY of a polyhedron is calculated using the following general formula derived from the divergence theorem [71] where n k is the outward-pointing unit normal vector of the k face, A k is its area, and P k,0 is the position vector of an arbitrary point of the k face (e.g. one of its vertices). Similarly, the area A k of the k polygonal face is computed with the formulas [72] A k = 1 2 n k · j P k, j × P k, j+1 , (22) where P k,0 , . . . , P k,N−1 are the position vectors of the N vertices of the k face, traveled according to the ordered face list.
Having a face orientation is an essential feature for computing the correct sign volume of polyhedral donating regions with a negative volume flux (including self-intersecting polyhedra) using Eq. (21).
A polyhedral truncation is computed by cutting each polygonal face by the truncating plane. Each time an intersection is detected, the intersection point between the plane and the edges of the face is added into the vertex table, and the ordered vertex list of the face is updated accordingly. Then, the capping faces of the cut polyhedron are formed by traveling through adjacently truncated faces and joining the successive intersection points. When a polyhedral truncation produces disconnected regions, the algorithm automatically creates several overlapping capping faces with opposite orientations connecting the disjoined regions, such that their suppositions yield the correct area. The detailed procedure that creates topologically consistent capping faces is described in [62].
The geometric algorithms described in [62] provide robust and efficient subroutines to perform the polyhedral operations. Unlike other VOF schemes described in [51,[53][54][55][56][57], the current method does not require the decomposition of polyhedra into tetrahedra. As discussed in [57], the algorithms relying on the tetrahedral decomposition do not always return the correct signed volume, if non-convex faces are not correctly triangulated. The proposed method is immune of such problems and naturally handle non-convex and self-intersecting polyhedra without need of any specific procedure.

Benchmark numerical tests
The 3D-CCU advection schemes (backward and forward) were coupled with Youngs', the MYC and LSF interface reconstruction methods. The VOF schemes were implemented on a structured Cartesian mesh for the purpose of validating the proposed methods. Nevertheless, the 3D-CCU schemes could as well be implemented on unstructured meshes of convex cells using the two geometric routines (volume calculation of arbitrary polyhedra and truncation by an infinite plane) discussed in the previous section.
The proposed 3D-CCU VOF schemes were applied to benchmark surface capturing tests with prescribed velocity fields in translation, rotation, shear, and deformation flows. In all benchmark tests, the liquid volume was initialized as a sphere. The liquid volume fraction of the sphere inside each Cartesian cell was calculated using the exact formulas derived in [73]. At the end of the tests, the volume conservation error is calculated as the absolute difference of total liquid volume between the initial time T i = 0 and final time T f of the simulations Moreover, since the exact shapes of liquid volume at the end of the benchmark tests are also spheres, the geometric error is monitored as where α exact C (T f ) is the exact liquid volume fraction distribution at the end of the tests, calculated using the formulas from [73]. Geometric errors appear from the combination of advection and interface reconstruction errors. Thus, geometric errors are unavoidable when a non-planar interface is represented by a PLIC.
All tests are performed inside the unit-length cubic domain [0, 1] 3 , for several successively refined Cartesian meshes, with N grid cells along each direction. The order of convergence of the geometric error with respect to successive mesh refinements is calculated as The backward and forward schemes returned similar numerical results, with at least three equal significant digits in all the tests (and up to six in some cases). Thus, the numerical results reported in the next subsections cover both the backward and forward 3D-CCU schemes.  Fig. 6. For CFL = 1, the translation tests produce geometric errors near the machine precision. In uniform translation flows, the 3D-CCU schemes calculate cells' image exactly. Therefore, the convergence rate of the geometric errors with mesh refinements solely reflects the accuracy of the interface reconstruction method (and not the advection scheme).

Translation tests
Moreover, when CFL = 1, the cells' image coincide with diagonally adjacent cells, thus the 3D-CCU schemes advect the entire liquid volume from a cell to another and the interface reconstruction has no influence on the advection. This explains why the geometric errors for CFL = 1 are near the machine precision. On the finer meshes, the convergence rate in the geometric error approaches a second order accuracy, for the CCU-MYC and CCU-LSF schemes, whereas the CCU-Youngs scheme exhibits a sublinear convergence rate, for CFL = 0.25 and CFL = 0.5. The geometric errors also decrease with a larger time step size, as it reduces the number interface reconstructions (that introduces the reconstruction errors) performed during the test. On the coarser mesh, Youngs' method somehow produces a lower interface reconstruction error than the MYC and LSF methods, when CFL = 0.5. However, these numerical results with the coarser mesh cannot be used to draw a conclusion about the accuracy of the reconstruction methods because the sphere is underresolved (only 6 grid cells over the sphere diameter) and a fortunate interface reconstruction in a few cells can favorably reduce the overall geometrical errors of the test.

Rotation tests
The rotation tests consist in rotating a sphere of radius 0. In all rotation tests, the volume conservation errors have a magnitude order of 10 −16 . In simple rotation flows, Lagrangian tracing of markers with the 4th order Runge-Kutta method calculates the exact images of cell's vertices. The approximation of streaktubes (that are second order surfaces) by polygonal surface introduces small positioning errors of the face-centered vertices in the cells' images, but most of the advection errors come from the PLIC reconstruction. The geometric errors of the rotation tests are reported in Table 1, together with data of other unsplit geometric VOF schemes available in the literature [52,55,74]. The convergence order of the geometric error approaches the second order accuracy for the CCU-MYC and CCU-LSF schemes, while the CCU-Youngs scheme exhibits a first order accuracy. On the coarser mesh, the geometric errors of the CCU-MYC schemes are 11% larger than the geometric errors of CCU-LSF schemes. However, their relative differences reduce to 6% when N = 128.
On the finest mesh, the CCU-Youngs scheme outperforms the other schemes that use Youngs' interface reconstruction method. Moreover, the CCU-LSF scheme is more accurate than the FMFPA-3D scheme using the improved ELVIRA interface reconstruction [74] and Jofre's scheme using the LVIRA interface reconstruction [55] (with the Cartesian mesh), at CFL ≈ 0.5 and CFL ≈ 1.0, respectively. Note that the iterative LVIRA interface reconstruction and its "efficient" non-iterative variant, the  ELVIRA method, are also second-order accurate interface reconstruction methods, like the LSF algorithm. The FMFPA-3D [52] scheme is more accurate than the CCU-LSF scheme when coupled to the CLCIR interface reconstruction method [74]. The reconstructed PLIC of the rotated sphere at the end of the tests are represented in Fig. 7 for the LSF reconstruction method and CFL = 1. PLIC polygons in light blue are superimposed to the exact sphere position in dark blue. The interface reconstruction algorithm applies the LSF procedure to all interfacial cells, as the sphere is properly resolved.

Shear tests
In the shear tests, a sphere with radius 0.15 is initially centered at (0.5, 0.75, 0.25). The sphere is deformed through a non-uniform vortical shear flow with the following velocity components  Table 2, together with data of other geometric VOF schemes available in the literature [51,55,57]. On the finest mesh, the CCU-Youngs schemes have lower geometric errors than the directional-split, PCFSC and Jofre's schemes when they are also coupled with Youngs' interface reconstruction method. On the coarser mesh, the CCU-LSF scheme is more accurate than the PCFSC-CVTNA [51] and UFVFC-Swartz [57] schemes, whereas, the PCFSC-CVTNA and UFVFC-Swartz schemes perform better than the CCU-LSF on the finer mesh. It is also noticed that the CCU-LSF scheme has a degraded convergence order about 1.37 on the finer mesh, whereas the CCU-MYC maintains a quadratic convergence order (but with larger geometric errors than the CCU-LSF).
The reconstructed PLIC of the sheared sphere are represented in Fig. 8 at T f /2 and T f , for the tests using the LSF reconstruction method. The PLIC polygons oriented with the LSF procedure are colored in blue. Green and red polygons indicate the PLIC where the normal vectors were computed with centered column differences and Youngs' method, respectively. The reconstructed PLIC at the sharp "edges" of the sheared sphere were excluded from the LSF procedure because of underresolution (they did not pass the condition of having at least one full and one empty cell among their adjacent neighbors).

Deformation tests
In the deformation tests, a sphere with radius 0.15 is initially centered at (0.35, 0.35, 0.35). A non-uniform threedimensional vortical flow field with the following velocity components is prescribed as where T f = 3 is the final time of the test. The computational domain is discretized with N grid cells along each direction, where N = 32, 64, 128, 256. As in the shear tests, the velocity field is varying over time. After the half time T f /2, the flow is reserved, such that the deformed sphere is advected back to its initial shape, at the end of the test. The time step size is adjusted at each time step to maintain a maximum Courant number CFL = 0.5. The analytical formula to calculate the time step size t according to the prescribed CFL number in the accelerating flow is provided in Appendix. In all the deformation tests, the volume conservation errors have a magnitude order between 10 −15 and 10 −16 . The geometric errors of the deformation tests are reported in Table 3, together with data of other VOF schemes available in the literature [33,36,40,51,52,54,55,57,59,60,75,76]. On all mesh resolutions, the CCU-Youngs scheme outperforms the other split and unsplit geometric schemes using Youngs' interface reconstruction method. The CCU-Youngs scheme is also more accurate than the MULES (interFoam v1.5) [33] and UMTHINC/QHQ [36] algebraic VOF schemes, as well as the CLSVOF scheme [76] using a directional-split Eulerian advection. The CCU-MYC outperforms the directional-split Lagrangian VOF scheme [40] that uses the MYC interface reconstruction method. Both the CCU-MYC and CCU-LSF schemes are more accurate than the isoAdvector [60] scheme. In addition, the CCU-LSF scheme has lower geometric errors than the conservative unsplit VOF schemes of Owkes and Desjardins [54], Jofre et al. [55] and Ivey and Moin (NIFPA-1) [59] using the second order accurate ELVIRA, LVIRA and EHF [77] interface reconstruction methods, respectively. However, the PCFSC-CVTNA [51] and UFVFC-Swartz [57] schemes are more accurate than the CCU-LSF scheme, on the fine mesh resolutions for the PCFSC-CVTNA, and on all mesh resolutions for the UFVFC-Swartz. Interestingly, the CCU-MYC and CCU-LSF schemes are more accurate on the fine meshes than the directional-split Eulerian advection using the Moment-of-Fluid (MOF) interface reconstruction [76], whereas the MOF method [68,78] is presumably superior to the LSF method. This illustrates the importance of using accurate interface reconstruction in combination with accurate advection schemes.
The reconstructed PLIC of the deformed sphere are represented in Fig. 9 at T f /2 and T f , for the tests using the LSF reconstruction method. The PLIC polygons oriented with the LSF procedure are colored in blue. The green and red polygons indicate that the normal vector of the PLIC was computed with centered column differences and Youngs' method, respectively, because of underresolution. Notice how the deformation of the sphere produces a very thin liquid film that remains underresolved despite the successive mesh refinement (central region in red, reconstructed with Youngs' method). The superquadratic convergence rates seen on the finer meshes with the CCU-MYC and CCU-LSF schemes are attributed to the reduction of the central area reconstruction by Youngs' method with mesh refinements. On the coarser mesh resolution, the PLIC is unable to properly capture the thin liquid film, which results in an artificial interface fragmentation. This type of numerical errors has been named "numerical surface tension" and "numerical dispersion" in the literature [43,79]. It is initiated by orientation errors in the PLIC [43], which explains why Youngs' interface reconstruction method (which is only first-order accurate) is more prone to those errors, as illustrated in [41].

Discussion and conclusions
Two three-dimensional cellwise conservative unsplit geometric VOF schemes have been developed, based on backward and forward advection. The proposed 3D-CCU advection schemes use the semi-Lagrangian advection approach with a 4th order Runge-Kutta method to compute the images of grid cells. The donating regions of cells' faces and the cells' images are represented as polyhedral volumes. The 3D-CCU schemes conserve liquid volume to the machine precision, via pyramidal corrections of the polyhedral donating regions (provided that the discrete face-centered velocities are divergence-free). Moreover, the cellwise advection approach ensures the temporal consistency of the schemes, which maintains the physical boundedness of liquid volume fractions.
The implementation of the 3D-CCU schemes uses robust geometric algorithms and a polyhedral data structure that are suited for arbitrary non-convex and self-intersecting polyhedra [62]. This polyvalent polyhedral representation removes the need of the tetrahedral decomposition that was used in previous unsplit geometric VOF schemes and simplifies the computation of polyhedral truncations and polyhedral volumes. As such, the 3D-CCU schemes are applicable to unstructured meshes, although the current implementation is limited to Cartesian grids. The 3D-CCU advection schemes were coupled to Youngs' interface reconstruction method, the MYC scheme and the LSF algorithm.
The proposed VOF schemes were tested in benchmark advection cases with prescribed velocity fields, including translation, rigid-body rotation, shear and deformation flows. Both 3D-CCU schemes using the backward and forward advection approaches produced nearly identical results, with a least three equal significant digits, at an equivalent computation cost. Therefore, we recommend the choice between a backward or forward advection to be guided by other considerations, such as the convenience to integrate the VOF schemes into an existing flow solver.
In the tests with the finest meshes, the CCU-Youngs schemes outperformed all the other unsplit advection schemes that also used Youngs' method, as well as the algebraic VOF schemes. The CCU-LSF schemes compare favorably to other conservative unsplit VOF scheme using the second-order LVIRA, ELVIRA and EHF interface reconstruction methods. However, the PCFSC-CVTNA [51] and UFVFC-Swartz [57] schemes produce lower geometric errors than the CCU-LSF in the shear and deformation tests. Those two VOF methods use interface reconstructions based on 3D extensions of the Mosso-Swartz algorithm [80]. The Mosso-Swartz algorithms (and its variants) are probably more accurate than the LSF method, due to the iterative nature of their interface orientation corrections, whereas the LSF interface correction is applied only once. Moreover, the discrimination strategy adopted in our work, which excludes cells in underresolved regions from the LSF algorithm, also has some limitations. Numerical experiments show that applying Youngs' method or the centered column differences in underresolved regions with a high interface curvature (like the "sharp edge" of the liquid volume at the half time of the shear tests shown in Fig. 8) is beneficial. However, Youngs' method deteriorates the accuracy of interface reconstruction in underresolved regions with a thin liquid film but a low interface curvature (like the central area of the liquid volume at the half time of the deformation tests shown in Fig. 9). The issue could be avoided by using the MOF interface reconstruction method [68,78], which produces second-order accurate PLIC without using information from the neighboring cells. Moreover, the MOF method was shown to be consistently more accurate than the LVIRA and Swartz methods in [68,78]. Nevertheless, the favorable comparison of the CCU-MYC and CCU-LSF schemes to a directional-split VOF scheme using the MOF interface reconstruction illustrates that highly accurate interface reconstruction methods also need to be coupled with accurate advection schemes. Finally, the semi-Lagrangian cellwise nature of the proposed 3D-CCU schemes makes them suitable to advect other scalar, vectorial or geometric quantities, such as the center of mass of liquid volumes inside mesh cells used in the MOF reconstruction method.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.