 # THINC-scaling scheme that unifies VOF and level set methods

We present a novel interface-capturing scheme, THINC-scaling, to unify the VOF (volume of fluid) and the level set methods, which have been developed as two completely different approaches widely used in various applications. The THINC-scaling scheme preserves at the samectime the advantages of both VOF and level set methods, i.e. the mass/volume conservation of the VOF method and the geometrical faithfulness of the level set method. THINC-scaling scheme allows to represent interface with high-order polynomials, and has algorithmic simplicity which eases its implementation in unstructured grids.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

VOF (volume of fluid) and level set are the two most popularly used methods in capturing moving interface, and find their applications in diverse fields, such as numerical simulation of multiphase fluid dynamics, graphic processing, topological optimization and many others. Historically, they were independently developed based on completely different concepts and solution methodologies, and have their own superiority and weakness.

VOF methodHirt and Nichols (1981); Youngs (1982); Lafaurie et al. (1994); Rider and Kothe (1998); Scardovelli and Zaleski (2000) uses the volume fraction of one out of multiple species in a control volume (mesh cell) to describe the distribution of the targeted fluid in space. The VOF function by definition has a value between 0 and 1, and the interfaces can be identified as isosurfaces (3D) or contours (2D) of a fractional VOF value, say 0.5 for example. A more accurate way to represent the interface is using geometrical reconstruction, such as the PLIC (Piecewise Linear Interface Calculation) schemes which are currently accepted as the main-stream VOF methodology. Rigorous numerical conservation can be guaranteed if a finite volume method is used to transport the VOF function, which is found to be crucial in many applications. The VOF function is usually characterized by large jump or steep gradient across the interface. So, directly using the VOF function to retrieve the geometrical information of the interface, such as normal and curvature, may result in large error.

The level set methodOsher and Sethian (1988); Sethian (1999); Osher and Fedkiw (2003), on other hand, defines the field function as a signed distance function (or level set function) to the interface, which possesses a uniform gradient over the whole computational domain and thus provides a perfect field function to retrieve the geometrical properties of an interface. However, the level function is not conceptually nor algorithmically conservative. The numerical solution procedure, including both transport and reinitialization, does not guarantee the conservativeness in numerical solution. It may become a fatal problem in many applications, like multiphase flows involving bubbles or droplets.

Efforts have been made to combine the VOF method and level set method, which lead to the coupled level set/VOF methods (CLSVOF) Sussman and Puckett (2000); Ménard et al. (2007); Yang et al. (2006); Sun and Tao (2010). The PLIC type VOF method is blended with the level set method so as to improve both conservativeness and geometrical faithfulness in numerical solution.

In this paper, we propose a new scheme that unifies the VOF and level set methods, based on the observation that the VOF field can be seen as a scaled level set field using the THINC (Tangent of Hyperbola Interface Capturing) function, which has been used in a class of schemes for capturing moving interface Xiao et al. (2005, 2011); Ii et al. (2012); Xie and Xiao (2017); Qian and Xiao (2018). The resulting scheme, so-called THINC-scaling scheme, converts the field function from level set to VOF by the THINC function, and converts the VOF field back to the corresponding level set field via an inverse THINC function. So, it can make use of the advantages of both VOF and level set at different stages of solution procedure, which eventually realizes the high-fidelity computation of moving interface regarding both numerical conservativeness and geometrical representation.

## 2 The connection between level set and VOF functions

We consider an interface separating two kinds of fluids, fluid 1 and fluid 2, occupying volumes and respectively in space. We introduce the following two indicator functions to identify the different fluids and the interface. Figure 1: Two indicator functions to identify multi-materials and interface: (a) The level set function; (b) The THINC function (a VOF function).
• Level set function (Fig.1(a)):
The level set function is defined as a signed distance from a point in three dimensions to the interface by

 ϕ(x)=⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩infxI∈∂Ω∥x−xI∥if x∈Ω10if x∈∂Ω−infxI∈∂Ω∥x−xI∥if x∈Ω2, (1)

where represents any point on the interface, and referred to as interface point.

• VOF function (Fig.1(b)):
The VOF function in the limit of an infinitesimal control volume is the Heaviside function. Assuming the VOF function as the abundance of fluid 1, we have the Heaviside function in its canonical form as

 Hc(x)=⎧⎪ ⎪⎨⎪ ⎪⎩1if x∈Ω112if x∈∂Ω0if x∈Ω2. (2)

Recall that

 Hc(x)=12limβ→∞(1+tanh(βx))), (3)

we get a continuous Heaviside function

 H(x)=12(1+tanh(βx))) (4)

with a finite steepness parameter .

Given a computational mesh with cells of finite size, the cell-wise VOF function is defined by

 ¯HΩl=1|Ωl|∫ΩlH(x)dΩ, (5)

where is the target cell with a volume . In practice, we recognize the interface cell where an interface cuts through, in terms of the VOF velue, by with being a small positive, e.g. .

From the definition of level set function (1), we know that (4) scales a level set function to a Heaviside function.

Now, we establish a connection between the level set function and the Heaviside function.

• THINC scaling (:

We convert the level set field by

 H(x,t)=12(1+tanh(β(P(x% )))), (6)

where is a polynomial

 P(x)=p∑r,s,t=0arstxryszt (7)

whose coefficients are determined from the level set function through the following constraints.

 ∂DP(x)∂xdx∂ydy∂zdz=∂Dϕ(x)∂xdx∂ydy∂zdz,  (dx,dy,dz)=0,1,2⋯ and  dx+dy+dz=D. (8)

In practice, we calculate the coefficients of via numerical approximations using the discrete level set field available in the computational domain.

We refer to (6) as the THINC scaling formula, and (7) as the level set polynomial.

• Inverse THINC scaling (:

Given the THINC function, we can directly compute the corresponding level set function by

 ϕ(x)=1βtanh−1(2H(% x)−1)  or  ϕ(x)=1βln(H(x)(1−H(x)). (9)

Formula (9) is referred to as the inverse THINC scaling that converts the Heaviside function to the level set function.

Remark 1. Formulae (6) and (9) provide analytical relations to uniquely convert between the level set and the THINC functions, which unifies the two under a single framework handleable with conventional mathematical analysis tool, and more importantly allows us to build interface-capturing schemes that take advantages from both VOF and level set methods.

Remark 2. The interface is defined by

 P(x)=0, (10)

where the level set polynomial defined in (7) enables to accurately formulate the geometry of the interface. In principle, we can use arbitrarily high order surface polynomial to represent the interface straightforwardly without substantial difficulty.

Remark 3. A finite value of the steepness parameter modifies the Heaviside function to a continuous and differentiable function (6), which serves an adequate approximation to the VOF function with adequately large steepness parameter .

## 3 The THINC-scaling scheme for moving interface capturing

We assume that the moving interface is transported by a velocity field u. Thus, the two indicator functions are advected in the Eulerian form by the following equations, i.e.

 ∂H∂t+∇⋅(uH)=H∇⋅u (11)

for the THINC function , and

 ∂ϕ∂t+u⋅∇ϕ=0 (12)

for the level set function .

Next, we present the THINC-scaling scheme to simultaneously solve (11) and (12). The computational domain is composed of non-overlapped discrete grid cells of the volume , which can be either structured or unstructured grids. For any target cell element , we denote its mass center by , and its surface segments of areas by with . The outward unit normal is denoted by .

Assume that we know at time step () the VOF value

 ¯Hni=1|Ωi|∫ΩiH(x,tn)dΩ, (13)

for each cell, and the level set value

 ϕni=ϕ(xic,yic,zic,tn) (14)

at each cell center, we use the third-order TVD Runge-Kutta schemeShu (1988) for time integration to update both VOF and level set values, and , to the next time step ().

We hereby summarize the solution procedure of the THINC-scaling scheme for one Runge-Kutta sub-step that advance and at sub-step to and at sub-step .

Step 1. Compute the level set polynomial of th order for cell from the level set field using the constraint condition (8),

 Pi(x)=p∑r,s,t=0arstXrYsZt (15)

where is the local coordinates with respect to the center of cell , i.e. , , . The coefficients

are computed by Lagrange interpolation or least square method using the level set values

in the target and nearby cells.

Step 2: Construct the cell-wise THINC function under the constraint of VOF value by

 1|Ωi|∫Ωi12(1+tanh(β(Pi(x)+ϕΔi)))dΩ=¯Hmi. (16)

With a pre-specified and the surface polynomial obtained at step 1, the only unknown can be computed from (16). In practice, we use the numerical quadrature detailed in Xie and Xiao (2017), and the resulting nonlinear algebraic function of is solved by the Newton iterative method. See Xie and Xiao (2017) for details.

We then get the THINC function

 Hmi(x)=12(1+tanh(β(Pi(x)+ϕΔi))), (17)

which satisfies the conservation constraint of the VOF value, and is a correction to the interface location due to the conservation. The piece-wise interface in each interface cell is defined by

 ψi(x)≡Pi(x)+ϕΔi=0. (18)

We refer to (18) as the Polynomial Surface of the Interface (PSI) in cell .

Step 3: Update the VOF function by solving (11) through the following finite volume formulation,

 ¯Hm+1i=¯Hmi−Δt|Ωi|J∑j=1(∫Γij((u⋅n)Hmi(x)iup)dΓ)+¯Hmi|Ωi|J∑j=1(u⋅n)ij∣∣Γij∣∣)Δt, (19)

where the upwinding index is determined by

 iup={=i, for (u⋅n)ij>0;=ij, otherwise, (20)

and denotes the index of the neighboring cell that shares cell boundary with target cell . The integration on cell surface is computed by Gaussian quadrature formula.

Step 4: Update level set value at cell center using a semi-Lagrangian method as follows.

Step 4.1: We first find the departure point for each cell center by solving the initial value problem,

 ⎧⎪⎨⎪⎩dxdτ=−u(x,tn+τ)x(0)=xic (21)

up to , which leads to

. We use a second-order Runge-Kutta method to solve the ordinary differential equation (

21).

Step 4.2: Update the level set value at cell center using the Lagrangian invariant solution,

 ϕm+1i=~ϕmid(xid), (22)

where is the level function on cell where the departure point falls in. Using the inverse THINC-scaling formula (9), we immediately get the level set function from (17),

 ~ϕmid(x)=1βtanh−1(2Hmid(x)−1), (23)

which gives the level set value everywhere in cell that includes the departure point .

Step 5: Reinitialize the level set field. We fix the level set values computed from (23) for the interface cells which are identified by with and being small positive numbers. The level set values at the cell centers away from the interface region are reinitialized to satisfy the Eikonal equation,

 |∇ϕ|=1. (24)

We use the Fast Sweeping Method (FSM) Zhao (2005) on structured grid and the iteration method Dianat et al. (2017) on unstructured grid for reinitializing level set values.

Step 6: Go back to Step 1 for next sub-time step calculations.

Remark 4. The THINC-scaling scheme shown above unifies the VOF method and level set method. Eq.(17) retrieves the VOF field from the level set field, while (23) retrieves the level set field from the VOF field with numerical conservativeness.

Remark 5. The inverse THINC-scaling (23) facilitates a semi-Lagrangian solution without any spatial reconstruction or interpolation, such as those used in Strain (1999). This step essentially distinguishes the present scheme from the coupled THINC/level set method in Qian and Xiao (2018), where the level set function is updated by a fifth-order Hamilton-Jacobi WENO scheme with a 3rd-order TVD Runge-Kutta time-integration scheme.

Remark 6. The interface is cell-wisely the PSI defined by (18), , within the interface cells, which provides the sub-cell interface structure with geometrical information, such as position, normal direction and curvature to facilitate the computation of so-called sharp-interface formulation. It distinguishes the present scheme with superiority from any other algebraic interface-capturing methods.

Remark 7. As discussed in Xiao et al. (2011), the steepness parameter

can be estimated by the thickness of the jump transition across the interface using

 β=1ηtanh−1(1−2ε) (25)

where denotes the normalized half thickness of the jump with respect to the cell size, and is a small positive number to define the range of interface transition layer in terms of the VOF value, i.e. , we use in this work. In order to keep a 3-cell thickness for the interface, can be set as , which approximately results in . Our numerical experiments show that the THINC method can resolve sharply interfaces within one or two cells by using larger .

## 4 Numerical tests

We verify the THINC-scaling scheme to capture moving interfaces using some advection benchmark tests. We focus on the numerical tests in two dimensions on both structured (Cartesian) grid and unstructured (triangular) grid. The numerical errors in terms of the VOF field are quantified via the error norm .

 E(L1)=∑i|¯Hi−¯Hei||Ωi|, (26)

where and are respectively the numerical and exact VOF values.

### 4.1 Solid body rotation test

In this test, so-called Zalesak’s slotted disk test Zalesak (1979), initially a circle with a radius of 0.5 centered at in a unit square computational domain is notched with a slot defined by . The slotted circle is rotated with the velocity field given by .

The steepness parameter is set , and a quadratic polynomial

 Pi(x)=2∑r,s=0arsXrYs (27)

is used in this test. Figure 2: VOF 0.5 contour line in Zalesak solid body rotation test after one revolution on (a) 50×50, (b) 100×100 and (c) 200×200 meshes. The black dashed line stands for the exact solution, and the red solid line for the numerical solution. Figure 3: The numerical results of Zalesak solid body rotation test after one revolution showing the reconstructed interfaces on (a) 50×50, (b) 100×100 and (c) 200×200 meshes. The PSI (red solid line) of the interface cells are plotted against the exact solution (black dashed line).

We compute this test on structured grid for different grid sizes with 50, 100 and 200 vertices evenly distributed on each edge of computational domain. Fig.2 shows the numerical results of interface identified by VOF 0.5 contour line. As observed from these results, the interface in the slot region is well captured. As demonstrated in Xie and Xiao (2017) and Qian and Xiao (2018), the quadratic polynomial representation of the interface preserves the geometrical symmetry of the solution, which deteriorates significantly if a linear function (straight line) is used as commonly observed in the results of VOF method using PLIC reconstructions. We also show the PSI of the interface cells in Fig. 3. The interface is retrieved and represented by the cell-wise quadratic curves in the interface cells.

### 4.2 Vortex deformation transport test

The THINC-scaling scheme is further assessed by single vortex test Rider and Kothe (1998), in which a circle initially centred at in a unit square domain is advected by time dependent velocity field given by the stream function as follows,

 Ψ(x,y,t)=1πsin2(πx)sin2(πy)cos(πtT), (28)

where is specified in this test. This test, as one of the most widely used benchmark tests, is more challenging to assess the capability of the scheme in capturing the heavily distorted interface with stretched tail when transported to . From to , the reverse velocity field restores the interface back to its initial shape. In case of , the spiral tail becomes so thin that it can not be resolved by the resolution of a coarse grid. Figure 4: Numerical results for Rider-Kothe single vortex test showing the PSI for interface cells at t=T/2 and the VOF 0.5 contour lines at t=T on (a) 64×64, (b) 128×128 and (c) 256×256 meshes Figure 5: The distorted interface at T/2 on 128×128 grid (a) where the boxed part is enlarged as the close-up of the thin tail in panel (b) which shows that the film of tail is under the resolution of grid cell, but can be still resolved by the THINC-scaling scheme.

We tested this case on Cartesian grids with different resolutions of , and respectively. Numerical results on different grids at and are shown in Fig. 4. The THINC-scaling scheme can capture the elongated tail even when the interface is under grid resolution, and can restore the initial circle with good solution quality.

The PSI at on grid is shown in Fig. 5. The tail tip is stretched into a thin film with a thickness smaller than the cell size. This sub-cell structure can still be reconstructed by the THINC-scaling scheme with quadratic or higher order polynomial representation. Consequently, the pieces of flotsam generated from the PLIC VOF methods are not observed here.

A quantitative comparison with other sophisticated geometrical VOF methods Maric et al. (2018); Owkes and Desjardins (2014); Scheufler and Roenby (2018) is given in Table LABEL:r-k-comparison. It reveals the appealing accuracy of the present scheme.

We solved the Rider-Kothe shear flow test case on an unstructured grid with triangular cells. In order to compare with the results in Scheufler and Roenby (2018), we set nodes on the domain boundaries in and directions respectively. We used , where represents the cell size and is defined as the hydraulic diameter with and being the area and perimeter of the triangular cell element respectively. Figure 6: Numerical results for Rider-Kothe single vortex test on a triangular unstructured grid with 64 nodes on each boundary of computational domain. (a): VOF 0.5 contour line at t=T and PSI of interface cells at t=T/2; (b): PSI of interface cells at t=T/2 with the stretched tail highlighted by the black dash-line box; (c): Enlarged view of PSI in the interface cells highlighted in (b).

We plot the numerical results in Fig. 6. The PSIs in the interface cells accurately present the reconstructed interface. The flotsams in the PLIC reconstruction Scheufler and Roenby (2018) are avoided in the present results. A quantitative comparison is given in Table LABEL:r-k-unstructured-comparison. The THINC-scaling scheme shows superiority in accuracy on unstructured grid.

## 5 Conclusions

We propose a novel scheme to unify the solution procedures of level set and VOF methods that are two interface-capturing methods based on completely different concept and numerical methodology. The underlying idea of the scheme, THINC-scaling scheme, is to use the THINC function to scale/convert between the level set function and a continuous Heaviside function which mimics the VOF field.

The THINC function uses the level set polynomial to accurately retrieve the geometrical information of the interface from the level set field with high-order polynomials, and is constructed under the constraint of VOF value to rigorously satisfy the numerical conservativeness. Being a function handleable with conventional calculus tools, the THINC function facilitates efficient and accurate computations in interface-capturing schemes, which take the advantages from both VOF and level set methods.

We verified the THINC-scaling scheme with advection benchmark tests for moving interfaces on both unstructured and unstructured grids, which demonstrate the super solution quality and the great potential of the proposed scheme as a moving interface-capturing scheme for practical utility.

## 6 Acknowledgments

This work was supported in part by the fund from JSPS (Japan Society for the Promotion of Science) under Grant No. 18H01366.

## References

• M. Dianat, M. Skarysz, and A. Garmory (2017) A coupled level set and volume of fluid method for automotive exterior water management applications. International Journal of Multiphase Flow 91, pp. 19–38. Cited by: item.
• C. W. Hirt and B. D. Nichols (1981) Volume of fluid (vof) method for the dynamics of free boundaries. Journal of computational physics 39 (1), pp. 201–225. Cited by: §1.
• S. Ii, K. Sugiyama, S. Takeuchi, S. Takagi, Y. Matsumoto, and F. Xiao (2012) An interface capturing method with a continuous function: the thinc method with multi-dimensional reconstruction. Journal of Computational Physics 231 (5), pp. 2328–2358. Cited by: §1.
• B. Lafaurie, C. Nardone, R. Scardovelli, S. Zaleski, and G. Zanetti (1994) Modelling merging and fragmentation in multiphase flows with surfer. Journal of Computational Physics 113 (1), pp. 134–147. Cited by: §1.
• T. Maric, H. Marschall, and D. Bothe (2018) An enhanced un-split face-vertex flux-based vof method. Journal of computational physics 371, pp. 967–993. Cited by: §4.2, Table 1.
• T. Ménard, S. Tanguy, and A. Berlemont (2007) Coupling level set/vof/ghost fluid methods: validation and application to 3d simulation of the primary break-up of a liquid jet. International Journal of Multiphase Flow 33 (5), pp. 510–524. Cited by: §1.
• S. Osher and R. Fedkiw (2003) Implicit functions. In Level Set Methods and Dynamic Implicit Surfaces, pp. 3–16. Cited by: §1.
• S. Osher and J. A. Sethian (1988) Fronts propagating with curvature-dependent speed: algorithms based on hamilton-jacobi formulations. Journal of computational physics 79 (1), pp. 12–49. Cited by: §1.
• M. Owkes and O. Desjardins (2014) A computational framework for conservative, three-dimensional, unsplit, geometric transport with application to the volume-of-fluid (vof) method. Journal of computational physics 270, pp. 587–612. Cited by: §4.2, Table 1.
• L. Qian and F. Xiao (2018) Coupled thinc and level set method: a conservative interface capturing scheme with high-order surface representations. Journal of Computational Physics 373, pp. 284–303. Cited by: §1, item, §4.1.
• W. J. Rider and D. B. Kothe (1998) Reconstructing volume tracking. Journal of computational physics 141 (2), pp. 112–152. Cited by: §1, §4.2.
• R. Scardovelli and S. Zaleski (2000) Analytical relations connecting linear interfaces and volume fractions in rectangular grids. Journal of Computational Physics 164 (1), pp. 228–237. Cited by: §1.
• H. Scheufler and J. Roenby (2018) Accurate and efficient surface reconstruction from volume fraction data on general meshes. Journal of computational physics 383, pp. 1–23. Cited by: §4.2, §4.2, §4.2, Table 1, Table 2.
• J. A. Sethian (1999)

Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science

.
Vol. 3, Cambridge university press. Cited by: §1.
• C. Shu (1988) Total-variation-diminishing time discretizations. SIAM J. Sci. Stat. Comput. 9, pp. 1073–1084. Cited by: §3.
• J. Strain (1999) Semi-lagrange methods for level set equations. Journal of Computational Physics 2, pp. 498–533. Cited by: item.
• D. Sun and W. Tao (2010) A coupled volume-of-fluid and level set (voset) method for computing incompressible two-phase flows. International Journal of Heat and Mass Transfer 53 (4), pp. 645–655. Cited by: §1.
• M. Sussman and E. G. Puckett (2000) A coupled level set and volume-of-fluid method for computing 3d and axisymmetric incompressible two-phase flows. Journal of computational physics 162 (2), pp. 301–337. Cited by: §1.
• F. Xiao, Y. Honma, and T. Kono (2005) A simple algebraic interface capturing scheme using hyperbolic tangent function. International Journal for Numerical Methods in Fluids 48 (9), pp. 1023–1040. Cited by: §1.
• F. Xiao, S. Ii, and C. Chen (2011) Revisit to the thinc scheme: a simple algebraic vof algorithm. Journal of Computational Physics 230 (19), pp. 7086–7092. Cited by: §1, item.
• B. Xie and F. Xiao (2017) Toward efficient and accurate interface capturing on arbitrary hybrid unstructured grids: the thinc method with quadratic surface representation and gaussian quadrature. Journal of Computational Physics 349, pp. 415–440. Cited by: §1, item, §4.1, Table 2.
• X. Yang, A. J. James, J. Lowengrub, X. Zheng, and V. Cristini (2006) An adaptive coupled level-set/volume-of-fluid interface capturing method for unstructured triangular grids. Journal of Computational Physics 217 (2), pp. 364–394. Cited by: §1.
• D. L. Youngs (1982) Time-dependent multi-material flow with large fluid distortion. Numerical methods for fluid dynamics. Cited by: §1.
• S. T. Zalesak (1979) Fully multidimensional flux-corrected transport algorithms for fluids. Journal of computational physics 31 (3), pp. 335–362. Cited by: §4.1.
• H. Zhao (2005) A fast sweeping method for eikonal equations. Mathematics of computation 74 (250), pp. 603–627. Cited by: item.