Optimization-based Level-Set Re-initialization: A Robust Interface Preserving Approach in Multiphase Problems

In spite of its overall efficiency and robustness for capturing the interface in multiphase fluid dynamics simulations, the well-known shortcoming of the level-set method is associated with the lack of a systematic approach for preserving the regularity of the distance function. This is mainly due to the stretching (or compressing) effect of the strain rate especially in the vicinity of the liquid-gas interface. Level-set re-initialization is an effective treatment for this issue. However, the traditional approach based on the hyperbolic Hamilton-Jacobi equation is a computationally expensive procedure. Crucially, due to the hyperbolic nature of the formulation, the accuracy of the results hinges significantly on the specialized handling of blind spots near the liquid-gas interface intersecting the substrate. The present work proposes a two-step elliptic level-set re-initialization approach that strictly preserves the location of zero level-set via incorporation of an element splitting process. The primary initialization step helps remove any non-smoothness in the to-be regularized level-set function dramatically improving the efficiency of the secondary optimization step. Geometric representation of the boundary conditions is utilized in the initialization step, while the optimization step minimizes the reliance of the results on the treatment of the blind spots. The performance of the proposed method is examined for free and sessile three-dimensional droplets.


Introduction
Since its first introduction in the context of computational fluid dynamics (CFD) [1], the level-set method [2,3] has extensively been used for multiphase flow problems [4,5,6,7].By capturing the phase interface evolution on a fixed computational mesh, the level-set method circumvents the need for complicated deformation and regeneration of the mesh.While the alternative technique, volume of fluid (VOF) [8], is usually preferred in the context of the finite volume method [9,10], the level-set method is the most widely used technique in the context of the finite element simulation of multi-phase flows [11,12].This is mainly due to the ease of calculating the major geometrical properties of the phase interface, i.e., the normal vector and curvature.
The main drawbacks of the level-set method are the lack of the mass conservation property [13] and its proneness to irregularities [11].The former is usually alleviated by using higher-order level-set convection approaches [14,15] and introducing corrective terms to the formulation [16,17].On the other hand, the irregularities are mainly caused by the local expansion and contraction of the distance between the convected surfaces of iso-level-set functions [18].This per se, is an expected result of a non-zero strain rate corresponding to the convection velocity field.In order to address such irregularities, one can either prevent (or minimize) their occurrence during the convection of the level-set [19,20], or, alternatively, frequently apply a levelset (distance) re-initialization process [21].While the former treatment has not yet been widely assessed, the latter option has become the common way of dealing with the irregularities in the convected level-set function.It is worth mentioning that in order to minimize the aforementioned drawbacks, some works have also been dedicated to the incorporation of the (semi-)Lagrangian interface tracking approaches, e.g., the particle level-set method [22].
The preservation of the configuration of the zero level-set surface is of utmost importance for a robust level-set re-initialization method [23], as it defines the phase interface in the context of two-phase flow.Failure to maintain the zero level-set during distance re-initialization not only introduces geomet-rical inaccuracies but also undermines mass conservation [14,24].Therefore, the so-called "anchoring" of the interface [25] is of essential importance in the development of distance re-initialization techniques.
The conventional approach used level-set re-initialization approach relies on the iterative solution of the Eikonal equation [4,26,27], characterized by a hyperbolic mathematical formulation governing the regular propagation of the iso-level-set surfaces [28].While widely used, this method exhibits poor efficiency when the initial state deviates significantly from the signed distance function [29].Moreover, for severely irregular cases, the preservation of the configuration of the zero level-set cannot be guaranteed [30].A potential remedy to the above-mentioned issue would be to rely on the geometrically evaluated distance to the fixed zero level-set.[31].Nevertheless, these approaches that basically reproduce the minimum distance would suffer from blind spots in the vicinity of the cut in the zero level-set, leading to failure in scenarios such as droplet dynamic simulations involving a contact line [32].In an alternative approach, Elias et al. [33] proposed to march over the elements (computational cells) and analytically calculate the unknown nodal value of the level-set function.This process that starts from the vicinity of the interface (zero level-set), however, occasionally disturbs the smoothness of the level-set function.Furthermore, sweeping the whole domain, far from the interface, makes the procedure computationally expensive.
Another class of potentially effective but infrequently utilized level-set re-initialization methods is based on solving an optimization problem [34,35,36].Governed by an elliptic equation, these methods do not encounter blind spots in the domain.By employing such a method, preserving the configuration of the zero level-set can be achieved through the introduction of appropriate penalty terms.Nevertheless, the non-linearity of the governing equation increases the sensitivity of the optimization procedure to the initial state (of the level-set function).Improper treatment of a severely non-regular initial state can result in high computational costs or even lead to the failure of the optimization procedure.
To address the limitations of optimization-based level-set re-initialization approaches, this work proposes a novel strategy, which begins by employing a reconstruction scheme, ensuring a smooth initial state for the subsequent optimization process.In essence, our approach consists of an elliptic initialization in the first step, followed by optimization of the initial guess in the second step.A geometric representation of the level-set gradient featuring a sophisticated blind spot treatment is proposed to handle the Newmann boundary condition required by the primary elliptic re-initialization step.The optimization process, per se, is also enhanced by introducing a modified potential function.Moreover, the proposed method utilizes an element splitting procedure, which facilitates the accurate incorporation of the interfacial penalty terms and consequently, ensures the preservation of the zero level-set configuration during the re-initialization process.
In the following section, the optimization process along with the proposed potential function are elaborated.Then, the elliptic reconstruction scheme and the imposition of the boundary conditions are detailed, focusing on the geometrical approach proposed to deal with the blind spots.Before closing this section, the element splitting technique is also briefly described.The method is verified through a series of static droplet simulations, and finally, its performance is shown for some dynamic cases.

Formulation
Since its first application in the context of the multi-phase methods [4,5], the level-set method has been known as a competent phase-interface capturing approach [3].Its application is based on the convection of smooth field variable ϕ(x) representing the signed distance to the interface (|∇ϕ| = 1).This reads ∂ϕ ∂t where u(x) is the associated velocity field.Here, x is the location vector within the problem domain of Ω ⊂ R d with d denoting the spatial dimensions.The level-set method works well while the regularity of the level-set field is preserved.In most cases, u deviates from a uniform velocity field characterized by a zero strain rate, i.e., ∇u = 0).This explains the expansion and contraction of the distance between iso-ϕ surfaces in the area with positive and negative strain rates, respectively.In the mathematical terms, n • (n • ∇u) > 0 and n • (n • ∇u) < 0, with n being the normal vector to the local iso-ϕ surface, correspond to the expansive and contractive behaviors, respectively.

Distance minimization problem
Seeking for a regularized level-set field, one ideally wants to enforce |∇ϕ| = 1.As first proposed by Li et al. [34], this process can be based on the minimization of the potential function P (|∇ϕ|) with the property of arg min s P (s) = 1. ( Having P defined locally, the total virtual energy of the whole domain is obtained as and the sought for minimization problem becomes

Variational formulation
The problem presented by Eq. ( 4) for energy functional R, is equivalent to finding ϕ that satisfies where δϕ denotes the variation in the level-set function.In order to calculate dR, it is essential to take into account the directional variation of ∇ϕ with respect to δϕ.To this end, the definition of the Gateaux differential can be used as Substituting the definition of the energy functional (3) in Eq. ( 6), one obtains where Consequently, the variational form of the minimization problem becomes with test function ψ = δϕ representing the admissible variations.

Maintaining zero level-set
In order to retain the configuration of the surface defined as Γ(ϕ 0 ) := {x ∈ Ω|ϕ(x) = ϕ 0 (x)} , the total energy functional can be augmented with an external energy function, E, whose minimum occurs at Γ(ϕ 0 ) [34].As pointed out by Basting and Kuzmin [35], a simple but effective choice for E is This external energy can be introduced into the minimization problem, Eq. ( 4), as a penalty term where α is the penalty coefficient.The resulting variational form of the minimization problem reads Alternatively, one can acquire the Lagrange multiplier approach [37] to enforce the constraint on the zero level-set surface.
It is worth mentioning that Eq. ( 11) characteristically presents a nonlinear elliptic problem with the diffusion coefficient of Nevertheless, solving Eq. ( 11) does not require the introduction of any boundary conditions and can readily be linearized and further discretized to obtain the linear system of equations for nodal values of ϕ as the unknowns.Following Basting and Kuzmin [35], the linearization can be based on a fixed-point algorithm, for which the m-th iteration step reads ( 13) Moreover, it should be noted that in most of the applications, the preservation of the configuration of the zero level-set surface is of main importance, i.e. ϕ 0 = 0 in the above equations.Therefore, in the following for more simplicity, Γ represents the zero level-set surface.

Objective potential
The satisfactory performance of the described minimization problem for the re-initialization of the level-set function is centrally determined by the proper definition of the potential function.As elaborated in Eq. ( 2), P (|∇ϕ|) has the essential property of exhibiting a (local) minimum at |∇ϕ| = 0.According to this property, one can characterize the diffusion coefficient (12) to tend to zero as |∇ϕ| → 1. Accordingly, one trivial option for P is which gives It is easily seen that while the problem is well defined for |∇ϕ| ≥ 1 choosing this potential function, d (|∇ϕ|) → −∞ as |∇ϕ| → 0. Considering the computational difficulties associated with a large negative (anti-)diffusion coefficient, it is favorable to modify the definition of P for |∇ϕ| < 1.
To cope with this issue, trigonometric [34] and polynomial [35] doublewell functions are proposed that pose a second local minimum at |∇ϕ| = c, where 0 < c < 1.However, this treatment is not suitable for general levelset applications without any guarantee for c < |∇ϕ|.As a more practical treatment, one can introduce a tailor-made polynomial that keeps the antidiffusivity finite for the whole range of |∇ϕ| < 1. Satisfying the essential conditions of d(1) = 0, and keeping d < 0 for |∇ϕ| < 1, a viable option is to define This coefficient leads to the definition of the potential function proposed by Adams et al.
[37] as In this work, a further relaxed definition for the anti-diffusive coefficient in the range of |∇ϕ| < 1 is proposed reading The rate of the variation of d by changing |∇ϕ| is the same on both sides of the inflection point at |∇ϕ| = 1 (see Fig. 1).This is obtained by reflecting d(s) for 1 ≤ s ≤ 2 around (1, 0).Numerical experiments showed that this would slightly improve the rate of convergence of the problem presented in Eq. ( 13).

Elliptic reconstruction scheme
A major drawback of the optimization-based level-set re-initialization approach is its low rate of convergence.This issue is exacerbated in cases with a significant local reduction of |∇ϕ| with respect to the optimal value of unity.Therefore, the performance of the method can be dramatically enhanced if ϕ was initially recalculated in a way that |∇ϕ| remains reasonably close to unity, e.g., 0.5 < |∇ϕ| < 2, within the computational domain.
In this work, the proposed scheme for the reconstruction of ϕ is based on the solution of with a prescribed source term q.This equation is subject to Neumann conditions defined on the outer boundaries of the domain, ∂Ω.The corresponding weak form of Eq. ( 19) that includes the penalization term for preserving zero level-set reads ( 20) Considering Eq. 19 and taking into account that the ideal solution would feature with |∇ϕ| = 1, one can obtain the admissible definition of q (x) = ∇ • (∇ϕ/|∇ϕ|).In other words, source term q (x) in Eq. 19 is equal to the local curvature of the iso-ϕ surface passing through x.

Boundary condition
As presented in Eq. ( 20), the successful application of the proposed levelset reconstruction scheme is subject to the natural imposition of the (estimated) gradient of the level-set function on the boundaries of the computational domain.Nonetheless, n • ∇ϕ should be consistently evaluated not to disturb the configuration of the contact line.In this work, this boundary term is calculated using the geometrical definition of the regularized levelset function.The technique discussed below can also be readily utilized to improve any geometrical level-set re-initialization technique (e.g., see [31]).
Knowing that for an idea level-set function, the normal vector to the phase interface reads n Γ = ∇ϕ, one can project n Γ onto the outer boundaries to estimate the corresponding boundary condition.As illustrated in Fig. 2, at each boundary node, the projected normal vector originates from the nearest point on the phase interface (the preserved zero level-set).Even though this simple procedure works fine in the absence of contact line, once the zero levelset contour is cut by an outer boundary, a specific treatment is necessary for estimating the boundary condition in the blind spots shown in Fig. 3.Such a condition frequently occurs in different applications, for example, in the case of a sessile droplet on a solid substrate.
In order to treat the blind spot, the associated boundary nodes are divided into two different groups; those lying on the cut plane (see Fig. 4) and those lying on the boundaries perpendicular to the cut plane (see Fig. 5).The treatment here is based on hypothetically extending the zero level-set contour assuming a constant curvature behind the cut plane.The procedure is illustrated in Figs. 4 and 5 for both groups of the boundary nodes in the blind spot.In this way, one can approximate the normal gradient of the   level-set function within the blind spot as and for the first and second groups of the boundary nodes, respectively.It is worth emphasizing that while solving Eq. ( 20) with the local curvature as the source term results in a relatively regularized ϕ, the absence of the aforementioned minimization problem (Eq.( 10)) offers no assurance of achieving |∇ϕ| ≈ 1.The approximations made for estimating the boundary condition term can also potentially increase the deviation from the ideal level-set function.On the other hand, the minimization procedure governed by Eq. ( 10) does not require the introduction of any boundary conditions, and therefore, can be seen as a robust corrector for the proposed level-set re-construction scheme with a low sensitivity to the possible errors in the estimation of the boundary condition term.

Element splitting and preserving contact line
As elaborated earlier, the imposition of the zero level-set constraint is one of the most important requirements for the success of the level-set reinitialization process.In this regard, for the described class of elliptic approaches, one needs to accurately calculate the corresponding penalty terms.To this end, not having ready access to the zero level-set surface on the fixed computational mesh, one can acquire a mechanism to project the terms into the nodal contributions [21].However, the more versatile and straightforward approach would be based on the element-splitting process described in this section.
The splitting of the cut elements starts with the introduction of the corners of the phase interface (ϕ = 0) plane on the edges of the element.For the linear tetrahedral elements used in this work, depending on the configuration, this leads to either a quadrilateral or a triangular cut plane representing the elemental portion of the phase interface.In the case of the quadrilateral shape, as illustrated in Fig. 6, the cut plane is further split into two triangles.Gauss quadrature is then used to integrate the contribution of the penalty term ensuring the preservation of the location of the phase interface.As previously mentioned, in droplet dynamics simulations, accurate capturing of the contact angle plays an essential role.This requirement urges the preservation of the location of the contact line while the zero level-set contour is fixed.With this aim, one can simply extend the penalty terms appearing in Eqs. ( 13) and (20) to read ( 23) where ∂Γ denotes the contact line, which is represented by the cut in the zero level-set contour.Coefficients α Γ and α ∂Γ are chosen large enough to assure the fixation of the phase interface (including the contact line).In this work, we have presented our formulation according to the conventional penalty method to maintain generality and enhance readability.However, for the sake of improved numerical stability and efficiency, our numerical implementation is based on Nitsche's scheme [38].
In summary, our proposed re-initialization method is a two-step process: in the first step, the level-set function is fully reconstructed using Eq. 24, ensuring a smooth initial state.Subsequently, re-initialization is achieved through an iterative optimization procedure based on an elliptic equation Eq. 23, by incorporating a feasible objective potential (Eq.18), and considering element splitting approach of Section 2.3 within.Notably, special treatment, detailed in Section 2.2.1, is applied when handling blind spots in boundary condition calculations.This two-step approach ensures accurate preservation of the zero level-set configuration and mitigates the issues commonly encountered in traditional re-initialization methods.

Results
The method has been implemented within Kratos Multiphysics framework, an in-house Open Source C++ object-oriented Finite Element platform [39,40].To simulate the fluid dynamics, we employ the enriched Finite Element model introduced by the authors in previous works [11,41].Furthermore, the dynamics of the contact line is modelled using a combination of the linear molecular kinetic theory and hydrodynamic theory, as previously detailed in our earlier publication [6].
Unless specified otherwise, the properties of the liquid and gas phases correspond to those of water and air, respectively.These properties are set as follows: water density, ρ l = 1000 kg/m 3 ; water dynamic viscosity, µ l = 0.001 Pa s; air density, ρ g = 1 kg/m 3 ; air dynamic viscosity, µ g = 0.00001 Pa s; and the surface tension of the water-air interface, γ = 0.072 N/m.Gravity is consistently defined as g = 9.81 m/s 2 in all test cases, except for the Oscillating Droplet scenario (Sec.3.2.1),where it is neglected.The simulations are carried out within a cubic domain measuring 1 mm × 1 mm × 1 mm, with the exception of the Droplet in Channel (Sec.3.2.2),where the domain dimensions are different.
Three key capabilities and goals are aimed to be achieved by our approach: 1) maintaining the magnitude of the level-set gradient at unity, i.e., |∇ϕ i |= 1, 2) accurately preserving the liquid-gas interface, and 3) ensuring the preservation of the contact line.To test the method's performance in achieving these goals, appropriate error measures are defined for each objective.
For the first goal, ensuring |∇ϕ i |= 1, the L 1 -norm of the error in the magnitude of the level-set gradient, denoted as is utilized as a quantitative measure, where N Ω ′ is the total number of nodes in the subdomain Ω ′ ⊂ Ω.This error metric allows the assessment of the method's capability in maintaining the correct magnitude of the gradient throughout the optimization process.It's important to note that our evaluation is not limited to the entire domain or full-band, i.e., Ω ′ ≡ Ω, but extends to a narrow-band defined as Ω ′ := {x ∈ Ω|−3h < ϕ(x) < +3h}.This extension allows us to evaluate the method's performance in the close vicinity of the liquid-gas interface, where elements are intersected by the interface, necessitating element splitting.
To evaluate the second desired capability, namely, accurate preservation of the interface Γ, the Chamfer Distance (CD) between the reinitialized interface, ϕ reinit 0 , and the exact interface, ϕ exact 0 , is employed.The Chamfer Distance is calculated as where N Γ is the total number of nodes on the interface Γ, and x i represents the position of the i-th node on Γ.Here, Γ exact represents the set of nodes on the exact interface.
Similarly, the Hausdorff Distance (HD) is used as a metric to assess the method's ability to preserve the contact line ∂Γ, and can be computed as follows: HD(∂Γ) = max sup where ∂Γ reinit and ∂Γ exact represent the sets of nodes on the reinitialized and exact contact line ∂Γ respectively.

Verification
In this subsection, we verify the performance of the proposed level-set reinitialization method using two benchmark test cases: a) Spherical Droplet and b) Ellipsoid Droplet.Despite their apparent simplicity, these cases serve as fundamental building blocks to evaluate the method's generalizability to more complex scenarios.Additionally, the presence of analytical solutions enables direct comparison with numerical results, validating the accuracy and reliability of our proposed approach.

Spherical droplet
The first benchmark test case involves simulating a sessile spherical droplet within a fixed computational domain.This case is highly relevant as it allows us to explore the method's ability to handle sessile droplets and preserve the contact line, representing the three-phase (liquid-gas-solid) situation.The presence of the contact line is crucial, as it highlights the method's efficacy in dealing with complex interfacial dynamics during re-initialization while accurately preserving the zero level-set configuration.
Side views of three spherical droplets on a substrate with different contact angles (Obtuse angle, Right angle, and Acute angle) and fixed radius R = 250µm are presented in Fig. 7.The contours of iso-level-set are displayed for each case, providing a clear visualization of their regularity.
The evolution of the (average) L 1 -norm of the error in the gradient of the level-set function during the iterative optimization procedure is illustrated in Fig. 8.It reveals that the error rapidly decreases to its minimum within the initial iteration steps of the iterative optimization procedure.Notably, the plot displays consistent results for different mesh sizes, ranging from coarser to finer meshes.This robust convergence behavior demonstrates the effectiveness of the proposed method in achieving accurate and efficient optimization of the level-set function, regardless of the mesh resolution.
The left column of Fig. 9 reveals a comprehensive depiction of the convergence behavior of the proposed re-initialization method.This illustration is achieved by systematically reducing the mesh size (h) and computing the (average) L 1 -norm of the error in the gradient of the level-set function for the full-band.The plot demonstrates a convergence pattern that closely aligns with second-order convergence.
The right column of Fig. 9 provides an analogous analysis, yet the focus shifts to the narrow-band of −3h < ϕ < 3h.The convergence trends  observed here exhibit a remarkable alignment with the third-order convergence line.This observation indicates that the proposed optimization-based re-initialization approach is highly beneficial in terms of accuracy.Together, these convergence analyses, depicted in Fig. 9, collectively underscore the accuracy and robustness of our proposed method in maintaining the magnitude of the level-set gradient (|∇ϕ|= 1) even in the presence of intricate geometries and dynamic simulations.In Fig. 10 (left column), the Chamfer Distance (Eq.26) plots for the three spherical droplets with different contact angles are presented.The observed third-order convergence behavior indicates the method's accuracy in preserving the interface, highlighting the effectiveness of the element-splitting scheme and optimization approach in our method.
Similarly, the right column of Fig. 10 displays a first-order convergence trend in Hausdorff Distance (Eq.27), showcasing the method's capability in accurately preserving the contact line even in relatively coarser meshes.
These figures collectively underscore the effectiveness and accuracy of the proposed level-set re-initialization method in preserving both the interface and contact line for spherical interfaces.

Ellipsoidal droplet
Building upon the spherical droplet case, we introduce a little added complexity by simulating ellipsoidal droplets with three semi-axes lengths of a = 350µm, b = 250µm, and c = 200µm along the x, y, and z axes, respectively.This allows us to demonstrate the generalizability of the proposed method to handle non-axisymmetric shapes.In Fig. 11, side views of iso-level-set contours for each side view, demonstrating the regularity of the reinitialized level-set.Fig. 12 demonstrates second-order convergence for the L 1 -norm of |∇ϕ| in both the full-band and narrow-band cases of the ellipsoidal droplet.Furthermore, as illustrated in Fig. 13, the proposed method showcases secondorder convergence for accurately preserving the ellipsoidal liquid-gas interface (ϕ = 0 iso-surface), while maintaining first-order convergence in preserving the contact line.

Dynamic Test Cases 3.2.1. Oscillating Free Droplet
This dynamic scenario explores the evolution of the droplet's interface, driven solely by surface tension forces.This test case introduces additional complexity to the simulation, as the surface tension forces cause the droplet to deform and change shape over time (Fig. 14).It allows us to examine the method's capability to preserve the zero level-set surface accurately while the interface dynamically evolves during the simulation.The initial geometry of the ellipsoidal droplet is defined with semi-axes of a = 2.5µm, b = 2.5µm, 0 .0-0 .5 0 .5 -1 .0 1 .0-1 .4 1 .4-1 .9 1 .9-2 .4 2 .4-2 .9 2 .9-3 .4 3 .4-3 .9 3 .9-4 . of the droplet.Remarkably low error levels reveal our method's proficiency in preserving the ϕ = 0 iso-surface.Notably, as time elapses, the error diminishes, reflecting the transformation of the droplet's shape from ellipsoid to sphere during oscillations.This result aligns with our earlier findings, which showcased the method's superior performance in preserving spherical interfaces compared to ellipsoidal ones.In Fig. 15(b), the distribution of normalized Chamfer distance errors across various time steps is depicted in the form of a histogram bar plot.While the maximum normalized error remains comfortably below 0.05 per cent, an important observation emerges-the first bar, representing the range from 0 to 0.005 per cent, is notably larger than the rest.This observation underscores the robustness of our approach in consistently preserving the interface throughout dynamic simulations, particularly during the level-set re-initialization process.

Droplet in Channel
The second subsection presents the simulation results of a challenging test case involving a droplet confined in a channel.The droplet's interaction with the channel walls and the effects of flow conditions make this test case particularly demanding in terms of accurately capturing the interface evolution.Confined droplets exposed to airflow are found in many important practical applications, including the flow channels of fuel cells, where complex droplet-airflow interactions are crucial for the correct evacuation of residual liquid from the device, an essential aspect of its reliable functionality [6].Additionally, this phenomenon has relevance in various other areas such as drug delivery, material synthesis, and lab-on-chip devices, where micro-scale droplet behavior in confined spaces plays a pivotal role [42,43,44].
The geometry of the test case, as illustrated in Fig. 16(a), includes a channel with a height of h = 200µm, a width of w = 300µm, and a length of l = 800µm.The initial diameter and the contact angle of the droplet are set to D 0 = 214µm, and 90 degrees, respectively.For the inlet boundary condition, a prescribed velocity profile is applied in the x-direction, described by the following equation: At the outlet, a constant pressure boundary condition (zero pressure) is enforced.Additionally, the initial position of the droplet is set to be 1.5h away from the inlet, aiming to minimize the influence of the spatially uniform velocity applied at the domain's boundary.The time step is chosen as ∆t = 10 −6 s, and the computational domain is discretized into approximately 250,000 elements.Fig. 16(b) provides a snapshot of the 3D view of the simulation results, capturing the moment when the droplet is influenced by the airflow within the channel.We conducted two simulations for this scenario.In the first simulation, the Reynolds number was set to 140 (u 0 ≈ 5 m/s).Confident in our method's capability to handle higher Reynolds numbers effectively, we performed a second simulation with intensified airflow to reach a Reynolds number of 200 (u 0 ≈ 7 m/s).This increase in airflow intensity eventually causes the droplet to detach from the substrate and begin moving within the channel.For a comprehensive view of the complete temporal evolution of these simulations, readers may refer to the animation (see the video uploaded  as supplementary material).
To further illustrate the effectiveness of our proposed method, Fig. 17 presents iso-phi contours on three distinct plane cross-sections: xy, xz, and yz.Thanks to our approach, which fully reconstructs the level-set, these contours exhibit an exceptional smoothness that allows us to successfully simulate this challenging test case.
It is also worth mentioning that the authors attempted to simulate this test case using other level-set re-initialization methods, i.e. the approach proposed by R. N. Elias, M. A. Martins, and A. L. Coutinho [33], as well as the traditional elliptic re-distancing [35].However, these simulations were unsuccessful due to irregularities in the level-set, leading to spurious droplet shape evolution.

Summary and conclusions
In this work, a two-step level-set re-initialization technique was presented.In the first step, the level-set function was reconstructed based on a Poisson equation, providing an accurate initial state for the subsequent optimization process.In the second step, the level-set function was regularized using an optimization approach, aided by a well-defined potential.
It must be noted that, in the presented approach, boundary conditions are needed only for the reconstruction step.Therefore, compared to the often used hyperbolic approaches, the proposed algorithm is less sensitive to the chosen boundary condition, and inaccuracies at the initialization step are corrected by the consecutive elliptic level-set re-initialization.
To ensure accurate re-initialization, a novel geometrical approach was employed to define boundary conditions.Its essential benefit is that it allows to effectively circumvent blind spots, which are critical in maintaining a smooth level-set function, particularly near the contact line.
Additionally, incorporating the element splitting into re-initialization plays crucial role in preserving the location of the zero level-set, which is a fundamental aspect of maintaining interface accuracy.
The proposed method is tested in various scenarios, including static and dynamic droplets, both spherical and ellipsoidal, and droplets confined in channels.The numerical tests carried out confirmed the accuracy and effectiveness of the approach in preserving the level-set function and accurately capturing interface dynamics.
Overall, even though the present work concentrated on the FEM/Levelset formulation, the proposed strategy can be beneficial for other methods that employ level-set for interface tracking.

Figure 2 :
Figure 2: Reproduction of the boundary condition for a free droplet.

Figure 3 :
Figure 3: Blind spots in the level-set domain adjacent to a droplet cut by a substrate.

Figure 4 :
Figure 4: Reproduction of the boundary condition in the blind spot adjacent to the cut interface.

Figure 5 :
Figure 5: Reproduction of the boundary condition in the blind spot on the side boundaries.

Figure 6 :
Figure 6: Generation of the sub-element entities for the calculation of the interface and contact line preserving penalty terms.In this figure, the portion of the outer boundary of the domain that coincides with element Ω e is shaded and denoted by ∂Ω e = ∂Ω ∩ Ω e .The points associated with the Gauss quadrature are also shown as red solid dots.

Figure 7 :
Figure 7: Side views of three spherical droplets on a substrate with different contact angles (Obtuse angle, Right angle, and Acute angle) and same radius R = 250µm (right column), as well as, the corresponding contours of iso-level-set (left column)

Figure 8 :
Figure 8: Evolution of the (average) L 1 -norm of the error in the gradient of the level-set function by during the iterative optimization procedure.

Figure 11 :
Figure 11: Side views of a sessile ellipsoidal droplet from three different perspectives: the xy, xz, and yz planes (upper row), as well as, the corresponding contours of iso-level-set (lower row)

Figure 14 :
Figure 14: Temporal evolution of an oscillating droplet

Figure 15 :
Figure 15: a) Time evolution of the Chamfer distance error for the oscillating droplet.The error is normalized to the mean radius of the droplet.b) Distribution of the normalized Chamfer distance error during simulation time steps

Figure 16 :
Figure 16: a) The geometry of the test case.b) Snapshot of the 3D view of the simulation results of droplet influenced by the airflow within the channel.

Figure 17 :
Figure 17: Droplet in flow channel.Simulated using the proposed elliptic distance reinitialization.