An Improved Conservative Direct Re-Initialization Method (ICDR) for Two-Phase Flow Simulations

We introduce an improved conservative direct re-initialization (ICDR) method (for twophase flow problems) as a new and efficient geometrical re-distancing scheme. The ICDR technique takes advantage of two mass-conserving and fast re-distancing schemes, as well as a global mass correction concept to reduce the extent of the mass loss/gain in twoand three-dimensional (2D and 3D) problems. We examine the ICDR method, at the first step, with two 2D benchmarks: the notched cylinder and the swirling flow vortex problems. To do so, we (for the first time) extensively analyze the dependency of the regenerated interface quality on both time-step and element sizes. Then, we quantitatively assess the results by employing a defined norm value, which evaluates the deviation from the exact solution. We also present a visual assessment by graphical demonstration of original and regenerated interfaces. In the next step, we investigate the performance of the ICDR in three-dimensional (3D) problems. For this purpose, we simulate drop deformation in a simple shear flow field. We describe our reason for this choice and show that, by employing the ICDR scheme, the results of our analysis comply with the existing numerical and experimental data in the literature.


Introduction
The level set (LS) method has been a popular and efficient means to capture interfaces in different two-phase (or multi-phase) flow fields since being introduced by Osher and Sethian [1]. In this method, the signed distance function (SDF) is approximated using the level set function or value (Φ). Therefore, one can employ the LS values at any arbitrary time (t) to explicitly estimate the interface without applying rigorous reconstruction algorithms, as is the case for the volume of fluid (VOF) method [2][3][4][5][6]. Besides, the evaluation of the interface normal vectors and curvature values are easily attainable, owing to the smoothness of the signed distance function [7].
The level set values equal the signed distance function at the initiation time (t = 0) when the interface position is pre-defined. However, after some time steps (∆t) of the interface displacement, the LS function deviates from the SDF. Consequently, it is necessary to re-initialize the LS values to keep them as the signed distance function.
The re-initialization does not violate the overall mass conservation law since both phases' mass summation is constant. However, due to the interface displacement after the re-initialization, the area or the volume occupied by one phase may reduce. On the contrary, the content of the other phase increases. The unphysical mass loss/gain after the re-initialization step is a numerical error.
Sussman et al. [8], for the first time, proposed a re-initialization method. They employed the Eikonal relation to force the Euclidean norm of the LS gradient to be one via numerical and experimental data in the literature. Finally, in Section 5, the conclusion is provided.

Interface Advection Algorithm
The level set function (Φ) in the LS method approximates the signed distance function. Therefore, a curved line or surface with the property of Φ = 0 defines the interface. We arbitrarily assign the primary phase by Φ > 0, and Φ < 0 holds for the second phase. Since the LS values are essential in the vicinity of the interface, it is desired to employ a narrowband LS method introduced by Chopp [28]. This technique limits the LS calculations to the interface vicinity instead of the whole domain. Figure 1, in a structured mesh N 1 × N 2 , shows a typical interface (the solid line) and the elements in its vicinity. The elements in the vicinity of the interface are defined using the interface radius (R i ), which is twice or three times the greatest element diameter [7].
Fluids 2021, 6, x FOR PEER REVIEW 3 of 25 unstructured meshes. In this section, we have also employed the ICDR method to study drop deformation in a simple shear flow field and have compared our results with the numerical and experimental data in the literature. Finally, in Section 5, the conclusion is provided.

Interface Advection Algorithm
The level set function ( Φ ) in the LS method approximates the signed distance function. Therefore, a curved line or surface with the property of Φ 0 defines the interface. We arbitrarily assign the primary phase by Φ 0, and Φ 0 holds for the second phase. Since the LS values are essential in the vicinity of the interface, it is desired to employ a narrow-band LS method introduced by Chopp [28]. This technique limits the LS calculations to the interface vicinity instead of the whole domain. Figure 1, in a structured mesh , shows a typical interface (the solid line) and the elements in its vicinity. The elements in the vicinity of the interface are defined using the interface radius ( ), which is twice or three times the greatest element diameter [7]. In Figure 1, the (light and dark) gray elements demonstrate the vicinity of the interface (the thick black line). Moreover, the dark gray shows active elements of the grids, which contain an interface. Here, we normalize the signed distance function as stated in Equation (1), and name it nodal volume fraction.
In Equation (1), the parameters and show the signed distance function and the nodal volume fraction, respectively. On this basis, the interface advection equation would be as follows: In Equation (2), which is based on Einstein summation notations over a common index, 1, 2, 3 are the components of the velocity vector ( ⃗ , , ). The parameter 1, 2, 3 show the coordinates of any arbitrary point in a Cartesian reference frame. Several stable strategies to solve the interface advection equation are proposed in the literature: Streamline Upwind Petrov-Galerkin [24] and Taylor-Galerkin [26] methods. In this research, we have employed a least-squares finite element method [7,27]. For this In Figure 1, the (light and dark) gray elements demonstrate the vicinity of the interface (the thick black line). Moreover, the dark gray shows active elements of the grids, which contain an interface. Here, we normalize the signed distance function as stated in Equation (1), and name it nodal volume fraction.
In Equation (1), the parameters d and f show the signed distance function and the nodal volume fraction, respectively. On this basis, the interface advection equation would be as follows: In Equation (2), which is based on Einstein summation notations over a common index, u i (i = 1, 2, and 3) are the components of the velocity vector ( → u = (u 1 , u 2 , u 3 )). The parameter x i (i = 1, 2, and 3) show the coordinates of any arbitrary point in a Cartesian reference frame.
Several stable strategies to solve the interface advection equation are proposed in the literature: Streamline Upwind Petrov-Galerkin [24] and Taylor-Galerkin [26] methods. In this research, we have employed a least-squares finite element method [7,27]. For this purpose, the domain is divided into rectangular or hexahedral subdomains or elements for 2D and 3D cases, respectively. Then, a linear-based shape function (φ) is considered to evaluate the volume fraction ( f ) in each element. Therefore, Equation (2) is discretized as follows: where { f } and φ are the column-wise and row-wise vectors of the calculated volume fractions and shape function coefficients in each element, respectively. To solve Equation (3), the explicit values of the velocity components (u i ) from the current time step are used. Additionally, the Crank-Nicolson method is employed to evaluate the nodal values of the volume fraction, as stated in Equation (4).
In this equation, the parameters f (0) and f (1) are the nodal volume fraction values at the current and next steps, respectively. The term α is the Crank-Nicolson adjusting parameters and changes between 0 and 1. α = 1 and α = 0 turn the equation to fully implicit and explicit methods, correspondingly. Substituting the f from Equation (4) The parameter ∆t in Equation (5b,c) is the time-step size. The weighted residual form of Equation (5a) is: In Equation (6), ω is the weight function, which is equal to Φ 1 in our least-square finite element method [29]. Equation (6) results in a linear set of equations for a known velocity field.
The time step (∆t) is an essential parameter in the solution of Equation (6). To find a suitable time increment, we define a characteristic time for each element. For this purpose, we assume an arbitrary 3D element, as shown in Figure 2. In this element, three different vectors → ab, → cd and → e f can easily be defined by connecting the middle of the front faces. If the velocity at the center of the element is → u o = (u o1 , u o2 , u o3 ), the projection of the velocity vector on each vector would be:  In Equation (7a,b), the operator "." demonstrates the dot product of the vectors. Then, the element characteristic time (Δ * ) is defined as stated in Equation (8a-c).
In Equation (8a,b), the operators ‖ ‖ and | | stand for the Euclidean size of a vector and the absolute value of a scalar quantity, respectively. Moreover, the operator shows the minimum amount. Similarly, we evaluate the characteristic time in a 2D element.
Mathematically, there are no limitations for the time step, when an implicit solution is used ( 1). Moreover, the limitation Δ Δ * (i: the number of the active elements) must be considered for an explicit solution ( 0). In our Crank-Nickelson scheme, we have set 1/2 (as selected in [24]), because it can provide more acceptable results based on our numerical experiments. In this case, the limitation for time step would be as follows: From now on, we briefly use Δ * instead of Δ * , and name it the characteristic time of the two-phase flow domain After solving Equation (5a) for each time step, we can easily reconstruct the interface. For this purpose, we look for the active elements or the subdomains, which contain at least one node with 1 2 ⁄ and another node with 1 2 ⁄ . If 1 2 ⁄ and 1 2 ⁄ are the volume fractions of nodes 1 and 2 on an element edge, we can find the interface intersection point by a simple linear interpolation, as stated in Equation (10).
In Equation (10), ⃗ , ⃗ , and ⃗ are the spatial coordinates of the intersection point, the first node, and the second node, respectively.

Geometrical Re-Initialization
The solution of Equation (5a), using any stabilized finite element, finite volume, or finite difference, introduces an artificial diffusion to the problem. Therefore, it is expectable that the level set function (Φ) deviates from the signed distance function. Although employing a narrow-band technique, such as the local least-square finite In Equation (7a,b), the operator "." demonstrates the dot product of the vectors. Then, the element characteristic time (∆t * ) is defined as stated in Equation (8a-c).
In Equation (8a,b), the operators and | | stand for the Euclidean size of a vector and the absolute value of a scalar quantity, respectively. Moreover, the operator min shows the minimum amount. Similarly, we evaluate the characteristic time in a 2D element.
Mathematically, there are no limitations for the time step, when an implicit solution is used (α = 1). Moreover, the limitation ∆t ≤ min ∆t * i (i: the number of the active elements) must be considered for an explicit solution (α = 0). In our Crank-Nickelson scheme, we have set α = 1/2 (as selected in [24]), because it can provide more acceptable results based on our numerical experiments. In this case, the limitation for time step would be as follows: From now on, we briefly use ∆t * instead of min ∆t * i , and name it the characteristic time of the two-phase flow domain After solving Equation (5a) for each time step, we can easily reconstruct the interface. For this purpose, we look for the active elements or the subdomains, which contain at least one node with f > 1/2 and another node with f < 1/2. If f 1 > 1/2 and f 2 < 1/2 are the volume fractions of nodes 1 and 2 on an element edge, we can find the interface intersection point by a simple linear interpolation, as stated in Equation (10).
In Equation (10), → x int , → x 1 , and → x 2 are the spatial coordinates of the intersection point, the first node, and the second node, respectively.

Geometrical Re-Initialization
The solution of Equation (5a), using any stabilized finite element, finite volume, or finite difference, introduces an artificial diffusion to the problem. Therefore, it is expectable that the level set function (Φ) deviates from the signed distance function. Although employing a narrow-band technique, such as the local least-square finite element method (LLSFEM) [7], prevents the excessive interface diffusion, it does not remove the demand for the re-initialization step.
Different types of re-distancing or re-initialization are discussed in the introduction section. However, in this research, we concentrate on the geometrical re-initialization.

The Motivations and Preliminaries
The geometrical re-initialization comprises the methods which calculate the nodal distance function (DF) using the spatial coordinates of the nodes and the interface. The sign of the DF is determined by the sign of the nodal level set values.
By introducing a finite number of imaginary particles on the interface and searching for the shortest distance between each node and the particles, one can simply evaluate the DF (|d|). The particle-based method is a time-consuming task. However, it could be more efficient by limiting the search algorithm to the narrow band's nodes (e.g., Figure 1).
Taking the particle-based scheme to re-initialize the level set values will lead to a severe mass loss/gain, which is attributed to the interface displacement after re-initialization. To illustrate the movement of the interface, we provide an example. The problem consists of a circle with an arbitrary radius (e.g., R circle = 0.15) in a square (1 × 1) flow domain. The spatial position of the circle center is  element method (LLSFEM) [7], prevents the excessive interface diffusion, it does not remove the demand for the re-initialization step. Different types of re-distancing or re-initialization are discussed in the introduction section. However, in this research, we concentrate on the geometrical re-initialization.

The Motivations and Preliminaries
The geometrical re-initialization comprises the methods which calculate the nodal distance function (DF) using the spatial coordinates of the nodes and the interface. The sign of the DF is determined by the sign of the nodal level set values.
By introducing a finite number of imaginary particles on the interface and searching for the shortest distance between each node and the particles, one can simply evaluate the DF (| |). The particle-based method is a time-consuming task. However, it could be more efficient by limiting the search algorithm to the narrow band's nodes (e.g., Figure 1).
Taking the particle-based scheme to re-initialize the level set values will lead to a severe mass loss/gain, which is attributed to the interface displacement after reinitialization. To illustrate the movement of the interface, we provide an example. The problem consists of a circle with an arbitrary radius (e.g., 0.15 in a square (1 1) flow domain. The spatial position of the circle center is ⃗ 0,0.25 as shown in Figure 3a, and a 2D solid-rotational flow around ⃗ 0,0 is considered ( ⃗ , ). The circle rotates two revolutions around ⃗ 0,0 in a structured 50 50 mesh.  Figure 3a shows that the reconstructed circle after two revolutions (solid line) is considerably smaller than the original one (dashed line), which is due to the mass loss in the circle. The dashed line in Figure 3b quantifies the extent of the mass loss during the revolutions. The solid line in Figure 3b shows the mean deviation of the active element curvature ( ) from the original one (1 ⁄ ), which is mainly attributed to the mass loss and the change of the circle radius. The curvature ( ) and normal vector ( ⃗) in each active element are easily calculated using Equation (11a,b). In these equations, the symbols and are the gradient and divergence operators, respectively.
. ⃗ (11b)  Figure 3a shows that the reconstructed circle after two revolutions (solid line) is considerably smaller than the original one (dashed line), which is due to the mass loss in the circle. The dashed line in Figure 3b quantifies the extent of the mass loss during the revolutions. The solid line in Figure 3b shows the mean deviation of the active element curvature (κ) from the original one (1/R circle ), which is mainly attributed to the mass loss and the change of the circle radius. The curvature (κ) and normal vector ( → n ) in each active element are easily calculated using Equation (11a,b). In these equations, the symbols ∇ and ∇ are the gradient and divergence operators, respectively.
A direct re-initialization method is another approach which is proposed by Cho et al. [26] and Ngo and Choi [27]. In their scheme (the direct method), the DF is the shortest perpendicular distance up to the interface segments. If the perpendiculars' foot is not inside the interface segments, the DF would be the minimum distance of each node up to the interface/element intersections. The direct method estimates the distance function equal to the particle-based scheme. Nevertheless, it requires a shorter computing time because it limits the number of searches to find the minimum distance. This method excludes the nodes of active elements from the re-initialization (the nodes in the dark-grey area of Figure 1) to avoid the displacement of the interface. Therefore, the direct approach shows good mass-conservative properties, although it does not take any measure for mass correction. To investigate this problem, in Figure 4a, we can compare the reconstructed circle (solid line) after five revolutions with the original circle (dashed line). The extent of the mass loss and the mean curvature deviation, as shown respectively by the dashed and solid lines in Figure 4b, are smaller than their corresponding values in Figure 3b. A direct re-initialization method is another approach which is proposed by Cho et al. [26] and Ngo and Choi [27]. In their scheme (the direct method), the DF is the shortest perpendicular distance up to the interface segments. If the perpendiculars' foot is not inside the interface segments, the DF would be the minimum distance of each node up to the interface/element intersections. The direct method estimates the distance function equal to the particle-based scheme. Nevertheless, it requires a shorter computing time because it limits the number of searches to find the minimum distance. This method excludes the nodes of active elements from the re-initialization (the nodes in the dark-grey area of Figure 1) to avoid the displacement of the interface. Therefore, the direct approach shows good mass-conservative properties, although it does not take any measure for mass correction. To investigate this problem, in Figure 4a, we can compare the reconstructed circle (solid line) after five revolutions with the original circle (dashed line). The extent of the mass loss and the mean curvature deviation, as shown respectively by the dashed and solid lines in Figure 4b, are smaller than their corresponding values in Figure 3b.  The solid line in Figure 4b shows an increasing trend for the deviation of the curvature values, which are calculated using Equation (11a,b). These equations employ the nodal volume fractions (or the LS values) on each active element and its neighbor elements to evaluate the value of [7] on the nodes. When the nodes of active grids do not participate in the re-distancing procedure, inaccurate estimation of and are expected, due to unsmooth values of .
Mut et al. [24] introduced a two-step method for estimation of the DF, which comprises initialization and evaluation stages. The initialization step calculates each node's distance function. As described for the particle-based scheme, this step consists of a considerable mass loss/gain due to the interface displacement. Therefore, during the evaluation step, the distance function values are corrected. Mut et al. developed the method for triangular and tetrahedral grids, and we have adapted it to rectangular and hexagonal meshes.
Any interface displacement after the re-initialization step leads to an unwanted variation of the area (Δ ) occupied by the primary phase in each active element. Equation (12a-c) evaluate the mass change in an active element (e.g., element number ). In these equations, the superscripts 0 and 1 stand for the values before and after re-initialization, respectively. Moreover, the operator ⌊ ⌋ calculates the integer part of a real number and Ω depicts the element domain. The symbol is a function that calculates the area (or volume) of the primary phase in an element with nodal volume factions of . The solid line in Figure 4b shows an increasing trend for the deviation of the curvature values, which are calculated using Equation (11a,b). These equations employ the nodal volume fractions (or the LS values) on each active element and its neighbor elements to evaluate the value of ∇ f [7] on the nodes. When the nodes of active grids do not participate in the re-distancing procedure, inaccurate estimation of ∇ f and κ are expected, due to unsmooth values of f .
Mut et al. [24] introduced a two-step method for estimation of the DF, which comprises initialization and evaluation stages. The initialization step calculates each node's distance function. As described for the particle-based scheme, this step consists of a considerable mass loss/gain due to the interface displacement. Therefore, during the evaluation step, the distance function values are corrected. Mut et al. developed the method for triangular and tetrahedral grids, and we have adapted it to rectangular and hexagonal meshes.
Any interface displacement after the re-initialization step leads to an unwanted variation of the area (∆A) occupied by the primary phase in each active element. Equation (12a-c) evaluate the mass change in an active element (e.g., element number i). In these equations, the superscripts 0 and 1 stand for the values before and after re-initialization, respectively. Moreover, the operator calculates the integer part of a real number and Ω e depicts the element domain. The symbol A is a function that calculates the area (or volume) of the primary phase in an element with nodal volume factions of { f }.
To our knowledge, the values of ∆A i are non-zero just for the active elements. To compensate for the area changes, Mut et al. [24] suggested a mass correction value (δ i ) for each active element, which is the root of the area residual (R), as defined in Equation (13): They solved Equation (13) for each active element by employing the scant method, as summarized in Equation (14): In Equation (14), the superscripts (m) to (m + 2) denote the iteration steps. This equation converges after 2 or 3 iterations (δ . The value of δ i is attributed to the nodes of the i th element, which is zero for any non-active element. By averaging δ i values of active elements that share a node, the nodal mass correction values (δ * ) are evaluated. For example, if node j is attached to two active elements i 1 and i 2 , δ * j would be as follows: Then, the corrected volume fraction values ( f * ) are evaluated using Equation (16): In Equation (16), the β parameter is a scalar quantity, adjusted using a mass conservation criterion. Mut et al. [24] calculated β by considering Equation (17), which keeps the main phase area (or volume) unchanged, all over the elements (N), before and after re-initialization. Figure 5 presents the results of the solid circle revolution, employing Equations (12)- (17). Comparing the solid line (the circle after five revolutions) and the original circle (dashed line) in Figure 5a reveals the circle deformation in a solid rotational flow field. The mass loss shown in Figure 5b by the dashed line does not suggest any improvement compared to the direct re-initialization method (Figure 4b). We attribute this to the concept of Equation (17), in which the mass correction is not performed according to the content of the primary phase at the onset of the problem.

Improved Conservative Direct Re-Initialization (ICDR)
In this section, we introduce the improved conservative direct re-initialization (ICDR) method. This approach employs the direct technique [26,27], which is faster than the particle-based scheme, to find the nodal distance function (DF). To reduce the mass loss/gain, we take the Mut strategy [24] and modify it using a global mass conservation concept [12]. The ICDR method, adapted to rectangular and hexagonal elements, involves four steps, as follows:

Improved Conservative Direct Re-Initialization (ICDR)
In this section, we introduce the improved conservative direct re-initialization (ICDR) method. This approach employs the direct technique [26,27], which is faster than the particle-based scheme, to find the nodal distance function (DF). To reduce the mass loss/gain, we take the Mut strategy [24] and modify it using a global mass conservation concept [12]. The ICDR method, adapted to rectangular and hexagonal elements, involves four steps, as follows: Step 1: We evaluate the nodal distance functions of the active elements (e.g., just the nodes in the dark gray section of Figure 1) using the direct technique.
Step 2: By employing the nodal values of the volume fraction before and after redistancing, we calculate the mass correction value of the node ( * ), according to Equations (12)- (15).
Step 3: We replace Equation (17) with Equation (18), which guarantees global mass conservation for each phase. Equation (18) forces the content of the primary phase to be initial value (i.e., 0) plus the part, which is lost or gained during the re-initialization step (∑ ).
Step 4: To obtain a smoothed set of values for the SD or nodal volume fraction, we re-distance the nodes in the vicinity of the interface (e.g., the dark and light gray in Figure  1) using the direct method.
By modifying the mass correction criterion in the third step, we anticipate better mass conservation properties for the ICDR method. Moreover, with a final re-initialization at step (4), a smoother set of nodal volume fraction and a better approximation for curvature values is expectable. These two steps distinguish the ICDR method from the existing geometrical re-initialization methods. To investigate the expected improvements, we compare the results for the rotating circle's problem using the ICDR in Figure 6a  Step 1: We evaluate the nodal distance functions of the active elements (e.g., just the nodes in the dark gray section of Figure 1) using the direct technique.
Step 2: By employing the nodal values of the volume fraction before and after redistancing, we calculate the mass correction value of the node (δ * ), according to Equations (12)- (15).
Step 3: We replace Equation (17) with Equation (18), which guarantees global mass conservation for each phase. Equation (18) forces the content of the primary phase to be initial value (i.e., t = 0) plus the part, which is lost or gained during the re-initialization step ( Step 4: To obtain a smoothed set of values for the SD or nodal volume fraction, we re-distance the nodes in the vicinity of the interface (e.g., the dark and light gray in Figure 1) using the direct method.
By modifying the mass correction criterion in the third step, we anticipate better mass conservation properties for the ICDR method. Moreover, with a final re-initialization at step (4), a smoother set of nodal volume fraction and a better approximation for curvature values is expectable. These two steps distinguish the ICDR method from the existing geometrical re-initialization methods. To investigate the expected improvements, we compare the results for the rotating circle's problem using the ICDR in Figure 6a,b with the corresponding results in Figures 4 and 5.
The circle deformation after five revolutions in Figure 6a seems to be less in comparison to Figures 4a and 5a. Moreover, the dashed line in Figure 6b demonstrates the superior conservative properties of the ICDR method in comparison to two other methods, as shown in Figures 4b and 5b. Similarly, the deviation of the mean curvature values from its exact value (κ = 1/R circle =4) is less for the ICDR method (the solid lines in Figure 4b, Figure 5b, and Figure 6b).
By modifying the mass correction criterion in the third step, we anticipate better mass conservation properties for the ICDR method. Moreover, with a final re-initialization at step (4), a smoother set of nodal volume fraction and a better approximation for curvature values is expectable. These two steps distinguish the ICDR method from the existing geometrical re-initialization methods. To investigate the expected improvements, we compare the results for the rotating circle's problem using the ICDR in Figure 6a

Results and Discussion
In this section, we examine the viability of the ICDR method. For this purpose, we use two benchmark problems to study the interface reconstruction and mass-conservative properties in two-dimensional form. We show that the new re-initialization algorithm is applicable for structured and unstructured meshes. Finally, we use the ICDR method for three-dimensional drop deformation problems.

Comparison of the Geometrical Re-Initialization Methods in 2D
To compare different re-initialization strategies, we use the notched cylinder problem, which was initially introduced by Zalesak [30]. In this problem, a notched cylinder with pre-defined dimensions revolves in a solid-rotational flow field, as introduced in Section 3.1. We expect the notched cylinder to comply with its original shape after each revolution. However, there exists some deviation when using any advection and re-initialization method. A defined norm value of 1 in Equation (19) is the criteria for quantifying the mentioned deviation. In Equation (19), the interface length is shown by L int .
To have a comparison with the literature, we choose the domain with N 1 × N 2 = 200 × 200 rectangular elements, which is the closest one to the cases with 80,000 and 79,668 triangular elements in [24,27], respectively. Besides the time step ∆t = 0.5, we have considered ∆t = 1.0 to compare the effect of time steps in different methods. Notably, all of the results are generated using rectangular sub-domains or elements. Figure 7a,b show the rotated cylinder after the first and fifth revolutions, respectively. In these figures, the dotted and solid lines show the results for ∆t = 1.0 and ∆t = 0.5, correspondingly. As depicted in Figure 7a, the direct method after the first revolution shows some unevenness at the lower part of the notched cylinder. The jagged area, which is similarly reported in Figure 8 of [27], expands during the second to fifth revolutions, as shown in Figure 7b. We attribute this to the nodes of the active elements, which are not re-initialized. We do not expect such deviations for two other methods, as shown in Figure 7a,b, since all of the nodes in the narrow-band are re-initialized.  Besides the visual assessment of the different techniques in Figure 7a,b, the extent of mass loss/gain, and the norm values, are presented in Figure 8a,b to provide a quantitative evaluation. Figure 8a shows that the Mut method has more mass loss than the direct approach. Moreover, the ICDR method has nearly no mass loss/gain due to the global mass correction stated in Equation (18). Figure 8b reveals that the norm values are roughly similar for the ICDR and the Mut methods, and the direct approach has the lowest values. However, this does not mean that the direct method provides better results when considering the reconstructed shapes in Figure 7a,b. Notably, the value of the norm for the Mut approach using rectangular mesh is 1.00 10 , which is very close to the calculated value using triangular meshes [24] ( 1.2 10 ).    Besides the visual assessment of the different techniques in Figure 7a,b, the extent of mass loss/gain, and the norm values, are presented in Figure 8a,b to provide a quantitative evaluation. Figure 8a shows that the Mut method has more mass loss than the direct approach. Moreover, the ICDR method has nearly no mass loss/gain due to the global mass correction stated in Equation (18). Figure 8b reveals that the norm values are roughly similar for the ICDR and the Mut methods, and the direct approach has the lowest values. However, this does not mean that the direct method provides better results when considering the reconstructed shapes in Figure 7a,b. Notably, the value of the norm for the Mut approach using rectangular mesh is 1.00 10 , which is very close to the calculated value using triangular meshes [24] ( 1.2 10 ). Besides the visual assessment of the different techniques in Figure 7a,b, the extent of mass loss/gain, and the norm values, are presented in Figure 8a,b to provide a quantitative evaluation. Figure 8a shows that the Mut method has more mass loss than the direct approach. Moreover, the ICDR method has nearly no mass loss/gain due to the global mass correction stated in Equation (18). Figure 8b reveals that the norm values are roughly similar for the ICDR and the Mut methods, and the direct approach has the lowest values. However, this does not mean that the direct method provides better results when considering the reconstructed shapes in Figure 7a,b. Notably, the value of the norm for the Mut approach using rectangular mesh is 1 = 1.00 × 10 −3 , which is very close to the calculated value using triangular meshes [24] ( 1 = 1.2 × 10 −3 ).

The Notched Cylinder Problem
The displacement of a free surface, such as a notched cylinder's revolution, comprises interface advection Equation (2) and interface re-initialization. Decreasing the time-step size (∆t) improves the interface advection's accuracy; also, it increases the demand for the re-initialization step. The mass loss/gain during re-initialization and its correction can introduce some new sources of error. Therefore, small time-step sizes can decrease the accuracy due to the extra re-initialization steps.
To describe this problem in more detail, we examine different time steps in the allowable range (∆t ≤ 2∆t * ), based on Equation (8) for the notched cylinder problem. Moreover, we consider various structured mesh sizes N 1 × N 2 , as depicted in Figure 1, to investigate the mesh dependency of the regenerated interface. The results are summarized in Table 1.  Table 1 quantitatively show that the increase or decrease of the time increment can destroy the results of an interface advection problem. However, there exists an optimum choice for the time increment (∆t = 1.00∆t * ), which is highlighted by gray color in the tables. Notably, the solution is less sensitive to the time-step size when a finer mesh is employed. For example, for the cases with N 1 × N 2 = 300 × 300, 400 × 400, and 500 × 500, time-step sizes in the range of 0.5∆t * to 1.75∆t * provide acceptable results.
In Figure 9, we visually investigate the effect of the time-step size and mesh spacing by comparing the notched cylinder after one revolution. The thick lines with light grey and dark grey colors are for the cases with N 1 × N 2 = 50 × 50 and 100 × 100, respectively. The regenerated interfaces in this figure reveal that the reconstructed surfaces are nearly similar when N 1 × N 2 ≥ 200 × 200. That is why it is difficult to distinguish between dashed, dash-dotted, dotted, and solid lines, which are for 200 × 200, 300 × 300, 400 × 400, and 500 × 500 meshing systems, respectively. The mesh independence of the results is more considerable for the case with ∆t = 1.00∆t * . Notably, ∆t = 1.00∆t * is also the best choice for time-step size according to Table 1. and dark grey colors are for the cases with 50 50 and 100 100 , respectively. The regenerated interfaces in this figure reveal that the reconstructed surfaces are nearly similar when 200 200 . That is why it is difficult to distinguish between dashed, dash-dotted, dotted, and solid lines, which are for 200 200, 300 300, 400 400, and 500 500 meshing systems, respectively. The mesh independence of the results is more considerable for the case with Δ 1.00Δ * . Notably, Δ 1.00Δ * is also the best choice for time-step size according to Table 1. Figure 9. The reconstructed notched cylinder after one revolution for different time increments (∆t = 0.25∆t * , 1.00t * and 1.75∆t * ) and mesh spacing: Now we take the best choice for the time step (∆t = 1.00∆t * ), and visually investigate the mesh dependency of the regenerated notched cylinder after the first, third, and fifth revolution, in Figure 10a-f. The mesh dependency analysis in these figures shows that the interface prediction is improved by mesh refinement. Moreover, comparing the regenerated notched cylinders after the first, third and fifth revolutions shows an increasing deviation from the original shape over time (revolution), which is less significant for finner meshes.  To quantify the observations of Figure 10, we show the defined norm value ( 1 ) during the revolutions for different mesh sizes in Figure 11a. These figures reveal that the mesh spacings of 50 × 50 and 100 × 100 do not provide a good prediction for the interface displacement. Moreover, for the range of 200 < N 1 the regenerated notched cylinders are nearly independent of the mesh size. To quantify the observations of Figure 10, we show the defined norm value ( ) during the revolutions for different mesh sizes in Figure 11a. These figures reveal that the mesh spacings of 50 50 and 100 100 do not provide a good prediction for the interface displacement. Moreover, for the range of 200 the regenerated notched cylinders are nearly independent of the mesh size.
Although the main focus of this research is the geometrical re-initialization, we have provided some data based on the PDE re-initialization in Figure 11a. By comparing the norm values after the first revolution and PDE-based results, we can infer that the geometrical re-initialization does not always provide better results. However, the ICDR method leads to an acceptable solution (as discussed in the following sections) with less mathematical effort. Figure 11b also shows the convergency order values for different mesh sizes and revolutions.  Although the main focus of this research is the geometrical re-initialization, we have provided some data based on the PDE re-initialization in Figure 11a. By comparing the norm values after the first revolution and PDE-based results, we can infer that the geometrical re-initialization does not always provide better results. However, the ICDR method leads to an acceptable solution (as discussed in the following sections) with less mathematical effort. Figure 11b also shows the convergency order values for different mesh sizes and revolutions.

Swirling Flow Vortex
The swirling flow vortex is another benchmark problem used to evaluate the advection method and re-initialization scheme. This problem comprises a circle in a square flow domain, as introduced in Section 3.1. Similar to the notched cylinder problem, the flow field is independent of the interface shape. Equation (20a,b) define the velocity vectors' components as a function of time and spatial coordinates. In this flow field, the circle stretches like a snake up to t = T/2, and after that time, it returns to its original position. In this research, we assume T = 8.
Figure 12a-f provide a graphical evaluation for the interface calculation results at t = 4 and t = T for structured meshes with different spacing. The dashed line shows the interface's exact shape at the end of the stretching (t = 4) and coming back (t = 8) stages. These figures also contain the values of the defined norm for a quantitative comparison. The mesh dependency analysis in these figures shows that by decreasing the mesh sizes, better results for the interface prediction are achieved. However, we must consider the capacity of our computing system and calculation time before choosing a suitable grid spacing.

Swirling Flow Vortex
The swirling flow vortex is another benchmark problem used to evaluate the advection method and re-initialization scheme. This problem comprises a circle in a square flow domain, as introduced in Section 3.1. Similar to the notched cylinder problem, the flow field is independent of the interface shape. Equation (20a,b) define the velocity vectors' components as a function of time and spatial coordinates. In this flow field, the circle stretches like a snake up to 2 ⁄ , and after that time, it returns to its original position. In this research, we assume 8.
Figure 12a-f provide a graphical evaluation for the interface calculation results at 4 and for structured meshes with different spacing. The dashed line shows the interface's exact shape at the end of the stretching ( 4) and coming back ( 8) stages.
These figures also contain the values of the defined norm for a quantitative comparison. The mesh dependency analysis in these figures shows that by decreasing the mesh sizes, better results for the interface prediction are achieved. However, we must consider the capacity of our computing system and calculation time before choosing a suitable grid spacing.

Unstructured Meshes
Notwithstanding that the description of the ICDR method in the third section does not rely on a structured mesh, we provide some examples to evaluate the capability of the introduced method in unstructured meshes. In this section, we take the assumption of Δ Δ * . Note that the values of and are the integer part of the second root of the total element numbers.
We start with the notched cylinder problem. To make the results visually comparable with corresponding regenerated interfaces in Figure 10, we have performed a mesh dependency analysis, as shown in Figure 13a-e. Similar to Figure 10, these figures show the regenerated notched cylinder after the first, third, and fifth revolutions.

Unstructured Meshes
Notwithstanding that the description of the ICDR method in the Section 3 does not rely on a structured mesh, we provide some examples to evaluate the capability of the introduced method in unstructured meshes. In this section, we take the assumption of ∆t = ∆t * . Note that the values of N 1 and N 2 are the integer part of the second root of the total element numbers.
We start with the notched cylinder problem. To make the results visually comparable with corresponding regenerated interfaces in Figure 10, we have performed a mesh dependency analysis, as shown in Figure 13a-e. Similar to Figure 10, these figures show the regenerated notched cylinder after the first, third, and fifth revolutions.

Unstructured Meshes
Notwithstanding that the description of the ICDR method in the third section does not rely on a structured mesh, we provide some examples to evaluate the capability of the introduced method in unstructured meshes. In this section, we take the assumption of Δ Δ * . Note that the values of and are the integer part of the second root of the total element numbers.
We start with the notched cylinder problem. To make the results visually comparable with corresponding regenerated interfaces in Figure 10, we have performed a mesh dependency analysis, as shown in Figure 13a-e. Similar to Figure 10, these figures show the regenerated notched cylinder after the first, third, and fifth revolutions. As it was reported for the structured meshes, we observe in Figure 13a-e that the reproduced interfaces comply better with the original shape by increasing the mesh numbers. To quantify the results of these figures, the values of the defined norm during the first to fifth revolutions are reported in Figure 14a,b also shows the order of the convergency for the norm.  As it was reported for the structured meshes, we observe in Figure 13a-e that the reproduced interfaces comply better with the original shape by increasing the mesh numbers. To quantify the results of these figures, the values of the defined norm during the first to fifth revolutions are reported in Figure 14a,b also shows the order of the convergency for the 1 norm. As it was reported for the structured meshes, we observe in Figure 13a-e that the reproduced interfaces comply better with the original shape by increasing the mesh numbers. To quantify the results of these figures, the values of the defined norm during the first to fifth revolutions are reported in Figure 14a,b also shows the order of the convergency for the norm.   Figure 14a shows that for the range of 204 < N 1 the regenerated interfaces are nearly independent of the mesh size, similar to the cases with structure mesh in Figure 11. By comparing the data of Figures 11a and 14a, we can surprisingly find that the values of the 1 norm of the unstructured meshes is less than the equivalent norms for the structured meshes. Figure 15a-e illustrate the results of our mesh dependency analysis for the swirling flow vortex using unstructured meshing systems. As expected, better results go for the cases with finer mesh. Paralleling the norm values in Figures 12 and 15 reveals that for the swirling flow vortex problem, a structured mesh can provide better results.
Fluids 2021, 6, x FOR PEER REVIEW 18 of 25 Figure 14a shows that for the range of 204 the regenerated interfaces are nearly independent of the mesh size, similar to the cases with structure mesh in Figure 11. By comparing the data of Figures 11a and 14a, we can surprisingly find that the values of the norm of the unstructured meshes is less than the equivalent norms for the structured meshes. Figure 15a-e illustrate the results of our mesh dependency analysis for the swirling flow vortex using unstructured meshing systems. As expected, better results go for the cases with finer mesh. Paralleling the norm values in Figures 12 and 15 reveals that for the swirling flow vortex problem, a structured mesh can provide better results.

The Extension to Three-Dimensional (3D) Problems
Aside from the benchmark problems, the main target of a re-initialization method is to predict the interface position in a physical phenomenon accurately. Various applications involve two-phase flow problems. However, we focus on the deformation of a drop (the second phase) in a simple shear flow. We have two reasons for this choice. The first reason is that in the drop deformation problem, despite the benchmark problems, the interface's curvature influences the flow field via surface tension forces. Expectedly, an inaccurate interface displacement during re-initialization and mass correction can destroy the flow field calculations and the predicted drop shape. Therefore, a drop deformation problem is a good checkpoint for our re-initialization method. Secondly, there exist equivalent numerical and experimental data in the literature to compare our results with them and validate the ICDR method in 3D. Figure 16 shows a typical two-phase simple shear flow field. In this picture, a drop with the radius of R drop is assumed in the middle of a channel with the dimensions of W × H × L. The upper and lower plates are assumed to move in channel longitudinal direction with the velocities of +U/2 and −U/2; respectively. The extent of the drop deformation or breakup is a function of the viscosity ratio (λ) and Capillary number (Ca), which are defined in Equation (21a,b), where the parameters µ 1 and µ 2 show the viscosity of the matrix (the primary phase)and drop (the second phase), respectively. Moreover, the parameter γ stands for the surface tension coefficient.

The Extension to Three-Dimensional (3D) Problems
Aside from the benchmark problems, the main target of a re-initialization method is to predict the interface position in a physical phenomenon accurately. Various applications involve two-phase flow problems. However, we focus on the deformation of a drop (the second phase) in a simple shear flow. We have two reasons for this choice. The first reason is that in the drop deformation problem, despite the benchmark problems, the interface's curvature influences the flow field via surface tension forces. Expectedly, an inaccurate interface displacement during re-initialization and mass correction can destroy the flow field calculations and the predicted drop shape. Therefore, a drop deformation problem is a good checkpoint for our re-initialization method. Secondly, there exist equivalent numerical and experimental data in the literature to compare our results with them and validate the ICDR method in 3D. Figure 16 shows  By taking the assumptions of the creeping flow and incompressible fluids, the governing equations, including the equation of motion and mass conservation criterion, would be as follows: The Equation (22a,b) are based on Einstein summation notations over a common index. The components of the total stress tensor ( = σ) are denoted by σ ij , which are described in Equation (23a-c).
where τ ij and D ij are components of the viscous stress ( = τ) and deformation rate ( = D) tensors; respectively. P and µ are hydrostatic pressure and Newtonian viscosity, correspondingly. I ij is the unit tensor (I) component (Kronecker delta), such that I ij is 1 if i = j and is 0 if (21a) is the component of the surface tension force per unit volume ( → F σ ) and is easily calculated using Equation (23) [32].
The solution strategy is to employ the Galerkin finite element method (FEM) for Equation (22a,b). Moreover, we consider a periodic boundary condition on the front/back and the left/right sides. The surface tension force ( → F σ ) is considered as a source term, which is non-zero in the active elements. Mostafaiyan et al. have described the solution method's detail (for the two-dimensional case) in [7].
Our setup in this section is to employ a structured mesh in the flow domain with hexagonal elements. We assume 40 grids in the depth of the channel (element size is H/40), for which we have obtained mesh-independent results in our two-dimensional analysis [7]. We assume the time increment to be slightly smaller than the characteristic time (e.g., ∆t = 0.9∆t * ).
In this section, we examine the proposed method by juxtaposing our results with some numerical and experimental data in the literature. For this purpose, the extent of the drop deformation and its breakup are presented for unit and non-unit viscosity ratios.

Comparison of the Geometrical Re-Initialization Methods in 3D
We consider a case with Ca = 0.4 and 2R drop = 0.5H. Figure 17a shows the results of our calculations, which are re-initialized using the ICDR method. In this figure, a binary breakup is observed, which complies with the reported data in [33] (Figure 13). Notably, for this problem the Capillary number of 0.4 is a critical value, below which no break up will happen [33].
Simulation of the problem using the Mut method does not lead to a breakup, as shown in Figure 17b. The simulation time equals in Figure 17a,b. However, in the recent case, the drop reaches its equilibrium before any breakup. We attribute the difference between the two figures to the second phase's mass loss during deformation. The drop losses part of its mass and converts to a smaller one (i.e., a drop with a smaller R drop ) in the simulation course. Thus, according to Equation (21b), the Capillary number drops below the critical value and will not break up, as a result.
We have also employed the direct re-initialization algorithm to estimate the extent of the drop deformation. In this case, the interface is destroyed during the deformation, as shown in Figure 17c, and we cannot continue simulation up to the breakup time in Figure 17a, due to the numerical errors. The unsuccessful simulation's main reason is the inaccurate prediction of the surface tension forces due to the volume fraction's unsmooth values in the interface's vicinity. It is worth reminding that the active elements' nodal values are not re-initialized in the direct algorithm.  Figure 17a shows that our analysis complies with the numerical results in the literature [33]. In this section, we use our two-phase flow strategy to simulate an experimental task of Sibillo et al. [34] Figure 17a shows that our analysis complies with the numerical results in the literature [33]. In this section, we use our two-phase flow strategy to simulate an experimental task of Sibillo et al. [34] with 2R drop = 0.7H and Ca = 0.46. Their result (a quinary breakup) are reproduced by our simulation and compared in Figure 18a,b.

Case Study for Drop Deformation with Non-Unit Viscosity Ratios
Since the flow field contains two different fluids with the viscosity values of µ 1 and µ 2 (µ 1 = µ 2 ), an appropriate viscosity estimation (µ) in Equation (23b) is required. We have estimated the viscosity using the area portion of an element occupied by the primary and second phases, as shown in Equation (25a,b).  To evaluate our re-initialization method for problems with non-unit viscosity ratio, we evaluate the steady shape a drop with 2 0.9 and 0.3 for different viscosity ratios (i.e., 10 and 1 10 ⁄ ) and compare it with similar data in the literature. Figure 19a,b show the shape of the drops at equilibrium, generated using the Boundary-integral method [35]. Figure 19c,d present the corresponding results while using our finite element strategy. For drop deformation with a non-unit viscosity ratio, the breakup is possible. Vananroye et al. [36] have experimentally reported such a breakup for a drop with 2 0.75 , 0.6, and 11.7, as shown in Figure 20b. Figure 20a depicts the reproduced results using our FEM formulation and the ICDR technique. The similarity between our calculations and different experimental and numerical data in the literature ( 1 or 1 ) reflects the ICDR technique's capability in conjunction with our FEM in analyzing two-phase flow problems.
It is worth mentioning that we have performed several successful comparisons for the drop deformation problems (with the unit and non-unit viscosity ratios). However, they are not reported here for conciseness. To evaluate our re-initialization method for problems with non-unit viscosity ratio, we evaluate the steady shape a drop with 2R drop = 0.9H and Ca = 0.3 for different viscosity ratios (i.e., λ = 10 and λ = 1/10) and compare it with similar data in the literature. Figure 19a,b show the shape of the drops at equilibrium, generated using the Boundaryintegral method [35]. Figure 19c,d present the corresponding results while using our finite element strategy.  For drop deformation with a non-unit viscosity ratio, the breakup is possible. Vananroye et al. [36] have experimentally reported such a breakup for a drop with 2R drop = 0.75H, Ca = 0.6, and λ = 11.7, as shown in Figure 20b. Figure 20a depicts the reproduced results using our FEM formulation and the ICDR technique. The similarity between our calculations and different experimental and numerical data in the literature (λ = 1 or λ = 1) reflects the ICDR technique's capability in conjunction with our FEM in analyzing two-phase flow problems.
It is worth mentioning that we have performed several successful comparisons for the drop deformation problems (with the unit and non-unit viscosity ratios). However, they are not reported here for conciseness.

Conclusions
We have introduced an improved mass-conservative direct re-initialization method (ICDR) based on two existing geometrical re-initialization methods (the Mut and the direct method) and a global mass correction concept. This method has excellent mass-conservative properties and provides a smooth set of nodal volume fractions in the interface's vicinity. We have shown that the ICDR method is applicable for 2D and 3D flow problems.
To assess the validity of the ICDR method in 2D, we have employed two benchmark problems. Then the dependency of the results on the element and time-step size has been extensively examined. A visual assessment using graphical demonstrations of the regenerated interfaces, as well as a quantitative evaluation via the definition of a norm value, have been provided.
Aside from the benchmark problems, we have decided to use an actual physical phenomenon to examine the ICDR method's applicability in 3D. For this purpose, a drop deformation in a simple shear flow field has been selected. Comparing our results with various numerical and experimental data in the literature has shown that, using the introduced strategy, one can predict the droplet equilibrium or break up for unit and non-unit viscosity ratios.
The 3D flow fields in this research are limited to creeping flow. However, the ICDR algorithm's evaluation by considering the inertial terms (e.g., bubble rising and dam-break problems) can be a subject of further studies.