An Edge-based Interface Tracking (EBIT) Method for Multiphase-flows Simulation with Surface Tension

We present a novel Front-Tracking method, the Edge-Based Interface Tracking (EBIT) method for multiphase flow simulations. In the EBIT method, the markers are located on the grid edges and the interface can be reconstructed without storing the connectivity of the markers. This feature makes the process of marker addition or removal easier than in the traditional Front-Tracking method. The EBIT method also allows almost automatic parallelization due to the lack of explicit connectivity. In a previous journal article we have presented the kinematic part of the EBIT method, that includes the algorithms for interface linear reconstruction and advection. Here, we complete the presentation of the EBIT method and ∗Corresponding author Email addresses: yi.pan@sorbonne-universite.fr (Jieyun Pan), tian.long@dalembert.upmc.fr (Tian Long), leonardo.chirco@framatome.com (Leonardo Chirco), ruben.scardovelli@unibo.it (Ruben Scardovelli), popinet@basilisk.fr (Stéphane Popinet), stephane.zaleski@sorbonne-universite.fr (Stéphane Zaleski) Preprint submitted to Elsevier September 4, 2023 ar X iv :2 30 9. 00 33 8v 1 [ ph ys ic s. fl udy n] 1 S ep 2 02 3 combine the kinematic algorithm with a Navier–Stokes solver. To identify the reference phase and to distinguish ambiguous topological configurations, we introduce a new feature: the Color Vertex. For the coupling with the Navier–Stokes equations, we first calculate volume fractions from the position of the markers and the Color Vertex, then viscosity and density fields from the computed volume fractions and finally surface tension stresses with the Height-Function method. In addition, an automatic topology change algorithm is implemented into the EBIT method, making it possible the simulation of more complex flows. A two-dimensional version of the EBIT method has been implemented in the open-source Basilisk platform, and validated with five standard test cases: (1) translation with uniform velocity, (2) single vortex, (3) capillary wave, (4) Rayleigh-Taylor instability and (5) rising bubble. The results are compared with those obtained with the Volume-of-Fluid (VOF) method already implemented in Basilisk.


Introduction
Multiphase flows are ubiquitous in nature and engineering, and their numerical simulation still represents a formidable challenge, especially when a wide range of scales is involved, as in breaking waves on the sea surface, in some industrial processes or in atomizing liquid jets. Scales from meters to microns are typically seen. Large Reynolds number turbulence, as well as mass and heat transfers are the main causes for the introduction of such a wide range of scales. The CO 2 transfer and heat exchange between the oceans and the atmosphere is tightly linked to multiphase flows, as it takes place through the production and dispersion of small bubbles and droplets.
Similar small structures are observed in technology, for example in the heat exchange and transport in nuclear reactors, the atomization of liquids in combustion and other settings, and in most of the synthesis processes in chemical engineering. For many of these natural and engineering problems, numerical modelling is extremely desirable, albeit a monumental challenge.
However, in the general area of multiphase flow simulation, much progress has recently been achieved in the interrelated issues of multiple scales, discretization, and topology change. These issues are: i) the long-standing one of how to represent numerically (i.e., simulate) a dynamic or moving curve or surface, ii) the more difficult or challenging problem of how to efficiently model and discretize problems at multiple scales, and finally iii) the still open issue of how to take advantage of hierarchical data structures and grids, such as quadtrees and octrees, to address the "tyranny of small scales" [1].
The first, moving curve issue is divided into two problems: the kinematic problem, where the motion of the curve or surface separating the phases must be described by knowing the fluid velocity field and the rate of phase change, and the dynamic problem, where the momentum and energy conservation equations must be solved for given fluid properties. Methods available for the kinematic problem are sometimes separated into Front-Capturing and Front-Tracking methods. In the first kind, Front-Capturing methods, a tracer or marker function f (x, t) is integrated in time with the knowledge of an adequate velocity field u(x, t) The tracer function can be a Heaviside function, leading to Volume-of-Fluid (VOF) methods [2,3] or a smooth function, leading to Level-Set (LS) methods [4,5].
In Front-Tracking methods, the interface or "front" is represented by a curve discretization, for example a spline, which evolves with a prescribed normal velocity V S , see [6]. Compatibility between the two formulations is achieved if u · n = V S . An introduction to the most popular methods may be found in [7]. Another point of view on the methods, that is fruitful in connection with issues ii) and iii) above, can be obtained by investigating how the data structures are tied to the underlying Eulerian grid. The data structures are local when their components have little or no "knowledge" of the overall connections among pieces of interface or connected regions of a given phase (Fig. 1b). Thus a discretization of the tracer function f is a local data structure, tied to the Eulerian grid. On the other hand a global data structure contains information not only about individual points on the interface, but also about their connections with the entire interface of a given object. For Front-Tracking methods this is achieved by linked lists, or pointers, that allow to navigate the data structure along the object ( Fig. 1a). It is then obvious and efficient, using the data structure, to find all the connected pieces of the interface. The local data structures are not limited to VOF or LS methods. For example, an unconnected marker or particle method, such as Smoothed Particle Hydrodynamics (SPH), or other methods, simply tracks particles of position x i by integrating without any concept of connection between the particles [8]. Such an unconnected marker method, despite its similarity with Front-Tracking is in fact The local methods are easy to parallelize on a grid, are free of constraints that require the solutions of large linear or nonlinear systems (as when constructing interfaces by high-order splines) and are generally computationally efficient. The global methods allow to track and control the topology, so they are also useful for the second issue, the handling of multiscale problems. Moreover, when dealing with complex multiscale problems, it is often useful to investigate geometrical properties such as skeletons [9] which are branched manifolds most naturally represented by Front-Tracking. Global methods also allow to naturally distinguish between continuous slender objects such as unbroken thin ligaments or threads and strings of small particles or broken ligaments. This latter distinction is of great importance when an-alyzing statistically highly-fragmented flows [10]. It is however possible to represent slender objects with local methods in some cases: Chiodi [11] proposed an enhancement of the VOF interface reconstruction method, involving two nearly parallel, closely located planes in three dimensions, in order to capture slender sheets thinner than the mesh size.
A method that combines some of the properties of global methods, such as Front-Tracking, and local methods, such as VOF or LS, seems desirable.
There have indeed been some prior attempts at such a combination.
Aulisa et al. [12,13]  We suggest a similar method in the present work, based on a purely kinematic approach developed by two of us [20], called the Edge-Based Interface-Tracking (EBIT) method. In that method the position of the interface is tracked by marker points located on the edges of an Eulerian grid, and the connectivity information is implicit. The basic idea and the split interface advection were discussed in [20], here we improve the mass conservation of the original EBIT method by using a circle fit to reconstruct the interface during advection. The topology change mechanism is also introduced. Compared to the LFRM of Shin, Yoon and Juric [19], markers in the EBIT method are bound to the Eulerian grid. These marker are obtained by a reconstruction of the interface at every time step, thus the Eulerian grid and Lagrangian markers can be distributed to different processors by the same routine as the one used in the parallelization of the Navier-Stokes solver. Second, a new feature called Color Vertex, which amounts to describe the topology of the interface by the color of markers at the vertices of the grid, is discussed.
Such a scheme is trivial on a simplex and just slightly more complicated on a square grid. It was proposed in a similar context by Singh and Shyy [18]. Third, we combine the EBIT method with a Navier-Stokes solver for multiphase flow simulations. The coupling is simply realized by computing volume fractions from the position of the markers and the Color Vertex, then we use the volume fractions as in typical VOF solvers [7].
The paper is organized as follows: the kinematics of the EBIT method is described in Section 2. This includes the interface advection algorithm, the updating rules for the color vertices and the automatic topology change algorithm. Then the coupling algorithm between the EBIT interface description and the multiphase fluid dynamics is presented. A two-dimensional flow solver based on a Cartesian or quadtree grid is implemented inside the opensource platform Basilisk [21,22]. In Section 3, the EBIT method is validated by the computation of typical examples of multiphase flow simulations. The results obtained by the combined EBIT and VOF methods are presented and compared with those calculated by the pure VOF method.

The EBIT method
In the EBIT method, the interface is represented by a set of marker points placed on the grid lines. The advection of the interface is done by moving these points along the grid lines. The equation of motion for a marker point at position x i is which is discretized by a first-order explicit Euler method where the velocity u n i at the marker position x n i is calculated by a bilinear interpolation.
For a multi-dimensional problem, a split method is used to advect the interface (see Fig. 2), which is similar to that described in [20], but with some

Color Vertex
In the EBIT method, the connectivity of the markers is implicit. When there are only two markers on the cell boundary, the interface portion is given by the segment connecting these two points. However, when there are four markers on the boundary, there are two possible alternative configurations, as shown in Fig. 4. In order to select one of the two configurations without any ambiguity, we consider a technique called Color Vertex, which was first proposed by Singh and Shyy [18] to implement an automatic topology change.
The value of the Color Vertex indicates the fluid phase in the corresponding region within the cell, and five color vertices (four in the corners and one in the center of the cell) are enough to select one of the two configurations of Fig. 4. In other word, we can establish a one-to-one correspondence between the topological configuration and the value of the color vertices within each cell, and then reconstruct the interface segments with no ambiguity.
Furthermore, the direction of the unit normal to the interface can also be it is compared to the data structure that is used for storing the connectivity in traditional Front-Tracking methods.
As the interface is advected, the value of a Color Vertex should also be updated accordingly to ensure that the implicit connectivity information is corresponding segment through points 1' and 2', before the advection step.
These two points are on opposite sides in Fig. 6a and on consecutive sides in

Topology change
The topology change is controlled by the Color Vertex value distribution in an automatic manner. In the present implementation of the EBIT method, a marker point is present on a cell edge only if the value of the Color Vertex at the two edge endpoints is different.
Furthermore, only one marker is allowed on a cell edge. When two markers move into the same edge, as shown in Fig. 7, they both will be eliminated breakup merge Figure 7: Topology change mechanism automatically, because the value of the Color Vertex at the corresponding endpoints is the same, hence there cannot be any marker on that edge.
Moreover, the "surviving" markers within the cell will be reconnected automatically, see again Fig. 7. This reconnection procedure enables ligament breakup or droplet merging in an automatic way during the interface advection. As a direct consequence of this procedure, the volume occupied by the reference phase will decrease or increase. In particular, it tends to remove droplets or bubbles which are smaller than the grid size.
This topology change mechanism only affects interfaces that are approaching each other along a direction parallel to the grid lines. Because of the presence of a Color Vertex in the cell center, interfaces approaching along a diagonal direction do not induce any topology change, as long as the four markers remain on a different edge, as shown in Fig. 4. Thus, this mechanism does not suffer the problem of orientation-dependent topology change, which takes place during the tetra-marching procedure in Yoon's [23] LCRM and Shin's [19] LFRM, along the two diagonal directions.
However, in our method approaching interfaces along diagonal directions behave differently from those approaching along parallel ones. A possible solution to this problem would be removing the restriction on the number of markers per edge, which would also allow us to capture the subscale interfacial structure and to control the topology change based on a physical mechanism.

Governing equations
The Navier-Stokes equations for incompressible two-phase flow with immiscible fluids written in the one-fluid formulation are where ρ and µ are density and viscosity, respectively. The gravitational force is taken into account with the ρg term. Surface tension is modeled by the term f = σκnδ S , where σ is the surface tension coefficient, κ the interface curvature, n the unit normal and δ S the surface Dirac delta function.
The physical properties are calculated as where H is the Heaviside function, which is equal to 1 inside the reference phase and 0 elsewhere.
Since the markers are located on the grid edges, we consider a simple strategy to couple the EBIT method with the Navier-Stokes equations. From The numerical implementation of the EBIT method has been written in the Basilisk C language [21,22], which adopts a time-staggered approximate projection method to solve the incompressible Navier-Stokes equations on a Cartesian mesh or a quad/octree mesh. The Bell-Colella-Glaz (BCG) [25] second-order scheme is used to discretize the advection term, and a fully implicit scheme for the diffusion term. A well-balanced Continuous Surface Force (CSF) method is used to calculate the surface tension term [22,26].

Adaptive mesh refinement (AMR)
An efficient adaptive mesh refinement (AMR) technique, which is based on a wavelet decomposition [27] of the specified field variables, is implemented To avoid this inconsistency, we consider a simple strategy, that is to refine the cells within the 3 × 3 stencil of each interfacial cell to the maximum allowable level. With this assumption and the timestep limitation due to the CF L number, the interface will not be advected between two grid cells at different resolution levels. In principle, this refinement strategy should be less efficient than that based on curvature, however the numerical results of the next section show that the efficiency is still comparable to that based on curvature when other criteria of refinement, such as velocity gradients, are introduced into the AMR strategy.

Translation with uniform velocity
In this first test a circular interface of radius R = 0.15 and center at We use this case to test the new EBIT method, where the unaligned markers are computed by a circle fit. The accuracy and mass conservation of the method are measured by the area and shape errors. The area error E area is defined as the absolute value of the relative difference between the area occupied by the reference phase at the initial time t = 0 and that at The shape error, in a L ∞ norm, is defined as the maximum distance between any marker x i on the interface and the corresponding closest point on the analytical solution where (x c , y c ) are the coordinates of the circle center and R its radius.
To initialize the markers on the grid lines we first compute the signed distance (10) on the cell vertices, then we use a root-finding routine, when The area errors E area and the shape errors E shape are listed in Table 1 and are shown in Fig. 11 for the different methods. For the new EBIT method we observe a second-order convergence rate for the area error and approximately a first-order convergence rate for the shape error, as shown in Fig. 11. For the first version of the EBIT method (SL), both the area and shape errors are much larger. Furthermore, this method loses all the reference phase at the lowest grid resolution, N x = 32 (this fact is denoted by "NA" in Table 1. For the PLIC-VOF method, the mass conservation is accurate to machine error and it is not shown in the figure, while the shape error is evaluated with the two endpoints of the PLIC-VOF reconstruction in each cut cell. For this error we observe in Fig. 11a first-order convergence rate. Due to the fact that in the new EBIT method we are fitting a circle with a circle, the shape error obtained with the PLIC-VOF method is larger than that of the new EBIT method.

Single vortex
The single vortex test was designed to test the ability of an interface tracking method to follow the evolution in time of an interface that is highly stretched and deformed [28].    The area error E area and the shape error E shape are listed in Table 2 and are shown in Fig. 14 for the different methods here considered. For the new EBIT method that implements a circle fit, a convergence rate between first-order and second-order is observed for both errors. The shape errors  The interface line at maximum deformation and back to its initial position, for the test with period T = 8 is shown in Fig. 15 for different mesh resolutions. In this test, the interface line at maximum deformation at half- For the EBIT method with a straight line fit, there is an even more pronounced mass loss. Furthermore, the interface does not return to its  Fig. 16a at maximum deformation.
The area error E area and the shape error E shape are listed in Table 3 and are shown in Fig. 17 In order to demonstrate the feasibility of the integration of the new EBIT method with AMR, we run again the single-vortex test case on a quadtree grid in Basilisk. The interface lines at halftime and at the end of the simulation and the corresponding meshes are shown in Fig. 18. The maximum level of refinement is N l,max = 7, which corresponds to the mesh resolution N x = 128, while the minimum level is N l,min = 4. All the cells near the interface are refined to the maximum level, to avoid an inconsistency like that shown in Fig. 8. Since the imposed velocity field is not affected by the mesh resolution and the same stencil is used for the velocity interpolation, the interface line obtained on the quadtree grid coincides with that on the fixed Cartesian mesh. A sinusoidal perturbation is applied to a plane interface between two fluids initially at rest. Under the influence of surface tension, the interface begins to oscillate around its equilibrium position, while the amplitude of the oscillations decay in time due to viscous dissipation. The exact analytical solution was found by Prosperetti [29] in the limit of very small amplitudes, and it is usually used as a reference. and Gerlach et al. [32]. The value of the other physical properties that are used in the simulation and of the Laplace number, La = (ρ 1 σλ) µ 2 1 , is listed in Table 4.

Capillary wave
The time evolution of the maximum amplitude of the interface is shown in Fig. 19, together with those of the analytical solution and of the PLIC-VOF method. Time is made dimensionless by using the normal-mode oscillation frequency ω 0 , which is defined by the dispersion relation where k = 2π/λ is the wavenumber and with λ = 1 in our simulation. The numerical results obtained with the new EBIT and PLIC-VOF methods agree rather well with the analytical solution.
The error between the theoretical solution [29] and the two numerical solutions can be further analyzed with the L 2 norm where h is the maximum interface amplitude obtained with a numerical method and h exact the reference value. The results are shown in Fig. 20 and a convergence rate close to second-order is observed for both methods, the error of the PLIC-VOF method being always somewhat smaller. In both simulations, the height function method has been used to calculate the curvature [22]. More particularly, since we are considering a very small amplitude of the oscillations, thus of the interface curvature as well, the straight line approximation of the new EBIT method in each cut cell provides a fairly good approximation of the volume fraction and hence of the calculation of the local height function.

Rayleigh-Taylor instability
In order to demonstrate the capability of the new EBIT method to deal with more complex flows, we investigate another classical test: the Rayleigh-Taylor instability at high Reynolds number, that involves a large deformation of the interface.
The Rayleigh-Taylor instability occurs when a heavy fluid is on top of a lighter one, with the direction of gravity from top to bottom. The density difference between the two fluids plays an important role in the instability and is present in the definition of the dimensionless Atwood number At where ρ 1 and ρ 2 are the densities of the heavy and light fluids, respectively.   This instability has been investigated in several studies [33,34,35], that consider an incompressible flow without surface tension effects, with At = 0.5 and different Reynolds numbers, Re = ρ 1 g 1/2 d 3/2 µ 1 .
In this study, we consider the rectangular computational domain [0, d] × [0, 4d], partitioned with N x × N y = 128 × 512 grid cells. The plane interface y 0 (x) = 2 d between the two fluids is perturbed by a sinusoidal wave y 1 (x) = 0.1 d cos(kx), so that the interface line at the beginning of the simulation is with λ = d = 1. A no-slip boundary condition is enforced at the bottom and at the top of the computational domain, and a symmetric boundary condition on the two vertical sides. The value of the other physical properties that are used in the simulation and of the Reynolds number Re is listed in Table 5. The results of the simulation are first presented in Fig. 21 with the representation of the interface line at several dimensionless times τ = t g At/d.
Overall, these results compare rather well with those presented in [35].
More specifically, in the early stages of the simulation, τ ≤ 1.75, the shape Guermond [34] and Ding [35]. The interface lines obtained with the new EBIT method combined with AMR are finally shown in Fig. 23. The maximum level of refinement is now N l,max = 9, while the minimum level is N l,min = 6. The refinement criteria are based not only on the position of the interface, but also on the velocity gradient, thus the cells within the roll-up region, characterized by a strong vorticity, are also refined to the maximum level, even if they are not very close to the interface (see Fig. 23b). The interface lines calculated with the new EBIT method on the quadtree grid agree rather well with those on the fixed Cartesian grid (see Fig. 23a), only minor differences are observed in the roll-up region.
The computational efficiency for this test of the new EBIT method with or without AMR is summarized in Table 6. For the PLIC-VOF method with AMR (VOF-AMR), the mesh refinement criteria are based on both curvature and velocity gradient. Without AMR, the dynamics with the new EBIT method is about 2 times slower than with the PLIC-VOF method. But when AMR is used, the wall times for the two methods are comparable.

Rising bubble
In this test we examine a single bubble rising under buoyancy inside a heavier fluid. This test case was first proposed by Hysing [36] and it provides a standard benchmark for multiphase flow simulations, since this configuration is simple enough to be simulated accurately. Nevertheless, the bubble shows a strong deformation and even complex topology changes in some flow regimes [37,38], thus giving an adequate challenge to interface tracking techniques. properties is that provided by Hysing [36] and is listed in Table 7, where Re = ρ 1 g 1/2 d 3/2 µ 1 is the Reynolds number and Bo = ρ 1 gL 2 σ the Bond number.  Hysing [36].
In the second test, the bubble lies somewhere between the skirted and the The computational efficiency for these two cases of the new EBIT method with or without AMR is summarized in Table 8. Without AMR, the dynamics with the new EBIT method is about 2 times slower than with the PLIC-VOF method, for test case 1 and the Cartesian problem. For test case

Conclusions
We present a novel Front-Tracking method, the Edge-Based Interface Tracking (EBIT) method, which is more suitable for parallelization due to the lack of explicit connectivity. Several new features have been introduced to improve the very first version of the EBIT method [20]. First, a circle fit has been implemented to improve the accuracy of mass conservation after the interface advection in the reconstruction phase. Second, a Color Vertex feature has been introduced to distinguish between ambiguous topological configurations and to represent the connectivity implicitly. Third, an automatic topological change mechanism has been discussed. Numerical results for various cases, including both kinematic and dynamical tests, have been considered and compared with those obtained with the VOF method. Good agreement is observed for all test cases.
In future work, we aim to remove the restriction on the number of markers on a cell side, to improve our control on topological changes, and to extend the EBIT method to three dimensions.