A conservative cell-based unsplit Volume of Fluid advection scheme for three-dimensional atomization simulations

The Volume of Fluid (VoF) method is widely used for capturing the interface motion in multi-phase ﬂow simulations. In particular, unsplit geometrical advection schemes have proved well-suited for ﬂows with complex topologies. In the cell-based approach, the computational cell is followed along its Lagrangian trajectories and provides a conceptually simple framework for advancing the interface. Thus, this is a semi-Lagrangian method, and enforcing conservation is difﬁcult. This paper presents a cell-based three-dimensional (3D) unsplit advection scheme that is conservative. The method relies on the Hybrid Lagrangian Eulerian Method (HyLEM) of Le Chenadec and Pitsch [3] but additionally ensures discrete conservation by introducing a correction of the projected cell, which is inspired by the 3D ﬂux-based method of Owkes and Desjardins [14] and the two-dimensional cell-based method of Comminal et al. [4]. While the projected cell vertices are evaluated as in HyLEM, additional vertices are introduced to modify the projected cell faces. The positions of those are obtained from conservative ﬂux volumes. The proposed method provides the same accuracy as the method of Owkes and Desjardins [14], but is more efﬁcient. The proposed method is tested in various benchmark cases and applied in an atomization case.


Introduction
Atomization of liquids is known to play a key role in many engineering problems. Although great research efforts have been spent, many aspects of the complex multiscale atomization process are poorly understood. In the last decades, numerical simulations became a powerful tool to study atomization. However, numerical simulations are still challenging, and the development of numerical methods is ongoing research. An important part of these numerical simulations is the accurate prediction of the interface motion. For this, the Volume of Fluid (VoF) method is widely used, and unsplit geometrical advection schemes, which can be classified into cell-based and flux-based schemes [12], have proved well-suited for predicting the interface motion in atomization processes [9,2,16]. The drawback of those schemes is the computational cost, and Maric et al. [12] pointed out in their comprehensive review on advection schemes that the question of how geometrical schemes can be made more efficient needs to be addressed by future research. Three-dimensional (3D) unsplit flux-based advection schemes have been proposed, which are conservative, bounded, and second-order accurate [14,8,13]. Although the cell-based approach provides a conceptually simple and efficient framework for advancing the interface location, enforcing conservation is difficult and has not been yet realized in three dimensions. This paper proposes a 3D unsplit cell-based advection scheme that is conservative. The proposed scheme can be seen as an improvement of the method of Le Chenadec and Pitsch (LP) [3] with the addition of discrete conservation by introducing a correction based on conservative fluxes. In two dimensions, this is similar to the two-dimensional (2D) unsplit cell-based scheme of Comminal et al. [4]. While being discretely equivalent to the method of Owkes and Desjardins (OD) [14], the computational costs are significantly reduced. This paper starts with a description of the theoretical and methodological background before the proposed correction is presented. Afterward, the scheme is verified in a 2D and 3D transport test and applied in an atomization case.

Problem Formulation
VoF methods aim for capturing the interface evolution by using an indicator function that is defined by Here, x and t indicate the position vector and the time, respectively. If phase-change is neglected and both phases are considered incompressible, the indicator function evolution can be described in an integral sense by where Ω m (t) denotes a material volume and D m the material derivative.
In numerical simulations, it is solved for a volume-averaged value for each computational cell, which are called Eulerian Cells (EC) in this work. Generally, the volume-averaged value is referred to as the volume fraction and is defined for an EC Ω E k by Though it is solved for a volume-averaged quantity, it is important to note that a sharp representation of the interface is required during advection to obtain meaningful results. Two ways of addressing this challenge are the unsplit cell-based and the unsplit flux-based advection approaches. Unsplit cell-based volume fraction advection relies on solving the integral indicator transport equation over a static Eulerian mesh. For updating the volume fraction, phase-specific volumes or the entire ECs are either tracked forward or backward in time and remapped on the Eulerian mesh. The proposed method relies on the backward scheme of LP, which projects the whole EC backward in time. The volume fraction can be advanced according to where Ω L k indicates a Lagrangian Cell (LC), which is obtained by moving the EC backward in time from t n+1 until t n . As the EC is assumed to be a material volume, the motion is described by Here,u is the velocity vector and x is the position of a point on the surface of the material volume with x 0 as its initial position. Mass conservation is ensured if the backward projection of the EC conserves the volume, i.e V L k = V E k . In the flux-based approach, the change of the volume fraction in an EC is described in terms of the liquid phase embedded in flux volumes. These flux volumes are formed through streaktubes that are emitted from the EC faces. Mass is conserved if the flux velocities, defined by are divergence-free. Here, ∆t is the timestep, V F k,i is the flux volume, and A k,i is the area of the i-th EC face.

Methodology
For the proposed scheme, the general solution procedure for discretely solving Eq. (4), including the steps interface reconstruction, LC estimation, and remapping, can be adopted from the method of LP and are described in this section. The novelty of this work is an LC correction that is applied before the remapping step and is presented in detail in the following section. Interface reconstruction: As the method of LP, the presented advection scheme is based on the Piecewise Linear Interface Calculation (PLIC) method [17,5] for approximating the indicator function. In the PLIC method, the interface inside a computational cell is represented by a plane. While, in general, any PLIC reconstruction scheme can be used for the presented method, this work relies on the Mixed-Youngs-Centered scheme proposed by Aulisa et al. [1]. Lagrangian cell estimation: The LC is defined by the material evolution of the EC. This results in an arbitrary curved shape as depicted in Figure 1 (left). LP proposed to discretize the LC by six tetrahedra, while the vertex connectivity from the EC is maintained. As a result, the Discrete Lagrangian Cell (DLC), is fully characterized by the position of the transported vertices, which can be obtained by numerically integrating Eq. (5) backward in time. For this, a fourth-order Runge-Kutta scheme is used in this work. Intermediate velocities are obtained from tri-linear interpolation in space and assuming the field as frozen. Remapping: The remapping step aims to compute the liquid phase volume embedded inside the DLC. As the DLC is described as an aggregate of tetrahedra, the task is to compute the reference phase volume inside each tetrahedron. Furthermore, the PLIC approximation is discontinuous so each tetrahedron first needs to be trimmed by the planes describing the EC before the resulting convex polyhedron can be further cut by the PLIC representation of the corresponding EC. Finally, the volume of the complex polyhedron can be computed, which corresponds to the liquid phase volume inside the tetrahedra. Computationally, these operations are performed using the robust polyhedral library proposed by LP.

Lagrangian Cell Correction
Due to the numerical integration of Eq. (5) and the triangulation, the DLC does not maintain the volume of the EC, which ultimately results in conservation errors. To provide degrees of freedom for correcting the DLC, additional correction vertices are introduced on the LC faces. In contrast to the transported vertices, the position of the correction vertices is not governed by Lagrangian trajectories, but are obtained from conservative flux volumes. The general correction idea, which was introduced by Comminal et al. [4] in their 2D cell-based method, is schematically shown for a 2D configuration in Figure 1 (right). The correction requires that the volume of the DLC is equal to the EC volume subtracted by the signed volume of each discrete flux volume. To achieve this, this work proposes to triangulate the LC using 18 tetrahedra, where six tetrahedra are used for the estimated LC, and two additional tetrahedra are added on each cell face to allow a correction. The proposed discretization provides the minimum number of tetrahedra to ensure a cell-face to cell-face compatible triangulation with adjacent DLCs, which is required to avoid gaps or overlaps between adjacent DLCs that ultimately would result in conservation errors. The minimization of the number of tetrahedra is desired to keep computational costs low. The discrete flux volumes are formed from a combination of selected LC and EC vertices. Eight tetrahedra are used to triangulate the flux volumes, which is illustrated for the EC face 1375 in Figure 2. The used vertex ordering ensures face-to-face compatibility of adjacent flux volumes and provides a consistent sign convention similar to OD. It is straightforward to show that the proposed discretization of the LC and the flux volumes span exactly the same volume. As the presented flux volume discretization based on a subset of LC and EC vertices is equivalent to the flux volume discretization proposed by OD, this proves the discrete equivalence to the method of OD. For this, it is assumed that the same indicator function approximation and the same vertex transport are used. A fundamental difference between both methods is, however, that the remapping is performed for the tetrahedra of the DLC instead of the discrete flux volumes. It can be shown that at least two times fewer tetrahedra are required for remapping, which is the main driver of the computational costs and scales with the number of tetrahedra. Consequently, it follows that the proposed method reduces computational costs. To obtain the positions of the correction vertices, the ideas of OD are followed. The equality of the EC face velocity and the discrete effective flux velocity is used to define a volume deficit for each EC face that needs to be corrected by the correction tetrahedra. The only degrees of freedom are the coordinates of the correction vertex. While the coordinates tangential to the EC face are constraint to the barycenter of the projected face, the coordinate normal to the EC face can be determined from an analytical relation.

Transport Tests
The proposed method is implemented into the interfacial solver of the in-house code CIAO, which has been successfully used for many multi-physics simulations [2,6]. The new method is verified using a 2D and a 3D transport test. In both tests, a periodic velocity field is prescribed to transport the interface instead of using a velocity field from a Navier-Stokes solver. The performance of the scheme is evaluated according to where F 0 k indicates the initial volume fraction and T the period of the flow. While the first metric, E mass , measures the change of liquid volume after one period, the second metric, E shape , quantifies unphysical interface deformations after one period. Zalesak's disk: In Zalesak's disk case [18], a notched disk is rotated in a solid body manner using the velocity field In all cases, a constant time step is used that corresponds to a CFL number of π/4. Figure 3 (left) shows the PLIC representation of the initial interface and after one rotation for the coarsest mesh used. Even for the mesh resolving the notch width by 2.5 grid points, the general shape of the disk is maintained. In Figure 3, the mass conservation error and shape error are presented for the proposed method. For comparison, the results without the proposed LC correction, which corresponds to the method of LP, and the reported values of OD normalized by the notched disk area are additionally shown. The normalization is required because OD used absolute definitions for the mass and shape error. Similar to this work, OD used a one-cell thick 3D computational domain without reporting the thickness. As a result, only the scaling is comparable. The mass conservation error is independent of the grid size and always close to machine precision for the proposed method and the method of OD. This agrees well with the expectation because both methods are discretely conservative. The method of LP shows a fifth-order scaling as reported by LP [3]. No significant differences in the shape error for the uncorrected and corrected . Left: PLIC representation of the interface before (red) and after (black) one rotation. Center/Right: Convergence of the mass conservation error (center) and shape error (right) after one rotation for the proposed method, the method of LP [3], and the method of OD [14]. Dotted lines show fifth-order (left) and second-order (right) scaling. The red dashed line indicates first-order scaling. method are observable. All presented values, including those of OD, show a convergence between first-and second-order. 3D deformation: In the 3D deformation test initially proposed by Leveque [11], a sphere is advected in a swirling flow. This test is used to demonstrate the transport properties and the robustness in under-resolved situations in a 3D configuration. The interface is transported with the velocity field  Figure 4, the PLIC representation of the interface at t/T = 0.5 and t/T = 1 is shown for three meshes. At t/T = 0.5, the sphere is at a maximum stretch, and a thin liquid phase sheet is formed. For the coarsest mesh, the grid resolution is too low to maintain the thin sheet, which results in artificial breakup. While for the 128 3 mesh an interface disturbance still can be observed, the method well captures the thin sheet on a 256 3 mesh. Furthermore, it can be observed that the interface at t/T = 1 deviates from the expected interface shape of a sphere. The discrepancies, however, reduce significantly with mesh refinement.  In Figure 5, the convergence of the conservation and shape error is displayed for the proposed method, the method of LP, and the method OD. In contrast to Zalesak's disk case, the values shown are comparable because the initial reference volume is known and is used to normalize the values of OD. Similar to Zalesak's disk case, the conservation error for the proposed method is always close to machine precision, which is also visible for the values reported by OD. The method of LP shows a third-order scaling. The shape errors shown to the right of Figure 5 reveal no significant difference, and convergence about second-order is found. Comparing the proposed method and the method of OD, minor differences are observable although the schemes are discretely equivalent. This can be related to the different vertex integration and interface reconstruction schemes that were used by OD.
Here, H is the liquid sheet height,Ū is the initial liquid bulk velocity, σ is the surface tension, ρ is the density, and η is the dynamic viscosity, where the indices l and g indicate the fluid properties of the liquid and gas phase, respectively. In contrast to the transport tests with a prescribed velocity field, the interface is advected using a velocity field from a two-phase flow Navier-Stokes solver. For this, the proposed method is coupled to the sharp interfacial flow solver of LP [10], which is available in CIAO. The interface curvature is computed using the height function methodology presented by Popinet [15], and the pressure jump is integrated into the Navier-Stokes equations using the ghost-fluid methodology [7,10]. The computational domain is a box of size 6H × 9H × 4.5H and is discretized using 256 × 384 × 192 grid points in directions x,y, and z, respectively. The time step is limited by a CFL number of 0.6. On the left of Figure 6, the PLIC representation of the interface after 30 time units t * = tH/Ū is shown, and a complex interface topology can be observed. The temporal evolution of the mass conservation error (Eq. (7)) is shown on the right of Figure 6 using the proposed method and the method of LP. While the method of LP results in a mass conservation error in the order of 10 −4 , the proposed method limits the error to 10 −12 . Note that the oberservable spikes are caused by the logarithmic plot at the zero crossing of the signed error. In addition to the mass conservation error, which is a global metric, the infinity norm of the relative LC compression or expansion, defined by E div = max(|V L k − V E k |/V E k ), is provided as a local conservation metric in Figure 6 (right). The method of LP shows E div values in the order of 10 −3 . In contrast, E div is bounded by 10 −8 with the proposed method, which is consistent with the used tolerance (10 −8 ) of the pressure solver.

Conclusions
This paper presents a conservative cell-based method for 3D unsplit VoF advection. The method relies on the backward scheme of LP that projects each EC backward in time along its Lagrangian trajectories. To ensure discrete conservation, an LC correction is introduced. For this, additional cell vertices are added, where the position of those is obtained from conservative flux volumes. This requires an equivalent LC and flux volume discretization, which is presented in this paper. The proposed method is discretely equivalent to the flux-based method of OD. However, the cell-based approach requires at least two times fewer tetrahedra for intersection and volume computation. The performance of the scheme is analyzed using benchmark transport tests. The test cases verified the discrete conservation property and its improvement over the method of LP while maintaining first-to second-order accuracy of the geometric error. The method was applied in an atomization case, and the ability to conserve mass in atomization simulations was demonstrated.