Computing volume fractions and signed distances from triangulated surfaces immersed in unstructured meshes
SSignificance and novelty
Successful simulations of surface tension-driven two-phase flows using the unstructured Volume-of-Fluid (VOF) method (cf. [27] for a recent review) depend on the accurate calculation of volumefractions. Erroneous initialization of volume fractions leads to critical instability since the VOFinitialization errors significantly amplify curvature errors [22]. Existing contributions address thisproblem [5, 16, 6, 40, 12, 23, 22] by modeling fluid interfaces using functions and combining higher-order quadratures with mesh adaptivity to ensure accuracy. Our numerical method significantlyextends existing contributions in terms of the admissible shape of the fluid interface. Such surfacescan be composed of disjoint parts, and they can be closed or open, admitting very complex shapesand configurations. We have recently successfully used the proposed methodology for initializingsimulations of experiments involving the breakup dynamics of capillary bridges on hydrophobicstripes [14], which was not possible using other contemporary methods. a r X i v : . [ phy s i c s . c o m p - ph ] F e b ighlights • Novel algorithm for computing fractions from triangulated surfaces immersed in unstructuredmeshes. • Admissible surfaces may have sharp edges and be composed of multiple disjoint parts, e.g.,given by 3D scans. • Accurate and second-order convergent results. • Efficient calculation of signed distances and inside/outside information from triangulated sur-faces on unstructured meshes. omputing volume fractions and signed distances from triangulatedsurfaces immersed in unstructured meshes
Tobias Tolle , Dirk Gr¨unding , Dieter Bothe , Tomislav Mari´c
Mathematical Modeling and Analysis Institute, Mathematics department, TU Darmstadt, Germany
Abstract
We propose a numerical method that enables, for the first time to the best of our knowledge,the calculation of volume fractions from triangulated surfaces immersed in unstructured meshes.The calculation of volume fractions is crucial for achieving numerically stable simulations of surfacetension-driven two-phase flows with the unstructured Volume-of-Fluid method (cf. [27] for a recentreview) and can be used as a discrete phase-indicator model for the unstructured Level Set / FrontTracking method [25, 45].Existing publications that address this problem [5, 16, 6, 40, 12, 23, 22] handle the complexity ofthe fluid interface shape using compositions of functions; our proposed numerical method extendsthe admissible shape of the fluid interface to a practically arbitrary shape, using triangulatedsurfaces that can be open or closed, disjoint, and model objects of technical geometrical complexity.Signed distances are calculated geometrically near the fluid interface, approximated as a trianglesurface mesh, while the inside/outside information is propagated throughout the solution domainby an approximate solution of the Laplace equation. Volume fractions are computed with second-order convergence using signed distances, either via geometrical intersections or by a polynomialapproximation. Adaptive tetrahedral decomposition of polyhedral cells and its subsequent local re-finement ensures a high level of absolute accuracy. Although primarily developed for two-phase flowsimulations and used in simulations of wetting phenomena [14], the proposed algorithm can poten-tially be used in other methods that require inside/outside information with respect to triangularsurfaces.The software implementation is available on GitLab [24].
Keywords: volume of fluid, triangular surface mesh, signed distances, unstructured mesh
1. Introduction
We present a new numerical algorithm that calculates initial conditions for simulations of two-phase flow problems for fluid interfaces of complex shapes. The initial conditions are calculated inthe form of signed distances and volume fractions from fluid interfaces approximated as arbitrarily
Email addresses: [email protected] (Tobias Tolle ), [email protected] (DirkGr¨unding ), [email protected] (Dieter Bothe ), [email protected] (Tomislav Mari´c ) haped triangular surfaces immersed into unstructured meshes. The signed distances are relevantas initial conditions for the Level Set method [42, 43] for multiphase flow simulation. Volume frac-tions on unstructured meshes are required for the unstructured Volume-of-Fluid (VOF) method (cf.[27] for a recent review). In fact, we have applied the proposed algorithms to model experimentalfluid interfaces from wetting experiments [14], which was not possible using available contemporaryapproaches that model fluid interfaces using (compositions of) implicit functions or parameter-ized surfaces. The proposed algorithm approximates the surfaces using triangle meshes that areomnipresent in Computer-Aided Design (CAD) because of their versatility: they can approximatebasic surfaces such as spheres and ellipsoids, but also surfaces of mechanical parts, disjoint surfacesin mechanical assemblies, or surfaces resulting from imaging scans.The overall simulation domain Ω ⊂ R is separated into two subdomains Ω = Ω + ( t ) ∪ Ω − ( t ),representing phase 1 and phase 2, respectively, as illustrated for a liquid drop on a surface in fig. 1.At the contact line Γ := ∂ Ω ∩ Ω + ∩ Ω − , the liquid-gas interface Σ encloses a contact angle θ withthe solid surface ∂ Ω wall . Furthermore, the normal vector n Σ of the interface Σ is oriented such thatit points into the gas phase. Typically, a continuum mechanical model is used for the descriptionΩ + ( t )Σ( t ) θ ∂ Ω ∂ Ω wall Ω − ( t )Γ n Σ Figure 1: The different domains for a liquid ( − ) drop on a solid surface surrounded by a gas (+) phase. of such fluid mechanical problems. This description is often based on a sharp interface model, asdepicted in fig. 1. With this model, the liquid-gas interface can be described using an indicatorfunction χ ( x , t ) := (cid:40) , x ∈ Ω − ⊂ R , otherwise . (1)An approximate solution of this model requires a decomposition of the solution domain into volumesthat have no volume overlaps, the closed cells Ω c , denoted byΩ ≈ ˜Ω = { Ω c } c ∈ C (2)where C = { , , , . . . , N c } is a set of indices to mesh cells. As can be seen in fig. 2, the mesh is aset of non-overlapping subsets ( cells ) Ω c ⊂ ˜Ω. With non-overlapping, we mean that the volume ofan intersection between any two cells is zero. Index sets represent the unstructured mesh data [13].We consider a set of cell corner-points P h where each point in P h is an element of R . Geometrically,2ach cell Ω c is a volume bounded by polygons, so-called faces . A global set of faces F h is defined,and each face is a sequence of indices of points in P h . In this context, we define a cell set C c as aset of indices of faces in the set of mesh faces F h . Therefore, when referring to a volume defined bythe cell, we use Ω c and its magnitude is then | Ω c | , and when we refer to the cell as an unorderedindex set, we use C c and its magnitude | C c | is the number of faces that bound the cell.Solutions of continuum mechanical problems in geometrically complex solution domains signifi-cantly benefit from unstructured meshes. For example, gradients of solution variables are resolved atgeometrically complex boundaries by employing mesh boundary layers, strongly reducing the num-ber of cells required to achieve specific accuracy. Hence, this approx reduces the overall requiredcomputational resources.As the phase indicator χ ( x , t ) given by eq. (1) contains a jump discontinuity, it poses difficultiesfor numerical simulations of two-phase flows. With Volume of fluid (VOF) methods, this non-continuous description is discretized by introducing the so-called volume fraction α c = 1 | Ω c | (cid:90) Ω c χ ( x , t ) dx. (3)The unstructured VOF methods [27] rely on the volume fraction field α c to track interface withthe advecting velocity obtained from the solution of two-phase Navier-Stokes equations in a single-field formulation. All multiphase flow simulation methods that utilize the single-field formulation ofNavier-Stokes equations approximate the phase-indicator function similarly to eq. (3). The phase-indicator approximation utilizes signed distances in the Level Set [42, 41, 43] method, the volumefractions approximate the phase indicator for the Volume-of-Fluid [8, 32, 15, 35] method.Various methods exist that compute the volume fraction α c based on the exact phase indicator χ ( x , t ). The majority of methods calculate the integral in eq. (3) numerically, as schematicallyshown in fig. 2, using numerical quadrature. Σ χ ( x , t ) = 1 χ ( x , t ) = 0ΣΩ c Ω − h ( x ) xy Ω − Figure 2: Calculating volume fractions of a circular interface by numerical integration.
Different approaches are below with increasing complexity in terms of admissible shapes of thefluid interface. The admissible shapes range from analytic descriptions of basic geometric shapessuch as spheres and ellipsoids to implicit functions (or their combinations) and more general shapesapproximated with volume meshes. 3trobl et al. [40] propose an exact intersection between a sphere and a tetrahedron, a wedge,or a hexahedron. The proposed algorithm is exact and fast, though it is limited to the sphericalinterface shape.Fries and Omerovi´c [12] represent the fluid interface as a level set and propose a higher-orderquadrature for the integral on the right-hand side of eq. (3). The parametrization of the surfaceuses roots of the implicit function found by the closest-point algorithm. Results are presented forhexahedral and tetrahedral unstructured meshes that may also be strongly deformed. Fries andOmerovi´c [12, fig. 52, fig. 53] also show results with higher-order ( >
2) convergence for the volumeintegration of an arbitrary non-linear function on hexahedral and tetrahedral meshes. However,the volume and area integration error is reported for a single function. While a relative globalvolume error between 1e −
08 and 1e −
06 is reported, no information about the required CPU timesis provided. In the approach proposed by Fries and Omerovi´c [12], fluid interfaces with complexshapes are modeled as a composition of implicit functions.Kromer and Bothe [22] propose an efficient third-order accurate quadrature for the eq. (3).Contrary to Jones et al. [19], who decompose cells into tetrahedrons, Kromer and Bothe [22] locallyapproximate the hypersurface by a paraboloid based on the principal curvatures. Applying theGaussian divergence theorem to eq. (3) then yields contributions from the cell boundary and theapproximated hypersurface patch. Using the surface divergence theorem, Kromer and Bothe [22]reformulate the contribution from the hypersurface patch into a set of line integrals, where theassociated integrand emerges from the solution of a Laplace-Beltrami-type problem. The method ofKromer and Bothe [22] is directly applicable to unstructured meshes. However, locally, i.e., withina cell, the fluid interface must be C2 and simply connected.Aulisa et al. [3] and Bn`a et al. [5, 6] calculate the volume fraction by representing the indicatorfunction as a height function inside cubic cells, using the structure of the underlying Cartesianmesh. Numerical integration of the height function is illustrated by fig. 2. However, extendingthis approach to unstructured meshes raises many questions. First, constructing a height functionin a specific direction is complex and computationally expensive [33]. Second, the orientation ofthe interface in the chosen coordinate system may easily make the problem ill-conditioned. Finally,required mesh-search operations are complicated as the face normals of polyhedral cells are typicallynot aligned with the coordinate axes.ΣΩ c ∩ Ω − Ω c ∩ Ω − Σ Figure 3: Polyhedral cell (left) and non-convex cell (right) for which the intersection volume (dark grey) has to becomputed. The light grey regions lead to cases that have to be identified and require special treatment increasingthe problem complexity far beyond a simple one dimensional integration. geometrical computes minimal signed distances from ˜Σ. This cal-culation is relatively straightforward on structured meshes [38, 39], but significantly more complexon unstructured meshes [25, 45]. Here we significantly extend the calculation of signed distancesfrom [25, 45] by introducing an efficient approximate propagation of the inside/outside informationfrom ˜Σ.Volume fraction calculation methods outlined so far model the fluid interface using exact func-tions and handle more complex interface shapes via combinations of these functions. A combinationof exact functions cannot accurately capture the shape of the fluid interface in many cases. For ex-ample, when the interface shape is prescribed experimentally Hartmann et al. [14].One approach exists that can handle arbitrarily complex interface shapes. In this approach, thefluid interface encloses a volumetric mesh as its boundary surface mesh. This mesh given by thefluid interface is intersected with a ”background” mesh that stores volume fractions. This approachis called volume mesh intersection . An example for such an intersection between ˜Ω and cells from˜Ω − is shown in fig. 4. In principle, this approach is relatively straightforward, provided an accurategeometrical intersection of tetrahedrons is available. However, geometrical operations based onfloating-point numbers are not stable and can lead to severe errors [47, chap. 45]. ˜ΣΩ c Ω c ˜Ω − l ˜Ω − l Figure 4: Calculating volume fractions from a circular interface by volume mesh intersection.
Ahn and Shashkov [1] have initialized volume fractions by volume mesh intersection as shown infig. 4. In this approach, the approximated phase ˜Ω − ( t ) is decomposed into volumes (an unstructuredmesh), equivalently to the decomposition ˜Ω given by eq. (2). The boundary ∂ Ω − is the fluid interfaceΣ( t ), and it is approximated as a polygonal surface mesh, leading toΩ − ≈ ˜Ω − := { ˜Ω − l } l ∈ L , (4)i.e. an approximation of Ω − . Generally, as shown in the detail in fig. 4, a cell Ω c of the backgroundmesh ˜Ω may overlap with multiple cells Ω l from the ˜Ω − mesh, and vice versa. We define a set ofindices l of cells ˜Ω − l in ˜Ω − that overlap with the cell Ω c : the so-called cell stencil of Ω c in ˜Ω − l ,namely S (Ω c , ˜Ω − ) = { l ∈ L : Ω c ∩ ˜Ω − l (cid:54) = ∅ , where Ω c ∈ ˜Ω , ˜Ω − l ∈ ˜Ω − } , (5)5here L is an index set, containing indices of cells from ˜Ω − . Volume fractions { α c } c ∈ C can then becalculated by performing the intersection α c = | ∪ l ∈S (Ω c , ˜Ω − ) Ω c ∩ ˜Ω − l || Ω c | . (6)Since each ˜Ω − l overlaps with at least a one cell from ˜Ω, and we can approximate the number of cellsfrom ˜Ω that intersect each cell from ˜Ω − as N ( ˜Ω − , ˜Ω) ≈ | ˜Ω − | mean l ∈ L ( |S ( ˜Ω − l , ˜Ω) | ) , (7)where | ˜Ω − | denotes the number of cells in the mesh ˜Ω − . The average number of cells Ω c overlapping˜Ω − l , mean l ∈ L | C ( ˜Ω − l , ˜Ω) | , depends on the mesh densities of both meshes, ˜Ω and ˜Ω − . However, we doknow that mean l ∈ L | C ( ˜Ω − l , ˜Ω) | >
1. Next, we know that | ˜Ω − | grows quadratically in 2 D and cubicallyin 3 D with a uniform increase in mesh resolution, taken as the worst case scenario. It grows linearlyin 2 D and quadratically in 3 D if ˜Ω − is refined only near the interface ˜Σ := ∂ ˜Ω − . Consequently, thecomputational complexity of the volume mesh intersection algorithm in terms of cell/cell intersec-tions is quadratic in 2 D and cubic in 3 D in the worst case, and linear in 2 D and quadratic in 3 D if local refinement is used to increase the resolution of ˜Σ. The quadratic complexity in 3 D is a seri-ous drawback of this algorithm, especially for large simulations where | ˜Ω − | easily reaches hundredthousand cells per CPU core. Menon and Schmidt [30] have extended the volume mesh intersectionalgorithm from Ahn and Shashkov [1] to perform a volume conservative remapping of variables inthe collocated Finite Volume Method (FVM) with second-order accuracy on unstructured meshes.Their results confirm the polynomial computational complexity in terms of absolute CPU times forthis volume mesh intersection algorithm [30, table 3].L´opez et al. [23] propose a volume truncation algorithm for non-convex cells and apply it tothe initialization of volume fractions from exact functions on unstructured meshes. Cell-subdivisionis introduced to handle cases for which the interface crosses an edge of a cell twice. Non-planartruncated volumes are triangulated [23, fig 18], and second-order accuracy is demonstrated in termsof the relative global volume error for a uniform resolution and a higher-order accuracy when locallyrefined sub-grid meshes are used.Ivey and Moin [16] initialize volume fractions on unstructured meshes using tetrahedral decom-position of non-convex cells and perform geometrical intersections with a similar approach as thefrom Ahn and Shashkov [1]. Unlike Ahn and Shashkov [1], Ivey and Moin [16] compute volumefractions of intersected tetrahedrons by intersecting them with exact signed distance functions thatare used to model the fluid interface. Therefore, this algorithm cannot directly utilize arbitrarilyshaped interfaces. However, their approach utilizes a linear interpolation of intersection points be-tween the tetrahedron and the signed-distance function and yields second-order accuracy. Accuracyis further increased using adaptive mesh refinement.The approaches reviewed so far require an exact representation of the interface using explicitanalytic expressions, which hinders the direct application of such algorithms to initial conditionsresulting from experiments as these are typically not available as function compositions. The volumemesh intersection algorithm [1] is flexible but computationally expensive, and it requires highlyaccurate and robust geometrical intersections.The following sections outline the proposed algorithm that uses an unstructured surface mesh˜Σ to compute signed distances and volume fractions on unstructured meshes. Relying on unstruc-6ured surface meshes retains the ability to handle arbitrary-shaped surfaces while avoiding com-putationally expensive cell/cell intersections. Of course, using surface meshes to approximate thefluid interface renders the proposed algorithm second-order accurate; however, the accuracy in theabsolute sense achieved using local mesh refinement [7, 11]. The proposed algorithm geometricallycomputes signed distances near the fluid interface. These signed distances (so-called narrow-band signed-distances) are then propagated throughout ˜Ω by an approximate solution of a diffusion equa-tion. The propagated signed distances determine the value of the phase indicator χ ( x , t ) in thosecells that are either completely empty ( α c = 0), or completely full ( α c = 1). Finally, second-orderaccurate volume fraction values are calculated in intersected cells (0 < α c <
2. Surface mesh / cell intersection algorithm
The calculation of volume fractions by the proposed Surface Mesh Cell Intersection/Approximation(SMCI/A) algorithm, outlined in fig. 5, requires signed distances to the interface at cell centres andcell corner points. As a naive computation is computationally expansive (section 2.2), we employan octree based approach to the calculation of signed distances. Starting point of the octree basedsearch is the calculation of search radii at the relevant points.
In the first step, a search radius r c and r p is calculated at each cell center and cell-corner point,respectively. This is illustrated in fig. 5a. Here, the cell search radius r c is defined by r c = λ s min f ∈ F c (cid:107) x f,O − x f,N (cid:107) , (8)where x c is the cell center, λ s > search radius factor detailed below and x f,O , x f,N arethe cell centers of two cells that share the face with index f of the cell Ω c ( O for owner cell witha smaller cell index than the neighbor cell N ). Here, the index set F c contains the indices of thosefaces that form the boundary of Ω c . Based on (8), the corner-point search radius r p is defined by r p = λ s min c ∈ C p ( x p ) r c , (9)where x p is the cell-corner point, while the point-cell stencil is the index set S ( x p , ˜Ω), that containsindices of all cells from ˜Ω whose corner-point is x p .The search radii introduced above are used to define search balls in 3 D (circles in 2 D ), whichare used to reduce the number of calculations to determine signed distances between the cell cornerpoints x p and the cell centers x c with respect to the provided surface mesh ˜Σ.7 p x p ∈ P h x c r c (a) Calculation of search radii. ˜Σ octree n ˜Σ (b) Octree sub-division of the sur-face mesh ˜Σ bounding-box. φ + c n ˜Σ φ − c φ + p φ − p (c) Narrow-band signed distancesfrom the search radii and the octree. + − − ++ − − ++ + + ++ + + ++ − − + − ++ + ++ + ++ + + − − + (d) Positive and negative sign diffu-sion throughout ˜Ω. + − − ++ − − ++ + + ++ + + ++ − − + − ++ + ++ + ++ + + − − +˜ χ c = 1˜ χ c = 0 (e) Phase indicator ˜ χ is 1 or 0 incells strictly inside/outside of ˜Σ, re-spectively. + − − ++ − − ++ + + ++ + + ++ − − + − ++ + ++ + ++ + + − − +˜ χ c = 1˜ χ c = 0 (f) Computing the approximatedphase indicator ˜ χ in intersectedcells. Figure 5: Steps of the Surface Mesh Intersection / Approximation (SMCI/A) algorithms.
In contrast to various other approaches for volume fraction initialization, the here interface isnot represented by some kind of function, but as a set of triangles. To define the interface ˜Σ, wefirst denote the convex hull of a set of n points P n = { x , . . . , x n } , x i ∈ R byconv( P n ) := (cid:40) x ∈ R : x = (cid:88) x i ∈ P n γ i x i , n (cid:88) i =1 γ i = 1 (cid:41) . (10)Using this, a triangle is defined as the convex hull of a point triple: T := conv( P ). Consequently,the surface mesh is defined as ˜Σ := {T , T , . . . , T n } . (11)With the structure of ˜Σ in mind, we want to emphasize why an octree based approach is the keyto obtaining reasonable computation times. Consider the case where a minimal distance betweena point x and ˜Σ would be calculated for each cell center x c and cell-corner point x p . The need8or the spatial subdivision and search operations becomes obvious, as this would require a distancecomputation between each point of the interface mesh and each cell centers and cell corner points ofthe background mesh. Consequently, this would require | C || ˜Σ | operations to compute the geometricsigned distances at cell centers and additional computations for evaluating signed distances atcell-corner points. For our computations below, the number | C | often reaches the order of 1 e | ˜Σ | is typically on the order of 1 e
04 per CPU core. Aiming at redistancingcomputations for a dynamic setting in multiphase flows where ˜Σ = ˜Σ( t ), such a large number ofdistance computations makes such a brute force redistancing approach prohibitively expensive.The first step of the signed distance calculation is the computation of an Axis-Aligned BoundingBox (AABB) from the surface mesh ˜Σ. The AABB is used to build an octree data structure,illustrated as a 2 D quadtree subdivision in fig. 5b, which is used to access ˜Σ. The octree datastructure enables fast search queries involving cell centers and cell corner-points that are close tothe surface mesh ˜Σ, with a logarithmic computational complexity with respect to the number ofvertices in ˜Σ [28, 29]. The structure of the octree depends on the ordering of vertices in ˜Σ: since˜Σ is an unstructured surface mesh, its vertices are generally sufficiently unordered, which makesthe octree well-balanced. Once the octree has been constructed, it can be used to find the closestpoints x ∈ ˜Σ to cell centres x c and cell corner points x p . Note that this is only true for those x c , x p which are sufficiently close to ˜Σ in terms of their search radius r c , r p . Thus, the search radii definea so-called narrow band around ˜Σ, where the nearest distances are calculated geometrically. Wedenote the narrow band of ˜Σ with N ( ˜Σ), and the closed ball B ( x ∗ , r ) := { x ∈ R | (cid:107) x − x ∗ (cid:107) ≤ r } with a radius r around a point x . Then N ( ˜Σ) := (cid:110) x ∈ R | ∃ T ∈ ˜Σ such that T ∩ B ( x , r ) (cid:54) = ∅ (cid:111) , (12)where r is either r p or r c .For a point x ∈ N ( ˜Σ), the octree provides the closest point x min ∈ T min for some T ∈ ˜Σand the corresponding triangle T min itself. While the absolute distance can be directly computedas (cid:107) x − x min (cid:107) , care must be taken when computing the sign with respect to the orientation of˜Σ. Directly using the triangle normals n T may lead to false signs and consequently, to erroneousvolume fractions. Thus, we follow the work of [44, 4] and compute angle weighted normals n x v = (cid:80) T ∈ ngh( x v ) β T n T (cid:80) T ∈ ngh( x v ) β T (13)at the vertices x v of ˜Σ. Here, ngh( x v ) denotes the set of all triangles containing x v , n T a trianglenormal and β T the inner angle of T at x v . Baerentzen and Aanaes [4] propose a classificationof the point x min whether it is located within a triangle, on an edge, or a vertex and base thechoice of the normal on this classification. While such a classification is simple in theory, a robustimplementation is difficult due to the limited precision of floating point arithmetic. Thus, we optfor a linear interpolation of n x v within T min to x min , denoted n I ( x min , T min ). With this normalcomputation, the signed distance between x and x min is calculated by φ g ( x , ˜Σ) = sign(( x − x min ) · n I ( x min , T min )) (cid:107) x − x min (cid:107) . (14)where the supindex g indicates a geometrical construction. This procedure is illustrated in fig. 5c.The robustness of this approach with regard to inside/outside classification is demonstrated insection 4.3. 9sing the spatial subdivision provided by the octree, the computational complexity for findingthe minimal distances between mesh points and ˜Σ is reduced severely, as the vast majority ofcell centers x c are not even considered for calculation as no triangle T ∈ ˜Σ exists within thecorresponding search ball. The closest triangles of those points x c , whose ball B ( x c , r c ) intersects˜Σ are found with logarithmic search complexity with respect to | ˜Σ | . This significant reductionof complexity can potentially enable a future application of the proposed algorithm on movinginterfaces ˜Σ( t ) as a geometrically exact marker field model for unstructured Front Tracking methods.Therefore, it is crucial to understand that the min T ∈ ˜Σ operation in eq. (14) throughout this textrelies on the octree spatial subdivision and search queries. After the calculation of geometric signed distances in the narrow band around ˜Σ, the signed dis-tances are propagated to the bulk of different phases, as shown in fig. 5d. In [25, 45], the geometricsigned distances are set to large positive numbers throughout the domain, and a graph-traversalalgorithm is used to iteratively correct the signs of signed distances using face-cell and point-pointgraph connectivity provided by the unstructured mesh. Graph-traversal is computationally expen-sive and complicated to implement in parallel. Here we propose a straightforward alternative thatinstantaneously propagates signs of signed distances through the solution domain and is parallelizedeasily. We rely on the diffusion equation for the signed distances, namely − ∆ φ = 0 , ∇ φ = 0 , for x ∈ ∂ Ω (15)and its discretization using the unstructured finite volume method in OpenFOAM [17, 21, 31],giving a linear system of equations. The key idea to sign propagation is to apply a few iterations( <
5) of an iterative linear solver to this system. In our case a Conjugate Gradient approach withan incomplete lower upper preconditioner has been used. With the initial field set to φ ( x ) = (cid:40) φ g ( x , ˜Σ) , if x ∈ N ( ˜Σ)0 , otherwise, (16)this small number of iterations suffices to properly propagate sign( φ ) with respect to the orientationof ˜Σ throughout ˜Ω. Prerequisite for this approach to work is that the narrow band has a certainminimum width in interface normal direction. At least four cells on each side of the interface arerequired to ensure a robust propagation. This is achieved by setting a global search radius factor λ s := 4 in eq. (8) used to calculate r c at cell centers. Note that increasing λ s beyond this valueonly increases computational costs, and does not impact the accuracy of the proposed algorithm, aswith a larger value of λ s the narrow band N (Σ) becomes wider and consequently the geometricalsigned distances are calculated at more points x c , x p , using eqs. (17) and (20), respectively.Two aspects have to be considered when solving the linear system of equations resulting fromthe discretization of eq. (15). First, cells for which x c ∈ N ( ˜Σ) have to be excluded from the vector ofunknowns as φ g ( x c ) is already known for those. Second, for cells away from N ( ˜Σ) the only relevantinformation is sign( φ c ) indicating Ω c ∈ Ω − or Ω c ∈ Ω + , respectively. A few iterations of a linearsolver suffice to reliably propagate sign( φ c ) to the entire domain. The resulting field is φ c = (cid:40) φ gc , if x c ∈ N ( ˜Σ) ,φ ac , otherwise, (17)10ith φ gc denoting geometric signed distances and φ ac approximate values from the solution of eq. (15)carrying inside/outside information but without geometric meaning.Once the cell-centered signed distances φ c are computed, they are used to calculate the signeddistances at cell corner-points via φ Ip = (cid:88) c ∈ C p w p,c φ c , (18)where C p is the index set of cells that contain the cell corner point x p and the supindex I indicatinginterpolation. Furthermore, w p,c is the inverse-distance weighted (IDW) interpolation weight w p,c = (cid:107) x c − x p (cid:107) − (cid:80) ˜ c ∈ C p (cid:107) x ˜ c − x p (cid:107) − . (19)As with φ c , the accuracy of φ p is irrelevant outside of the narrow band of ˜Σ, only the sign ofthe signed distance is important in the bulk. To correct for the error introduced by the IDW-interpolation in eq. (18), signed distances at cell-corner points of intersected cells are calculatedgeometrically φ p = (cid:40) φ gp , if x p ∈ N ( ˜Σ) ,φ Ip , otherwise. (20)Equations (17) and (20) define the final signed distances at cell centers and cell-corner points,respectively. These quantities will have the value of a geometrical distance to ˜Σ in the narrow band,while outside of the narrow band only the correct sign resulting from the approximative solution ofeq. (15) is of relevant. Once the signed distances at cell centers { φ c } c =1 , ,..., | ˜Ω | and cell corner points { φ p } p =1 , ,... | P h | arecalculated as outlined in the previous section, the SMCI algorithm calculates the volume fractionsin a straightforward way. The volume fraction calculation is shown schematically for the SMCIalgorithm in fig. 6b. Each cell is decomposed into tetrahedrons, using the cell centroid x c as thebase point of the tetrahedron, the centroid of the face x c,f , and two successive points from the cell-face, x c,f,i , x c,f,i +1 . The resulting tetrahedron has the distance φ c associated to the cell centroid, thedistance φ c,f associated to the face centroid, and and ( φ c,f,i , φ c,f,i +1 ) pair of distances associatedwith a pair of points that belong to the cell-face ( c, f ), as shown in fig. 6b. If all the distancesof the tetrahedron are negative, the tetrahedron lies in the negative halfspace with respect to ˜Σ,and its total volume contributes to the sum of the volume of phase 1 inside the volume Ω c . If apair of distances in a tetrahedron has different signs, the tetrahedron is intersected by the interfaceapproximated by the surface mesh ˜Σ. The volume of this intersection is calculated by geometricallyintersecting the tetrahedron with those triangles from ˜Σ, that have a non-zero intersection with aball B enclosing the tetrahedron. The center of the ball B c,f,i := B ( x c,f,i , Rc, f, i ) is the centroid ofthe tetrahedron x c,f,i = 0 . x c + x c,f + x c,f,i + x c,f, mod( i +1 , | F c,f | ) ), where i = 0 , . . . , | F c,f | −
1, and F f is the oriented set of indices of the points x (cf. fig. 6b) that belong to the face f of the cell Ω c .The radius of the tetrahedron-ball B c,f,i is then R c,f,i = max( (cid:107) x c − x c,f,i (cid:107) , (cid:107) x c,f − x c,f,i (cid:107) , (cid:107) x c,f,j − x c,f,i (cid:107) , (cid:107) x c,f, mod( j +1 , | F c,f | ) − x c,f,i (cid:107) ) , (21) j = 0 , . . . , | F c,f | −
1. This sub-set of ˜Σ is found using the octree data structure with logarithmiccomplexity with respect to ˜Σ, as outlined in the previous section. For the example tetrahedron11n the cell shown in fig. 6b, the resulting intersection between the approximated interface ˜Σ anda tetrahedron from the cell Ω c is shown as the shaded volume. The magnitude of this volume iscomputed by applying the Gauss divergence theorem using eq. (31). The phase-specific volumesfrom cell-tetrahedrons are summed up for the cell Ω c , into the total phase-specific volume of thephase 1 within the cell Ω c , and the volume fraction is therefore computed as α c = (cid:80) f =0 ,... | C c |− (cid:80) i =0 ,..., | F c,f |− | T ( x c , x c,f , x c,f,i , x c,f, mod( i +1 , | F c,f | ) ) ∩ ( B c,f,i ∩ ˜Σ) || Ω c | (22)with T := { x , x , x , x } denoting a tetrahedron. The SMCI algorithm is summarized by algo- x c (a) A cell Ω c intersected by ˜Σ. x c , φ ( x c ) φ ( x c,f,i +1 ) φ ( x c,f,i ) x c,f,i R c,f,i x c,f , φ ( x c,f ) x c,f,i +1 x c,f,i (b) Tetrahedral cell decomposition. Figure 6: Centroid decomposition of an interface cell into tetrahedra and calculation of α c using the SMCI/Aalgorithms. rithm 1.
3. Surface-Mesh / Cell Approximation algorithm
This section presents an alternative approach to the computation of volume fractions presentedin section 2.4. While section 2.4 details a method based on geometric intersections, this sectionintroduces an algorithm based on volumetric reconstruction by adaptive mesh refinement. Detrixheand Aslam [10] introduce a second order accurate approximation for the volume fraction of a triangle(2D) or a tetrahedron (3D). Their model is an algebraic expression taking the signed distances φ ofthe vertices as arguments. In contrast, we propose a volume fraction initialization algorithm thatemploys this model in combination with an adaptive tetrahedral cell decomposition and the octree-based signed distance calculation described in section 2. We term this algorithm Surface-Mesh/CellApproximation (SMCA) and it is outlined below.The SMCA-algorithm is based on the signed distance results of the SMCI-algorithm introducedin section 2. The steps depicted in fig. 5a - 5d of the SMCI/A are used to compute φ c , φ p in thenarrow band and propagate inside/outside information in the rest of the mesh points. Subsequentsteps for the computation of volume fractions are displayed in fig. 7. First, all cells intersected by ˜Σ12 lgorithm 1 The Surface-Mesh / Cell Intersection Algorithm (SMCI) α c = 0, φ c,p = 0 Compute search radius for cell centers r cc ∈ C using eq. (8). for cell centroids { x c } c ∈ C do Place the vertices of ˜Σ into an octree (section 2.2). Find the triangle T n ∈ ˜Σ nearest to x c within a ball B ( x c , r c ). Set φ gc := φ g ( x c , T n ) using eq. (14). end for Approximately solve eq. (15) to propagate sign ( φ c ). Compute search radius for cell corner points r pp ∈ P using eq. (9). Find all intersected cells I = { c, φ c φ p < } . Use eq. (17) to correct φ c within the narrow band. Compute φ p in the bulk using eq. (18). Use eq. (20) to correct φ p within the narrow band. for cells { Ω c } c ∈ C do if φ c ≤ φ p ≤ then (cid:46) Cell is inside the negative ˜Σ-halfspace. α c = 1 end if if cell Ω c is intersected, c ∈ I then (cid:46) Cell is intersected by ˜Σ. α c given by eq. (22). end if end for are identified to reduce computational costs, as only these cells have intermediate values 0 < α c < x c ∈ N ( ˜Σ) is checked with the bounding ballcriterion . We define a bounding ball (bb) for a point x bb ∈ Ω c using r bb = max x ∈ Ω c (cid:107) x − x bb (cid:107) .This ball is the smallest ball that contains all points of Ω c . We compare this bounding ball to B ( x bb , | φ ( x bb )). These balls are shown in fig. 8, where the bounding ball is illustrated by a dashedand the other ball by a continuous line. As a general observation, if the bounding ball is containedin the ball with the radius | φ ( x bb ) | , i.e. B ( x bb , r bb ) ⊆ B ( x bb , | φ ( x bb ) | ), then such a cell is guaranteedto be a bulk cell. This cell can then be removed from the set of cells in the narrow band to reducethe number of cells which are considered for decomposition in the next step. If the criterion is notsatisfied, the cell is considered an interface cell. Two remarks on this criterion: first, the existenceof such a x bb is not a necessary but a sufficient condition. Second, in a practical implementationevaluation of this criterion is only feasible for a small number of points when aiming to keepcomputational costs reasonable. Thus, the actual check is performed by evaluating f bb ( x , φ x , Ω c ) = (cid:40) , max x i ∈ Ω c (cid:107) x i − x (cid:107) ≤ | φ x | , , otherwise (23)with x ∈ Ω c . The evaluation of the max-operator is based on a comparison to the corner points x i of the cell Ω c . For example, in our implementation this function is only evaluated at cell centres x c (original mesh cells, see below) or cell corner points (tetrahedra resulting from decomposition). Asa consequence, a few bulk cells are considered as interface cells (fig. 8b). We deem this acceptable asthis only has a minor impact on the computational time, but not on the computed volume fractions.After identification of interface cells, the cell volume fractions are initialized according to the13 a) Identify potential interface cells(marked grey) using bounding ballcriterion. Shown are circles withradii | φ c | . (b) Adaptive, tetrahedral decompo-sition of interface cells. Compute φ at new vertices. (c) Compute the volume fraction α c using the model of Detrixhe andAslam [10] (detail view). Figure 7: Steps of the SMCA algorithm following signed distance computation and inside/outside propagation. sign of φ c , α c = (cid:40) , φ c ≤ , , otherwise . (24)This gives correct volume fractions for bulk cells, while the values of interface cells are updated asdescribed below. Each cell flagged as an interface cell by the method described above is decomposedinto tetrahedra using its centroid and cell face centroids as shown in fig. 6. Each resulting tetrahedronis further refined in an adaptive manner such that resolution is only subsequently increased wherea new tetrahedron is again intersected by the interface. To achieve this, a tetrahedron T is checkedwith the bounding ball criterion eq. (23). The criterion is only evaluated at the vertex x max ∈ T for which | φ ( x max ) | = max x ∈ T | φ ( x ) | . Only if f bb ( x max , φ, T ) = 0 (eq. (23)), T is considered forfurther decomposition. An obvious choice would be decomposition at the centroid of T . However,repeated application of this approach results in increasingly flattened tetrahedra. To avoid thisproblem, we apply the decomposition shown in fig. 9. First, from the vertices edge centres of thetetrahedron x ij = 12 ( x i + x j ) , i, j ∈ { , , , } , i (cid:54) = j (25)are computed (fig. 9a). By combining each vertex x i with the three edge centres of the adjacentedges, four new tetrahedra are created (fig. 9b). The remainder of the original tetrahedron is anoctahedron (fig. 9b grey dashed lines) constituted by the edge centres x ij . This octahedron is de-composed into four additional tetrahedra by choosing two opposite edge centres as shown by theblack line in fig. 9c. The indices of vertices of such a line are the numbers one to four. From theremaining four edge centres, point pairs are created such that { x mn , x mo } or { x mn , x on } , yieldingfour pairs. Combining each pair with { x ij , x kl } (e.g. black edge in fig. 9c) gives the aforementionedfour tetrahedra. Subsequently, φ is computed for the added vertices x ij . The decomposition is basedon the pair of edge centres that have the smallest distance between each other. Refinement is com-pleted when a maximum refinement level l max is reached. This can either be an arbitrary prescribedvalue or can be computed such that the edge length of the refined tetrahedra is comparable to the14 r bb | φ c | x c (a) Bulk cell: the ball B ( x c , | φ c | ) contains the cellbounding ball B ( x c , r bb ). Σ | φ c | x c r bb (b) False positive: a bulk cell which is not detected bythe bounding ball criterion as B ( x c , r bb ) (cid:42) B ( x c , | φ c | ) . Figure 8: Illustration of the idea of the bounding ball criterion in 2D for clarity. The solid grey line represents B ( x c , | φ c | ), the grey dashed one B ( x c , r bb ). edge length of surface triangles. In the latter case, l max = min l ∈ N (cid:18) L tet L tri < l (cid:19) (26)with L tet and L tri being cell specific reference lengths for tetrahedra and surface triangles, respec-tively. Different choices for L tet and L tri are possible. We choose L tet = 1 n t (cid:88) e ∈ E cdc | e | ,L tri = min e ∈ E ˜Σ ,c | e | with E cdc denoting the set of edges resulting from tetrahedral decomposition of a cell Ω c at itscentroid, n t the number of edges in E cdc and E ˜Σ ,c a subset of edges of ˜Σ. The set E ˜Σ ,c consists ofall edges of T ∈ ˜Σ for which
T ∩ B ( x cp , r cp ) (cid:54) = ∅ . Here, x cp = 1 | P cp | (cid:88) x i ∈ P cp ,P cp := { x ∈ ˜Σ : min x i ∈ Ω c (cid:107) x − x i (cid:107) } and the radius r cp = max x ∈ P cp (cid:107) x − x cp (cid:107) .Finally, after computing a tetrahedral decomposition of each interface cell, the volume fractionof a cell Ω c is calculated as α c = 1 | Ω c | (cid:88) T ∈ T c α ( T ) | conv( T ) | (27)15 x x x x x x x x x (a) Original tetrahedron withvertices ( x i , black) and edgemidpoints ( x ij , grey). (b) Four tetrahedra are created bycombining each vertex with itsconnected edge midpoints (indicatedby dashed lines). x x x x x x (c) Decompose octahedron into fourtetrahedra by combining each greyedge with the black line formed bytwo opposite points (here x , x ). Figure 9: Decomposition of a tetrahedron into eight tetrahedra using edge midpoints. where T c denotes the set of tetrahedra resulting from the decomposition of Ω c and | conv( T ) | thevolume of T . The volume fraction α ( T ) is computed with the approach of Detrixhe and Aslam [10](eq. 7), repeated here α ( T ) = , φ ≤ , − φ ( φ − φ )( φ − φ )( φ − φ ) , φ ≤ < φ , − φ φ ( φ + φ φ + φ ) + φ φ ( φ φ − ( φ + φ )( φ + φ ))( φ − φ )( φ − φ )( φ − φ )( φ − φ ) , φ ≤ < φ , − φ ( φ − φ )( φ − φ )( φ − φ ) , φ ≤ < φ , φ > , (28)where φ ≥ φ ≥ φ ≥ φ are the signed distances at the vertices x i of T . The overall approach issummarized in algorithm 2. Algorithm 2
The Surface-Mesh / Cell Approximation Algorithm (SMCA) Follow algorithm 1 up to step 13. Identify interface cells (eq. (23)) Set bulk α c (eq. (24)) Centroid decomposition of cells into tetrahedra (fig. 6) for l ∈ { , . . . , l max } do Flag tetrahedra for further refinement (eq. (23)) Decompose flagged tetrahedra (fig. 9) Compute φ for new points (eq. (14)) end for Compute α c for interface cells (eq. (27)) . Results We use the difference between the total volume given by the volume fraction calculated fromthe surface on the unstructured mesh, and the exact volume bounded by the surface, namely E v = 1 V e (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) V e − (cid:88) c ∈ C α c | Ω c | (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (29)as the measure of accuracy of the proposed algorithms. Here, V e is the volume given by the exactsurface function, or the volume that is bounded by a given surface mesh if an exact surface functionis not available, e.g. in sections 4.2 and 4.3. In these cases, we calculate V e using V e = 13 (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) V e ∇ · x dV (cid:12)(cid:12)(cid:12)(cid:12) = 13 (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) ∂V e x · n dS (cid:12)(cid:12)(cid:12)(cid:12) (30)where ∂V e is the surface that bounds V e . As this surface is triangluated, eq. (30) can be expandedfurther V e = 13 (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:88) t ∈ ..N ˜Σ (cid:90) T t x · n dS (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = 13 (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:88) t ∈ ..N ˜Σ (cid:90) T t ( x − x t + x t ) · n dS (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = 13 (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:88) t ∈ ..N ˜Σ x t · S t (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (31)where N ˜Σ is the number of triangles in ˜Σ, T t ∈ ˜Σ are triangles that form the interface mesh, and x t , S t are their respective centroids and area normal vectors. Computing architectureCPU vendor id : AuthenticAMDcpu family : 23model : 49model name : AMD Ryzen Threadripper 3990X 64-Core Processorfrequency : 2.90 GHzCompiler version : g++ (Ubuntu 10.2.0-5ubuntu1 20.04) 10.2.0optimization flags : -std=c++2a -O3
Table 1: Used computing architecture.
Table 1 contains the details on the computing architectures used to report the absolute CPUtimes in the result section. We have fixed the CPU frequency to 2.9GHz to stabilize the CPU timemeasurements.
Exact initialization algorithms for spheres are available on unstructured meshes [40, 22]. We usethe sphere and ellipsoid test cases to confirm the second-order convergence of SMCI/A algorithmsand their applicability as a volume fraction model for the unstructured Level Set / Front Trackingmethod [25, 45]. The sphere case consists of a sphere with a radius R = 0 .
15, and the ellipsoidhalf-axes are (0 . , . , . . , . , . .1.1. SMCI Algorithm Figure 10 shows the expected second-order convergence of the global error E v given by eq. (29)on cubic fig. 10a and irregular hexahedral fig. 10b unstructured meshes. In fig. 10, N c is the numberof cells used along each spatial dimension of ˜Ω and N T is the number of triangles used to resolvethe sphere. √ N T − − E v N c = 16 N c = 32 N c = 64 N c = 128 (a) Equidistant mesh. √ N T − − E v N c = 16 N c = 32 N c = 64 N c = 128 (b) Irregular hexahedral mesh. Figure 10: E v errors of the SMCI algorithm for the sphere. The grey dashed line indicates second order convergence. The CPU times reported in fig. 11 for the architecture A1 in table 1 show that the SMCIalgorithm is a promising candidate for a volume fraction model for the unstructured Level Set/ Front Tracking method. The complexity of the algorithm expressed in terms of the measuredCPU time remains, linear for a constant ratio √ N T /N c . The computational complexity increases toquadratic with an increasing number of triangles per cell √ N T /N c : this happens when a very finesurface mesh is used to compute volume fractions on a very coarse volume mesh. An intersectionbetween a highly resolved surface mesh and single cell of a relatively coarse mesh is shown in fig. 12a.This configuration is relevant for accurate initialization of volume fractions on coarse meshes,but irrelevant for calculating the phase indicator for Front Tracking, where only a small number oftriangles per multimaterial cell ( ≤
10) is present. Therefore, linear complexity of the SMCI algorithmfor small ratios √ N T /N c makes SMCI a potential candidate for a highly accurate geometrical volumefraction model for the unstructured Level Set / Front Tracking method. We will investigate thispossibility in our future work. When considering the absolute CPU times, it is important to notethat the SMCI algorithm has not yet been optimized for performance.The volume error E v for a sphere is shown in fig. 10b for a perturbed hexahedral mesh. Anexample perturbed mesh from this parameter study is shown in fig. 12b. The mesh is distorted byrandomly perturbing cell corner points, using a length scale factor α e ∈ [0 ,
1] for the edges e thatsurround the mesh point. We have used α e = 0 .
25, resulting in perturbations that are of the size of0 . × the edge length. This results in a severe perturbation of the mesh shown in fig. 12b, as well asnon-planarity of the faces of hexahedral cells. Still, as shown in fig. 10b, SMCI retains second-orderconvergence, which is also the case for the initialization of the ellipsoid on the equidistant fig. 13and perturbed hexahedral mesh fig. 13b. 18 √ N T C P U t i m e i n s ec o nd s N c = 16 N c = 32 N c = 64 N c = 128 Figure 11: CPU times of the SMCI algorithm for the sphere initialized on a cubic unstructured mesh. (a) SMCI: intersected cell. α c (b) SMCI: sphere and ellipsoid volume fractions. Figure 12: SMCI algorithm used with a sphere and an ellipsoid on an unstructured hexahedral mesh.
First, the effectiveness of the local adaptivity employed in the SMCA algorithm is examined witha spherical interface as described in section 4.1. Resolution of the volume mesh is fixed to N c = 16cells in each direction while the sphere is resolved with √ N T ≈
410 triangles. Maximum refinementlevels l max from 0 to 3 are manually prescribed. In fig. 14, the resulting global volume errors E v are displayed. This test case confirms the expected second-order convergence of E v with adaptiverefinement. An exemplary tetrahedral decomposition of a perturbed hex cell with a part of thethe surface mesh is displayed in fig. 15. It demonstrates that the adaptive refinement based on the19 √ N T − − − E v N c = 16 N c = 32 N c = 64 N c = 128 (a) Equidistant mesh. √ N T − − − E v N c = 16 N c = 32 N c = 64 N c = 128 (b) Irregular hexahedral mesh. Figure 13: E v errors of the SMCI algorithm for the ellipsoid. The grey dashed line indicates second orderconvergence. bounding ball criterion eq. (23) works as intended. Refinement is localized to the vicinity around theinterface. Yet, the approach ensures all tetrahedra intersected by the interface are actually refined.The effectiveness of the local adaptive refinement compared to a uniform one becomes apparent whencomparing the resulting number of tetrahedra. Our adaptive approach yields around 2247 tetrahedraper interface cell on average for the spherical interface with √ N T ≈ N c = 16 and l max = 3.A uniform decomposition, on the contrary, would result in M i × M l max r = 24 × ≈ . × tetrahedra, where M i denotes the number of tetrahedra from initial cell decomposition and M r thenumber of tetrahedra from refining a tetrahedron. Thus, the local adaptive refinement reduces therequired overall number of tetrahedra by a factor of 5 . E v (eq. (29)) are shown in fig. 16 for cubic cells (fig. 16a) and perturbed hexahedralcells (fig. 16b). Domain size, sphere centre and radius are identical to the SMCI setup as well as theperturbation factor α e = 0 .
25. The maximum refinement level is computed according to eq. (26).Both mesh types yield nearly identical results and show second-order convergence. Resolution of thevolume mesh N c has a minor influence for coarser surface meshes which vanishes for √ N T > E v are shown in fig. 17. The results are qualitatively andquantitatively similar to those of the spherical interface. Absolute computational times requiredfor the initialization of a sphere with the SMCA algorithm are displayed in fig. 18. Run times havebeen measured on the architecture listed in table 1. As the implementation SMCI algorithm, ourimplementation of the SMCA algorithm has not yet been optimized for performance. Some methods that are surveyed in section 1 can initialize volume fractions from exact im-plicit surfaces, such as a sphere or an ellipsoid, analyzed in section 4.1. One novelty of SMCI/A20 l max − − E v Figure 14: E v errors of the SMCA algorithm using different refinement levels l max for a sphere. Resolution of volumeand surface mesh are fixed to N c = 16 and √ N T ≈ α c . Tetrahedra from differentrefinement levels are shown in different colors (level 1: blue, level 2: grey, level 3: red). Due to adaptivity, thehighest refinement level is localized in the vicinity of the surface mesh.. √ N T − − E v N c = 16 N c = 32 N c = 64 N c = 128 (a) Equidistant mesh. √ N T − − E v N c = 16 N c = 32 N c = 64 N c = 128 (b) Irregular hexahedral mesh. Figure 16: E v errors of the SMCA algorithm for the sphere. The grey dashed line indicates second orderconvergence. √ N T − − E v N c = 16 N c = 32 N c = 64 N c = 128 (a) Equidistant mesh. √ N T − − E v N c = 16 N c = 32 N c = 64 N c = 128 (b) Irregular hexahedral mesh. Figure 17: E v errors of the SMCA algorithm for the ellipsoid. The grey dashed line indicates second orderconvergence. algorithms is their ability to compute volume fractions from arbitrary surfaces on arbitrary unstruc-tured meshes. For example, volume fractions given by an experimental surface were calculated bythe SMCI algorithm in Hartmann et al. [14] for studying breakup dynamics of a capillary bridge ona hydrophobic stripe between two hydrophilic stripes. In [14], the experimental setup involves a liq-uid bridge that is formed between two larger droplets across a hydrophobic stripe. The hydrophobicstripe drives the collapse of this liquid bridge, that is observed experimentally and in a simulationin [14]. The quantitative comparison of the simulation and the experiment from [14] is shown infig. 19a. The experimental surface from Hartmann et al. [14], used to initialize volume fractions, isshown in fig. 19b. The SMCI algorithm computes the volume vractions of the experimental fluidinterface from [14] with the volume error E v = 7 . e −
06. As shown in section 4.1, the accuracy22 √ N T C P U t i m e i n s ec o nd s N c = 16 N c = 32 N c = 64 N c = 128 Figure 18: CPU times of the SMCA algorithm for the sphere initialized on a cubic unstructured mesh. of the initialization depends on the quality of the surface mesh, not on the resolution of the volumemesh, that is chosen in this case to appropriately resolve the hydrodynamics in [14]. e x p e r i m e n t s i m u l a t i o n τ = Rd dd
250 µm (a) Qualitative comparison with experiment, image from[14]. f . . . . . . (b) Initialization of volume fractions f for the wettingexperiment, image adapted from [14]. Figure 19: Simulation of the wetting experiment with the fluid interface given as a triangular surface mesh [14].
To demonstrate that the SMCI/A algorithms are able to handle interfaces more complex thanshown in section 4.1 and section 4.2, the surface mesh from a CAD model displayed in fig. 20a isused. In contrast to the previous interfaces, this one features sharp edges and geometric featuresof distinctly different sizes. The mesh for this test case has been generated with the cartesianMesh a) Surface mesh from a CAD model. (b) Cross section of the volume mesh with part of the surfacemesh, colored by signed distance. Figure 20: Surface and volume mesh of the CAD model test case. tool of cfMesh [20]. Refinement is used in the vicinity of the interface. This meshing procedureis chosen to obtain a mesh that closer resembles that of an industrial application than a uniformcubic mesh. A cross section of the mesh is depicted in fig. 20b. Before examining the computedvolume fractions for this case, the signed distance calculation (section 2.2) and sign propagation(section 2.3) are verified. The presence of sharp edges (see fig. 20a) makes this test case more proneto false inside/outside classifications than the others shown so far. Yet our procedure yields thecorrect sign for the distance in all cells as shown in fig. 21a. The enclosed volume of the surfacemesh is considered as Ω + , thus φ > x ∈ Ω + . As displayed in fig. 21a and confirmed byfurther manual inspection of the results, the proposed signed distance calculation correctly classifiesall cells within the narrow band and robustly propagates this information to the entire domain. Thisis reflected in the volume fractions as computed, shown in fig. 21b. Bulk cells are assigned values ofeither 1 or 0, depending on whether they are located in Ω + or Ω − and mixed cells with 0 < α c < E v depicted infig. 22 have been obtained with the SMCA algorithm using different refinement levels. As for thespherical interface (see fig. 14), second-order convergence is achieved, even though the surface meshapproximates a non-smooth interface here.
5. Conclusions
The proposed Surface-Mesh Cell Intersection / Approximation algorithms accurately computesigned distances from arbitrary surfaces intersecting arbitrary unstructured meshes. Geometricalcalculations ensure the accuracy of signed distances near the discrete surface. The signed distances24 a) Cells for which φ c ≥ Figure 21: Inside/outside computation and resulting volume fractions for the CAD geometry. l max − − E v Figure 22: E v errors of the SMCA algorithm using different refinement levels l max for the CAD model with thereference volume V e computed by eq. (31). The grey dashed line indicates second order convergence. (actually their inside / outside information) are propagated into the bulk using the approximatesolution of a Laplace equation. Once the signed distances are available in the full simulation domain,25he SMCI algorithm computes volume fractions by intersecting arbitrarily-shaped mesh cells withthe given surface mesh, while the SMCA algorithm approximates volume fractions using signeddistances stored at cell corner points. Both algorithms are robust and show second-order convergencefor exact surfaces and arbitrarily shaped surface meshes. The SMCI algorithm scales linearly witha small number of surface triangles per cut-cell. Since a small number of triangles per cell is arequirement for Front Tracking, this linear-complexity makes SMCI an interesting candidate forcomputing volume fractions in the 3D unstructured Level Set / Front Tracking method Mari´c et al.[25], Tolle et al. [45], which will be the subject of future investigations.
6. Acknowledgments
Calculations for this research were conducted on the Lichtenberg high performance computer ofthe TU Darmstadt.Funded by the German Research Foundation (DFG) – Project-ID 265191195 – SFB 1194,Projects B02, B01 and Z-INF.We are grateful for the discussions on the phase-indicator calculation in the LCRM methodon structured meshes [39] with Prof. Dr. Seungwon Shin, Dr. Damir Juric, and Dr. Jalel Cherguiwithin the project ”Initiation of International Cooperations” MA 8465/1-1.
References [1] H. T. Ahn and M. Shashkov. Multi-material interface reconstruction on generalized polyhedralmeshes. Technical Report LA-UR-07-0656, 2007. URL http://cnls.lanl.gov/~shashkov/ .[2] E. Aulisa, S. Manservisi, and R. Scardovelli. A mixed markers and volume-of-fluid methodfor the reconstruction and advection of interfaces in two-phase and free-boundary flows.
J.Comput. Phys. , 188(2):611–639, 2003. ISSN 00219991. doi: 10.1016/S0021-9991(03)00196-7.URL https://doi.org/10.1016/S0021-9991(03)00196-7 .[3] E. Aulisa, S. Manservisi, R. Scardovelli, and S. Zaleski. Interface reconstruction with least-squares fit and split advection in three-dimensional Cartesian geometry.
J. Comput. Phys. ,225(2):2301–2319, 2007. ISSN 00219991. doi: 10.1016/j.jcp.2007.03.015. URL https://dx.doi.org/10.1016/j.jcp.2007.03.015 .[4] J. A. Baerentzen and H. Aanaes. Signed distance computation using the angle weightedpseudonormal.
IEEE Transactions on Visualization and Computer Graphics , 11(3):243–253,2005. doi: 10.1109/TVCG.2005.49. URL https://doi.org/10.1109/TVCG.2005.49 .[5] S. Bn`a, S. Manservisi, R. Scardovelli, P. Yecko, and S. Zaleski. Numerical integrationof implicit functions for the initialization of the VOF function.
Comput. Fluids , 113:42–52, 2015. doi: 10.1016/j.compfluid.2014.04.010. URL http://dx.doi.org/10.1016/j.compfluid.2014.04.010 .[6] S. Bn`a, S. Manservisi, R. Scardovelli, P. Yecko, and S. Zaleski. Vofi - A library to initializethe volume fraction scalar field.
Comput. Phys. Commun. , 200:291–299, 2016. ISSN 00104655.doi: 10.1016/j.cpc.2015.10.026. URL http://dx.doi.org/10.1016/j.cpc.2015.10.026 .267] S. J. Cummins, M. M. Francois, and Douglas B. Kothe. Estimating curvature from volumefractions.
Comput. Struct. , 83(6-7):425–434, 2005. ISSN 00457949. doi: 10.1016/j.compstruc.2004.08.017. URL http://dx.doi.org/10.1016/j.compstruc.2004.08.017 .[8] R. B. DeBar. Fundamentals of the KRAKEN code.
Tech. Rep. , pages UCID–17366, 1974.[9] S. S. Desphande, L. Anumolu, and M. F. Trujillo. Evaluating the performance of the two-phaseflow solver interFoam.
Comput. Sci. Discov. , 5(014016):1—-36, 2012. doi: 10.1088/1749-4699/5/1/014016. URL https://doi.org/10.1088/1749-4699/5/1/014016 .[10] M. Detrixhe and T. D. Aslam. From level set to volume of fluid and back again at second-order accuracy.
Int. J. Numer. Methods Fluids , 80:231–255, 2016. doi: 10.1002/fld. URL https://dx.doi.org/10.1002/fld .[11] M. M. Francois, S. J. Cummins, E. D. Dendy, D. B. Kothe, J. M. Sicilian, and M. W. Williams.A balanced-force algorithm for continuous and sharp interfacial surface tension models withina volume tracking framework.
J. Comput. Phys. , 213(1):141–173, 2006. ISSN 00219991. doi:10.1016/j.jcp.2005.08.004. URL http://dx.doi.org/10.1016/j.jcp.2005.08.004 .[12] T. P. Fries and S. Omerovi´c. Higher-order accurate integration of implicit geometries.
Int. J.Numer. Methods Eng. , 106(5):323–371, 2016. ISSN 10970207. doi: 10.1002/nme.5121. URL https://dx.doi.org/10.1002/nme.5121 .[13] S. Ghali.
Introduction to geometric computing . Springer Science & Business Media, 2008. URL .[14] Maximilian Hartmann, Mathis Fricke, Lukas Weimar, Dirk Gr¨unding, Tomislav Mari´c, DieterBothe, and Steffen Hardt. Breakup dynamics of capillary bridges on hydrophobic stripes.
International Journal of Multiphase Flow , page 103582, 2021.[15] C. W. Hirt and B. D. Nichols. Volume of fluid/VOF/ method for the dynamics of free bound-aries.
J. Comput. Phys. , 39(1):201–225, 1981. ISSN 00219991. doi: 10.1016/0021-9991(81)90145-5. URL https://dx.doi.org/10.1016/0021-9991(81)90145-5 .[16] C. B. Ivey and P. Moin. Accurate interface normal and curvature estimates on three-dimensional unstructured non-convex polyhedral meshes.
J. Comput. Phys. , 300:365–386,2015. ISSN 10902716. doi: 10.1016/j.jcp.2015.07.055. URL http://dx.doi.org/10.1016/j.jcp.2015.07.055 .[17] H. Jasak.
Error analysis and estimation for the finite volume method with applications to fluidflows.
PhD thesis, 1996.[18] L. Jofre, O. Lehmkuhl, J. Castro, and A. Oliva. A 3-D Volume-of-Fluid advection methodbased on cell-vertex velocities for unstructured meshes.
Comput. Fluids , 94:14–29, 2014.ISSN 00457930. doi: 10.1016/j.compfluid.2014.02.001. URL http://dx.doi.org/10.1016/j.compfluid.2014.02.001 .[19] B. W.S. Jones, A. G. Malan, and N. A. Ilangakoon. The initialisation of volume fractionsfor unstructured grids using implicit surface definitions.
Comput. Fluids , 179:194–205, 2019.ISSN 00457930. doi: 10.1016/j.compfluid.2018.10.021. URL https://doi.org/10.1016/j.compfluid.2018.10.021 . 2720] F. Jureti´c. The cfMesh library for polyhedral mesh generation. https://sourceforge.net/projects/cfmesh/ . Accessed: 2020-01-15.[21] F. Jureti´c.
Error analysis in finite volume CFD . PhD thesis, Imperial College London (Uni-versity of London), 2005.[22] J. Kromer and D. Bothe. Highly accurate computation of volume fractions using differentialgeometry.
J. Comput. Phys. , 396(July):761–784, 2019. ISSN 00219991. doi: 10.1016/j.jcp.2019.07.005. URL https://dx.doi.org/10.1016/j.jcp.2019.07.005 .[23] J. L´opez, J. Hern´andez, P. G´omez, and F. Faura. Non-convex analytical and geometricaltools for volume truncation, initialization and conservation enforcement in VOF methods.
J.Comput. Phys. , 392:666–693, 2019. ISSN 10902716. doi: 10.1016/j.jcp.2019.04.055. URL https://doi.org/10.1016/j.jcp.2019.04.055 .[24] T. Mari´c, T. Tolle, and D. Gr¨unding. The argo OpenFOAM module: the implementa-tion of Surface Mesh Cell Intersection / Approximation algorithms. https://gitlab.com/leia-methods/argo/-/tree/2021-02-17-SMCIA-SUBMISSION . Accessed: 2021-02-16.[25] T. Mari´c, H. Marschall, and D. Bothe. lentFoam – A hybrid Level Set/Front Tracking methodon unstructured meshes.
Comput. Fluids , 113:20–31, may 2015. ISSN 00457930. doi: 10.1016/j.compfluid.2014.12.019. URL https://dx.doi.org/10.1016/j.compfluid.2014.12.019 .[26] T. Mari´c, H. Marschall, and D. Bothe. An enhanced un-split face-vertex flux-based VoFmethod.
J. Comput. Phys. , 371:967–993, apr 2018. ISSN 10902716. doi: 10.1016/j.jcp.2018.03.048. URL https://doi.org/10.1016/j.jcp.2018.03.048 .[27] T. Mari´c, D. B. Kothe, and D. Bothe. Unstructured un-split geometrical volume-of-fluidmethods–a review.
Journal of Computational Physics , 420:109695, 2020. doi: 10.1016/j.jcp.2020.109695. URL https://doi.org/10.1016/j.jcp.2020.109695 .[28] D. Meagher. Geometric modeling using octree encoding.
Comput. Graph. Image Process. , 19(2):129–147, 1982. ISSN 0146664X. doi: 10.1016/0146-664X(82)90104-6. URL http://dx.doi.org/10.1016/0146-664X(82)90104-6 .[29] D. P. Mehta and S. Sahni.
Handbook of data structures and applications . CRC Press, 2004.[30] S. Menon and D. P. Schmidt. Conservative interpolation on unstructured polyhedral meshes:An extension of the supermesh approach to cell-centered finite-volume variables.
Comput.Methods Appl. Mech. Eng. , 200(41-44):2797–2804, 2011. ISSN 00457825. doi: 10.1016/j.cma.2011.04.025. URL http://dx.doi.org/10.1016/j.cma.2011.04.025 .[31] F. Moukalled, L. Mangani, and M. Darwish.
The finite volume method in computationalfluid dynamics , volume 113. Springer, 2016. URL .[32] W. F. Noh and P. R. Woodward. SLIC (Simple Line Interface Calculation) method.
Proc.Fifth Int. Conf. Numer. Methods Fluid Dyn. June 28–July 2, 1976 Twente Univ. Enschede ,pages 330–340, 1976. doi: 10.1007/3-540-08004-X 336. URL https://dx.doi.org/10.1007/3-540-08004-X_336 . 2833] M. Owkes and O. Desjardins. A mesh-decoupled height function method for computing inter-face curvature.
J. Comput. Phys. , 281:285–300, 2015. ISSN 10902716. doi: 10.1016/j.jcp.2014.10.036. URL https://dx.doi.org/10.1016/j.jcp.2014.10.036 .[34] M. Owkes and O. Desjardins. A mass and momentum conserving unsplit semi-Lagrangianframework for simulating multiphase flows.
J. Comput. Phys. , 332:21–46, 2017. ISSN 10902716.URL http://dx.doi.org/10.1016/j.jcp.2016.11.046 .[35] W. J. Rider and D. B. Kothe. Reconstructing Volume Tracking.
J. Comput. Phys. , 141(2):112–152, 1998. ISSN 00219991. doi: 10.1006/jcph.1998.5906. URL https://doi.org/10.1006/jcph.1998.5906 .[36] G. Russo and P. Smereka. A Remark on Computing Distance Functions.
J. Comput. Phys. ,163(1):51–67, 2000. ISSN 00219991. doi: 10.1006/jcph.2000.6553. URL https://dx.doi.org/10.1006/jcph.2000.6553 .[37] H. Scheufler and J. Roenby. Accurate and efficient surface reconstruction from volume fractiondata on general meshes.
J. Comput. Phys. , 383:1–23, apr 2019. ISSN 10902716. doi: 10.1016/j.jcp.2019.01.009. URL https://doi.org/10.1016/j.jcp.2019.01.009 .[38] S. Shin and D. Juric. Modeling Three-Dimensional Multiphase Flow Using a Level ContourReconstruction Method for Front Tracking without Connectivity.
J. Comput. Phys. , 180(2):427–470, 2002. ISSN 00219991. doi: 10.1006/jcph.2002.7086. URL https://dx.doi.org/10.1006/jcph.2002.7086 .[39] S. Shin, I. Yoon, and D. Juric. The Local Front Reconstruction Method for direct simulationof two- and three-dimensional multiphase flows.
J. Comput. Phys. , 230(17):6605–6646, 2011.ISSN 00219991. doi: 10.1016/j.jcp.2011.04.040. URL http://dx.doi.org/10.1016/j.jcp.2011.04.040 .[40] S. Strobl, A. Formella, and T. P¨oschel. Exact calculation of the overlap volume of spheres andmesh elements.
J. Comput. Phys. , 311:158–172, 2016. ISSN 10902716. doi: 10.1016/j.jcp.2016.02.003. URL http://dx.doi.org/10.1016/j.jcp.2016.02.003 .[41] M. Sussman and E. Fatemi. An efficient, interface-preserving level set redistancing algo-rithm and its application to interfacial incompressible fluid flow.
SIAM Journal on sci-entific computing , 20(4):1165–1191, 1999. doi: 10.1137/S1064827596298245. URL https://doi.org/10.1137/S1064827596298245 .[42] M. Sussman, E. Fatemi, P. Smereka, and S. Osher. An improved level set method for in-compressible two-phase flows.
Computers & Fluids , 27(5-6):663–680, 1998. doi: 10.1016/S0045-7930(97)00053-4. URL https://doi.org/10.1016/S0045-7930(97)00053-4 .[43] M. Sussman, A. S. Almgren, J. B. Bell, P. Colella, L. H. Howell, and M. L. Welcome. Anadaptive level set approach for incompressible two-phase flows.
Journal of ComputationalPhysics , 148(1):81–124, 1999. doi: 10.1006/jcph.1998.6106. URL https://doi.org/10.1006/jcph.1998.6106 .[44] G. Th¨urrner and W¨uthrich C. A. Computing vertex normals from polygonal facets.
Journal ofGraphics Tools , 3(1):43–46, 1998. doi: 10.1080/10867651.1998.10487487. URL https://doi.org/10.1080/10867651.1998.10487487 . 2945] T. Tolle, D. Bothe, and T. Mari´c. SAAMPLE: A Segregated Accuracy-driven Algorithm forMultiphase Pressure-Linked Equations.
Comput. Fluids , 200:104450, 2020. ISSN 00457930. doi:10.1016/j.compfluid.2020.104450. URL https://dx.doi.org/10.1016/j.compfluid.2020.104450 .[46] Tobias Tolle, Dirk Gr¨unding, Dieter Bothe, and Tomislav Maric. Computing volume fractionsand signed distances from arbitrary surfaces on unstructured meshes: diagram data, 2021-02-18. URL https://tudatalib.ulb.tu-darmstadt.de/handle/tudatalib/2581.3 .[47] C. D. Toth, J. O’Rourke, and Jacob E. Goodman.
Handbook of discrete and computationalgeometry . Chapman and Hall/CRC, 2017.[48] G. Tryggvason, B. Bunner, A. Esmaeeli, D. Juric, N. Al-Rawahi, W. Tauber, J. Han, S. Nas,and Y. J Jan. A front-tracking method for the computations of multiphase flow.
J. Comput.Phys. , 169(2):708–759, 2001. ISSN 00219991. doi: 10.1006/jcph.2001.6726. URL https://doi.org/10.1006/jcph.2001.6726 .[49] O. Ubbink.