Treatments of the singularity in the level set method for unstructured meshes

Finite Elements in Analysis and Design 162 (2019) 34–50

Contents lists available at ScienceDirect

Finite Elements in Analysis and Design journal homepage: www.elsevier.com/locate/finel

Treatments of the singularity in the level set method for unstructured meshes Long Cu Ngo a, Hyoung Gwon Choi b, *, Muhammad Sufyan a a b

Department of Mechanical Engineering, Seoul National University of Science and Technology, 232 Gongneung-ro, Nowon-gu, Seoul, 01811, South Korea Department of Mechanical and Automotive Engineering, Seoul National University of Science and Technology, 232 Gongneung-ro, Nowon-gu, Seoul, 01811, South Korea

A R T I C L E I N F O

A B S T R A C T

Keywords: Level set method Curvature Singularity Finite element Unstructured mesh

Straightforward computation of normal vectors and curvature of interface is a main advantage of the level set method. However, the level set method may have discontinuities (kinks) in its derivatives that may lead to erroneous estimation of the normal vector and curvature during topological changes. In this study, we extended two existing methods of singularity treatment for the level set method on structured mesh, i.e., the Salac and Lu method (SLM) and local level set extraction method (LOLEX), to those with unstructured mesh in the framework of finite element formulation. An iterative treatment for the dependent nodes that belong to two elements cut by different interfaces was newly proposed for the LOLEX approach on unstructured mesh. Three static simulations showed that the present SLM or LOLEX for unstructured mesh accurately computed the normal vector and the curvature in the regions with singularity. Two dynamic cases were further simulated to investigate the effect of singularity-treatment on the mass conservation as well as the merging mechanism of droplets. While SLM and LOLEX provided results consistent with existing experimental results, it was confirmed that LOLEX provided a better mass conservation and more accurate curvature fields than SLM after coalescence.

1. Introduction

1.1. Singularity in the level set method

The accurate description of interfacial flows offers a particular challenge in numerical simulations of multiphase flow. Among many approaches, popular schemes are front tracking [1], the body fitted coordinate method [2], front capturing [3,4], and hybrid [5,6]. The front tracking method and body fitted coordinate method perform well considering the pressure jump across an interface; however, treating the merging/breaking of the interface is quite complicated. Compared with the body fitted and front tracking method, both the volume of fluid (VOF) [3] and the level set (LS) method [4] are easy to implement, and the breaking/merging problem of the interface is more straightforward. VOF methods conserve mass well, but their weakness is poor estimates for normal vector and curvature because of the evaluation of second-order spatial derivatives of abruptly-varying volume fractions [7]. The LS method naturally handles the topological changes in the interface due to the implicit formulation [8,9]. Furthermore, the calculations of curvature and normal vector are straightforward using the LS method. However, the LS method often suffers from poor conservation of mass.

In the LS method, the normal vector and the curvature of an interface can be computed from the derivatives of LS function. However, when two interfaces are close to each other, the discontinuity in a derivative of LS function between the interfaces affects the estimation of normal vector and curvature. Thus, the derivatives of the LS function would be erroneous at equidistant points from the interfaces, as depicted in Fig. 1, leading to an erroneous normal vector and curvature estimation. This issue was first described by Smereka [10] and several approaches have been developed to solve the problem of singularity in the LS method. 1.2. Curve fitting approaches For instance, Macklin and Lowengrub [11] introduced a quality function to detect the kink in the LS function. They proposed a direction difference scheme to treat the singularity of LS function. Later, they presented an improved method that employed curve fitting to accurately calculate the curvature [12]. In this method, the curvature was calculated at the interfaces with the use of a least squared parameterization of the interface. The curve parameterization was used to create a local LS

* Corresponding author. E-mail address: [email protected] (H.G. Choi). https://doi.org/10.1016/j.finel.2019.05.005 Received 10 January 2019; Received in revised form 27 April 2019; Accepted 21 May 2019 Available online 31 May 2019 0168-874X/© 2019 Elsevier B.V. All rights reserved.

L.C. Ngo et al.

Finite Elements in Analysis and Design 162 (2019) 34–50

Fig. 4. Description of the direct re-initialization algorithm for a 2D unstructured grid. xf denotes the foot of the perpendicular projection from a node to the interface segment, xsm the middle point of interface segment xs1 xs2 and dx the distance from a node to the interface segment.

Fig. 1. Two droplets about to merge. The dotted line on the top figure indicates the region of singularity/kink. The bottom plots the signed distance function along the x axis given by the black thick line on the top figure and the red circles represent the points where derivative of ϕ is not defined [14].

Fig. 5. Heaviside function distribution over two bodies in close proximity. The red lines denote the interface and the black lines represent the outer boundary of smoothing zone.

Fig. 2. Schematics of (a) P2P1 (quadratic interpolation (N) for velocities, linear interpolation (M) for pressure and Heaviside function) and (b) P2-iso-P1 elements (linear interpolation of LS function on each sub-element).

Fig. 3. Schematic presentation of cut element i.e. an element with interface segment.

Fig. 6. Schematic representation of a cut element with regions of ϕ < 0 (blue) and a region of ϕ > 0 (white). Case 1 and case 2 are only two possible cases with linear interpolation of ϕ.

function from which the curvature was calculated using standard discretization techniques. They extended this method to calculate the normal vector in addition to the curvature at the interface [13]. Another treatment based on curve fitting was presented by Lervag et al. [14]. They also used quality function [11] to detect a kink in the LS function. Interface parameterization was used to create a local LS field from which the curvature was calculated on the grid points with kink. Monotone cubic Hermite splines were used to parameterize the interface. The proposed method was found to be accurate in the calculation of curvature and normal vector. They also reported that the direction difference

approach was not always sufficient to accurately calculate the normal vector. These methods [12,14] gave good results in 2D but were difficult to extend to 3D simulations due to the use of curve-fitting. 1.3. Extraction based treatment In general, the singularity problem can be solved by using multiple 35

L.C. Ngo et al.

Finite Elements in Analysis and Design 162 (2019) 34–50

Fig. 7. Schematic representation of the construction of separate LS function in SLM.

Fig. 9. Construction of a local domain surrounding a given singular node.

simply separate the bodies, and then constructs their own LS field from which geometric properties can be computed. This method for singularity treatment was first reported by Salac and Lu [16] and is considered simple as compared to the previous methods [11–14]. When a grid point is adjacent to more than one interface, a weighted average of curvature from corresponding local LS functions was used. In Ref. [16], an advection equation was also solved for each local LS field. This method is easily extendable to three dimensions. Recently, Ervik et al. [17] proposed a new method to treat the singularity problem on structured mesh. They named their method as local LS extraction, or LOLEX, method. The quality function [12] was used to detect a kinky (singular) node in the LS field. For each singular node, two separate local LS fields were temporarily constructed by using the

Fig. 8. Two droplets about to merge. (a) Each drop is represented by a separate level set field and (b) the plot of level set functions along the center line in the top figure.

level sets method where each phase is represented by a single level set as in Ref. [15]. Therefore, the normal and curvature field can be determined for each interface from the LS function that represents the interface. However, this method needs to solve the transport equation for each LS function, thus leading to the increase of the computational cost. Moreover, when an interface folds into itself, the singularity may still appears. If only one single LS function is used to represent the interface, one can

36

L.C. Ngo et al.

Finite Elements in Analysis and Design 162 (2019) 34–50

Fig. 10. Schematic of LOLEX procedure on an unstructured mesh. (a) Construction of local domain; (b) extraction of bodies on the local domain; (c) reinitialization for body 1 on the local domain; and (d) reinitialization for body 2 on the local domain.

numerical simulations with droplet merging to investigate the effect of singularity treatment on the mass error of the LSM as well as the calculation of the normal vector and curvature of interface. In the next section, the numerical method for solving the incompressible Navier-Stokes equations, the LS advection equation and LS reinitialization are introduced, followed by the description of singularity treatment by SLM and LOLEX methods for unstructured grids. In the fourth section, we provide static test cases and dynamic simulation results with discussion. The paper is finalized with a conclusion in the last section.

interfaces that are adjacent to a singular node. An inherited kink from the original LS field could be removed by the reinitialization of the local LS fields. Geometric properties were estimated for each singular node using the local LS fields. Only one advection equation for a global LS field was solved, unlike for SLM.

1.4. Motivation of the present study Accurate calculation of geometric properties is particularly important in the region with singular nodes. All the existing studies of singularity treatment of the LS method (LSM) are based on the finite difference method and uniform structured grids [10–14,16]. To our knowledge, no papers have been published on the singularity treatment for the LSM using unstructured grids. We implemented the two singularity treatment methods (SLM and LOLEX) for unstructured mesh in the context of the finite element method (FEM). Reinitialization was conducted using a direct approach [18,19]; it is quite efficient to reinitialize a local LS field. No study has investigated the effect of singularity treatment on the mass conservation property of LSM except the Smereka's observation of mass error during topological change [10]. Therefore, we performed

2. Numerical method 2.1. An integrated FEM for the incompressible Navier-Stokes equations The governing equations are the incompressible Navier-Stokes equations with surface tension coupled with the advection equation of the LS variable, which can be written as follows: r⋅u ¼ 0 ;

37

(1)

L.C. Ngo et al.

Finite Elements in Analysis and Design 162 (2019) 34–50

Fig. 11. Schematic of dependent node. (a) Seclection of LS value at a dependent node located near two interfaces and (b) incorrect representation of interface (circle) due to singularity.

Fig. 12. Interface after (left) 1st and (right) 5th iterations of re-initialization. Blue solid lines and dotted lines denote the interface after the reinitialization and the desired interface, respectively. The five numbers in the figures indicate the number assigned to each dependent node.

ρðϕÞ

Du ¼  rp þ r ⋅ ½μðϕÞðru þ ruT Þ þ ρðϕÞg ; Dt

∂ϕ þ u ⋅ rϕ ¼ 0; ∂t

where Γ denotes the interface, j½ jΓ the jump value across the interface, σ the surface tension coefficient, κðϕÞ the interfacial curvature, and nðϕÞ the outward unit normal vector. As shown in Fig. 2, a P2P1 element allocates the velocity variable at both the vertex node and mid-nodes for quadratic interpolation, and the pressure variable only at the vertex nodes for linear interpolation. A P2-iso-P1 element is employed to allocate the LS variable at both the vertex node and mid-node for linear interpolation on each of four sub-elements. The unit normal vector and curvature are calculated from LS function as follows:

(2)

(3)

where u is the flow velocity vector, ϕ the LS function, ρðϕÞ the fluid density, p the pressure, μðϕÞ the fluid dynamic viscosity, and g the gravity. The continuity of velocity and the balance between the normal stress jump and the surface tension across the interface are to be satisfied as follows:
 
j½ujΓ ¼ 0;
 p~I þ μðϕÞðru þ ruT Þ
Γ ⋅ nðϕÞ ¼ σ κðϕÞnðϕÞ;

nðϕÞ ¼

(4)

38

rϕ ; jrϕj

(5)

L.C. Ngo et al.

Finite Elements in Analysis and Design 162 (2019) 34–50

8 0 > > >  

< 1 ϕ 1 πϕ HðϕÞ ¼ 1 þ þ sin > 2 ε π ε > > : 1

if ϕ <  ε ; if jϕj  ε ;

(8)

if ϕ > ε ;

where ε is an interfacial thickness and set to 2h with h being the characteristic length of a linear P2-iso-P1 element (Fig. 2). Density ρðϕÞ and dynamic viscosity μðϕÞ are interpolated by HðϕÞ as follows:



rϕ jrϕj

 (6)

In 2D, the curvature can be computed by using the standard discretization as follows ϕ2x ϕxx þ ϕ2y ϕyy þ 2ϕx ϕy ϕxy ϕxx þ ϕyy κðϕÞ ¼  1=2   32 ϕ2x þ ϕ2y þ ε ϕ2x þ ϕ2y þ ε

(9)

μðϕÞ ¼ μ2 þ ðμ1  μ2 ÞHðϕÞ;

(10)

where the subscripts 1 and 2 denote outside and inside the interface, respectively. The LS values of Eq. (9) and Eq. (10) are computed from the linear interpolation based on the P2-iso-P1 element to avoid unboundedness that can be caused by quadratic interpolation [20]. Although surface tension term has been treated by continuum surface tension (CSF) model [21] in most of papers, it can be directly computed by using explicitly constructed interfaces [22–25]. This approach could be quite effective when an interface is reconstructed by some algorithms such as SLM/LOLEX. However, the present study employs the semi-implicit consistent CSF model following the work of Choi et al. [26]. Therefore, the same linear basis function of a P2P1 element is employed for HðϕÞ of the surface tension term (see Fig. 2a) [27]. To alleviate the time step limitation due to surface tension, the semi-implicit CSF model proposed by Hysing [28] is adopted in this study. The model is based on the concept of differential geometry [29]. Detailed discretization procedures of the semi-implicit CSF model can be found in the study of Raessi et al. [30]. For temporal discretization of governing equations, the Crank-Nicolson scheme is used. For the numerical stability, the simulations with the employment of CSF model must satisfy the following condition for the time step size Δt [21]:

Fig. 13. Convergence histories for five dependent nodes for the case reported in Fig. 11.

κðϕÞ ¼  r ⋅ nðϕÞ ¼ r ⋅

ρðϕÞ ¼ ρ2 þ ðρ1  ρ2 ÞHðϕÞ;

(7)

To evaluate the first derivatives of LS function at a node (e.g. ϕx ), we first compute the first derivatives of LS function for all P1 sub-elements surrounding the node which are constants, and then take the average of those values. The second derivatives of LS function are calculated in the same way. The interface is smoothed by using the Heaviside function HðϕÞ defined as follows:

Δt  Δtc ¼

rffiffiffiffiffiffiffiffi ρh3 2πσ

(11)

where ρ is the average density of two phases.

Fig. 14. Curvature fields computed by using LOLEX method. Thick solid lines denote the interface before correction and dashed line indicates the interface on the local domain after the treatment of dependent nodes. 39

L.C. Ngo et al.

Finite Elements in Analysis and Design 162 (2019) 34–50

the shape of the interface with the velocity field obtained by solving the incompressible Navier–Stokes equations coupled with the LS variable. For the discretization of Eq. (3), the Taylor–Galerkin method [31,32] has been adopted. The Taylor–Galerkin discretization of the advection equation can be found in Ref. [33]. 2.2.2. Direct reinitialization of level set method for an unstructured triangular mesh This work uses a direct method to reinitialize the LS function for an unstructured triangular mesh. This method provides an improved mass conservation and is also efficient in terms of computational overhead [18,19]. In this approach, the interface is constructed from a collection of interface segments (Fig. 3) that are determined by the coordinates and LS values at the vertices of the edges containing intersection points. LS values at nodes near the interface are then replaced by the shortest distance from the node to the interface segments without changing its sign. The calculation of the shortest distance is based on the position of the foot of the orthogonal projection from a node onto the interface segment. Particularly, the shortest distance can be either the orthogonal distance (Fig. 4, case I) or the minimum distance between the node and two vertices of the interface segment (Fig. 4, case II). We reinitialize the Fig. 15. Schematic of two static disks within close distance to each other. Table 1 Comparison of computational cost with and without singularity treatment.

2.2. A level set method for multiphase flow simulation 2.2.1. Advection equation of the level set function The advection equation of LS function (Eq. (3)) is solved to capture

Method

No model

SLM

LOLEX

CPU time (s)

0.078125 (1.00)

0.25 (3.18)

0.2125 (2.65)

Fig. 16. Normal vector fields computed by different approaches.

Fig. 17. A comparison of curvature calculation. (a) Without singularity treatment (b) with SLM and (c) with LOLEX method. 40

L.C. Ngo et al.

Finite Elements in Analysis and Design 162 (2019) 34–50

3. Description of singularity treatments 3.1. Salac and Lu method A method for singularity treatment was introduced by Salac and Lu [16], called the Salac and Lu method, or SLM. This method scans the computational domain and extracts all the separate bodies represented by the global LS field, and then the local level set fields (ϕL ) of each body are temporarily built. Each ϕL does not include another body that can induce kinks. Thus, each ϕL can be reinitialized and the curvature can be calculated without any problem as depicted in Fig. 7. 3.1.1. Implementation of SLM for unstructured mesh  Construct local LS fields ϕL (ϕ1 and ϕ2 in case of two bodies) that individually denote a single interface/body as depicted in Fig. 8. The respective ϕL has no issues of kink.  Compute the normal vector n and curvature κ for each ϕL .  If a node is shared by the smoothing zones of two interfaces, mark the node as a singular node. A geometrical weighted average of the curvatures from ϕ1 and ϕ2 is used for the singular points. Fig. 18. Schematic of the second geometric/static test.

κ¼

LS function for nodes within five layers of elements away from the interface. For more details, see Ref. [18].

κ 1 jϕ2 j þ κ 2 jϕ1 j jϕ1 j þ jϕ2 j

(13)

2.2.3. Geometric approach for area/mass error calculation In LSM, the area (A) is usually computed via a Heaviside function defined in Eq. (8) by using the following formula: Z A¼

HðϕÞdΩ

(12)

Ω

However, it should be noted that the area computed by the smooth Heaviside function is the same as the analytic area, which is defined as the area surrounding the interface, in the case of a straight interface, but different from the analytic area in the case of a curved interface. Fig. 5 shows that the area computed via Heaviside function when two bodies are in close proximity is different from the actual area of regions enclosed by two interfaces, causing an incorrect calculation of mass distribution. This is due to the overlapping of smoothing zones around the interfaces when two bodies come close to each other. In this paper, we compute the area by using a geometric-based approach which is a summation of areas of regions occupied by bodies only (regions enclosed red lines shown in Fig. 5). For element cut by interface, we use a linear interpolation to determine the interface segment of cut element and then compute the area of body-occupied regions, i.e., regions of ϕ  0, as shown in Fig. 6.

Fig. 20. The error profiles of average values of estimated curvatures along a portion of a circle where singularity may occur with different grid resolutions for the problem of a circular disk above a rectangle.

Fig. 19. Comparison of curvature fields (a) without singularity treatment (b) with SLM and (c) with LOLEX method. The solid lines represent the interface. 41

L.C. Ngo et al.

Finite Elements in Analysis and Design 162 (2019) 34–50

cannot treat the singularity that appears when an interface folds onto itself as illustrated in Fig. 10a. In this case, SLM computes curvature using the typical LSM. According to Ervik et al. [17], this situation often occurs in numerical simulations, e.g., the coalescence of a falling drop into a pool. Thus, using the SLM approach for interface folds may lead to erroneous estimation of normal vector and curvature. 3.3. LOLEX method As described in the previous section, the implementation of SLM method from structured to unstructured meshes is straightforward. However, SLM can be problematic with interfacial folds where singularity exists in a single body. The LOLEX method [17] has been proposed to deal with this problem. In this section, the extension of the LOLEX algorithm to unstructured mesh is described for the cases with and without the dependent node as defined in Refs. [16,17]. 3.3.1. LOLEX without dependent node (a) Marking a singular node In the first step, the following quality function estimated by using the gradient of LS function is employed to mark a singular node: Qi ¼ j1  jrϕi jj

(14)

A node is empirically marked as singular if Qi > η for a given node i as suggested by Ervik et al. [17]. A big value of η leads to the increase of computational cost, while a small values might ignore some singular nodes. A value of η ¼ 0:005 is used in this paper following the work of [17].

Fig. 21. Comparison of curvature fields for circles and straight interfaces computed without singularity treatment (top), SLM method (middle) and LOLEX method (bottom).

(b) Construction of local domain  Solve the incompressible Navier-stokes equations to get an updated velocity field.  Solve the advection equation and then conduct reinitialization for each ϕL .  Construct a global LS field by taking the minimum among local LS values at each grid point, i.e., ¼ minfϕ1 ; ϕ2 g.

A local domain for each singular node is constructed by using the information of connectivity list. Fig. 9 shows the schematic of a local domain construction (layers ¼ 3) for a singular node in a 2D unstructured triangular grid. The first layer of neighboring nodes is built using neighboring node information of the given singular node. Then, the same marking procedure is used for the each node of the first layer to build the second layer of neighboring nodes. Repeating the same procedure, a local domain with multilayer is constructed for each singular node. After creating a local domain, a global LS field is copied to local domain (ΩL ), as shown in Fig. 10a.

When two interfaces have a cut element in common, numerical coalescence is achieved. After the coalescence, local LS functions are discarded, and only one global ϕ is used. 3.2. Shortcomings of SLM

(c) Identification of bodies in local domain To construct a local LS field (ϕL ) inside each body, we start at an arbitrary point inside a body (ϕ < 0) and move toward neighbor nodes until all the connected nodes within the body are found. After the scan

The idea of SLM is simple as compared with the curve-fitting methods [12,14]. It is also relatively straightforward to extend to the 3D domain. However, SLM can be regarded as a multibody approach because it

Fig. 22. Configuration of two drops in shear flow (left) and sketch of initial conditions actual domain (right). 42

L.C. Ngo et al.

Finite Elements in Analysis and Design 162 (2019) 34–50

Fig. 23. The evolution of the velocity field and interfaces using LOLEX for the case of drop merging in shear flow.

that the positive/negative level set values of cut elements are fixed [18]. Based on these interface segments, positive level set values for the nodes outside each body are calculated by the direct reinitialization method of reference [18]. After building a local LS field for each body by using reinitialization as shown in Fig. 10c and d, spatial derivatives of LS function can be calculated without kinks for each local LS field. Lastly, curvature and normal vector at a singular node are calculated by using the local LS field of which interface is closer to the node.

process is completed, all the nodes inside the body are excluded in the next scan process. This process is repeated until all the negative local LS fields are constructed. Fig. 10b shows that two separate bodies are identified after the scan processess. (d) Reinitialization and calculation of curvature based on a local level set field After marking each body in ΩL by using a negative LS field, a separate local LS field for each body needs to be reconstructed in ΩL . For each body, interface segments are determined by using the global LS field such

Remark 1. Although the LOLEX method can be used to estimate the curvature field for the case of the interface folding into itself, it fails to

43

L.C. Ngo et al.

Finite Elements in Analysis and Design 162 (2019) 34–50

Fig. 24. Mass error histories for the case of two drops in the shear flow with grid resolutions: (a) G1: 20 nodes/D (b) G2: 30 nodes/D (c) G3: 40 nodes/D and (d) G4: 50 nodes/D.


mþ1

ϕi  ϕmi

mþ1
ϕ

accurately estimate the normal and curvature field at the corner of folding interface region where the local domain contains single interface.

ε ¼



3.3.2. Treatment of dependent nodes When a singular node belongs to two elements cut by different interfaces at a given instant, one of the interfaces may not be correctly described (for instance, the circular interface in Fig. 11b) because the singular node will store only one LS value: the shorter among the distances to the two interfaces (Fig. 11a). These nodes are termed as “dependent nodes” in the literature [16,17]. The interface segments constructed from the LS values of those nodes may generate an incorrect local LS field resulting in the wrong estimation of normal and curvature. Therefore, a special treatment is necessary to determine the correct values of LS function of dependent nodes in order to provide correct interface segments used to reinitialize the local LS field in the local domain. In the LOLEX method for a structured grid, the local LS value of a dependent node was explicitly reconstructed using the geometrical approach introduced by Adalsteinsson et al. [34]. However, this approach cannot be applicable for unstructured meshes because the edges of element of unstructured meshes do not parallel to the axis of coordinate system, making the assumption in computing distance with an element is not permissible. In the present algorithm for unstructured mesh, LS values at dependent nodes are iteratively updated by using the direct reinitialization method [18,19]. The iteratively updated LS field is assumed to be convergent when the following error (ε) is less than 1% for all dependent nodes.

denotes the updated local LS values at a dependent node i where ϕmþ1 i after the (m þ 1) iteration. During the iterative procedure, either the negative or positive local LS values of each body are assumed to be accurate and thus fixed. If the dependent node has positive sign of LS value, we fix the negative local LS values of each bodies and vice versa. The interfaces defined by local LS fields after one and five iterations are shown in Fig. 12. Dashed line indicates a desired interface (ϕ ¼ 0) based on initial configuration of two interfaces, that is, circular and flat interfaces. Fig. 13 shows the convergence histories of local LS values for five dependent nodes shown in Fig. 11. LS values at dependent nodes are correct with respect to the interface closer to them. Fig. 14 compares the two LOLEX methods with and without iterative treatment of dependent nodes for unstructured grids. As shown in Fig. 14a, dependent nodes appear when two interfaces are about to merge and the erroneous curvature around these points may lead to unphysical results.

i

 100 ð%Þ

(15)

max

Remark 2. The SLM and LOLEX methods for the singularity treatment are purely dependent on geometry and independent of time step size.

3.4. Extension to three dimensions Though we present the singularity treatment process for the twodimension unstructured meshes, the proposed algorithm can be 44

L.C. Ngo et al.

Finite Elements in Analysis and Design 162 (2019) 34–50

directly extended to the three dimensions with no more complexity. The processes of determining singularity nodes, constructing local domains and identification of bodies in local domain are similar for both two- and three-dimensional cases. The reinitialization of level set function for local domain and the treatment of dependent nodes in the three-dimensional case are carried out by using a direct reinitialization technique for 3D unstructured meshes proposed in Refs. [18,19]. Finally, normal and curvature for singular nodes can be obtained from the local level set

Fig. 25. Grid resolution versus maximum mass error for the case of shear merging.

Fig. 27. Schematic of initial configuration (left) and boundary conditions (right) for droplet merging in a pool.

Fig. 26. Comparison of curvature fields at some selected instants. 45

L.C. Ngo et al.

Finite Elements in Analysis and Design 162 (2019) 34–50

Fig. 28. Effect of singularity treatment on curvature field (left) and pressure field (right) after one time step of the simulation. Interface from the simulation without singularity treatment is portrayed with dashed lines. The solid line is used to represent the interface from simulation with SLM treatment.

Fig. 29. A comparison of experimental results [37] and the present numerical results (with SLM). Snapshots were taken 542 μs apart. The arrow in experimental results indicates the capillary wave and the horizontal lines highlights the top of the droplet in the first frame.

Fig. 30. A comparison of experimental results [37] and the present numerical results (with LOLEX).

function in the same manner used in the two-dimensional case.

singularity treatments with respect to normal vectors, curvature, and pressure fields. We provide the detail investigation of the effects of singularity treatment on normal and curvature computation as well as the additional cost due to the singularity treatment in the static tests. In the case of dynamic simulations, the effect of singularity treatment on mass conservation was also examined.

4. Results and discussion In this section, some static and dynamic problems of bubbles/drops are investigated to compare numerical results with and without 46

L.C. Ngo et al.

Finite Elements in Analysis and Design 162 (2019) 34–50

Fig. 31. Comparison of curvature fields. (a) Before coalescence (b) after coalescence. SLM fails in the accurate estimation of curvature in after coalescence as indicated by the arrow in (b).

Fig. 32. Comparison of mass error evolutions. (a) Coarse mesh (50 nodes/D) (b) Medium mesh (60 nodes/D) and (c) Fine mesh (70 nodes/D).

singularity treatment. According to Lervag et al. [14], grid refinement without singularity treatment increased the error in curvature calculations. To examine the additional computational cost due to singularity treatment, we measure the CPU time spent for SLM (Section 3.1), LOLEX (Section 3.2) algorithms, which mainly consist of the reinitialization for each body and curvature computation. Then, we compare them with the CPU time required to compute a curvature field without singularity treatment. Table 1 provides the CPU times and the normalized computational cost (the numbers in the parenthesis) by the CPU time of no model. It shows that the CPU times with SLM and LOLEX are 3.18 and 2.65 times bigger than that without singularity treatment. However, it should be noted that the percentage of singular nodes among the total number of nodes for which curvature need to be computed is about 6.2% and that the cost for reinitialization occupies only a small portion for twodimensional simulations of dynamic problems as reported in Ref. [18]. Therefore, the cost for the singularity treatment can be considered insignificant for the dynamic simulations.

4.1. Two disks within close distance to each other A benchmark problem of two disks within close proximity (Fig. 15) was tested to examine the effect of the singularity treatment on the accuracy of normal vector/curvature fields. The size of the computational domain is 1:5 m  1:5 m and an unstructured grid with 50 nodes along the x and y direction was used. Two disks of radius r ¼ 0:25 m were initially placed at a distance of h with h being the characteristic length of a P2-iso-P1 element. Lervag et al. [14] reported that the singularity treatment based on direction difference approach [11] yielded an erroneous normal vector field for this case. Fig. 16 compares normal vector fields calculated without singular treatment and with the two singular treatments, SLM and LOLEX methods. The normal vector field without singularity treatment is erroneous in the kink region. However, both methods with singularity treatments accurately calculated the normal vector in this kink region. Fig. 17 compares curvature fields calculated with SLM, LOLEX, and without any treatment. Computing curvature is necessary only within the smooth zone of interface because the CSF model is considered only in the smooth zone. Fig. 17 demonstrates that the curvature calculation without singular treatment yields an inaccurate result in the kink region whereas both singularity treatments provide smooth curvature fields around the interface (see Fig. 17). We also confirmed that the increase in grid resolution did not improve the accuracy of curvature field without

4.2. A circular disk above a rectangle A benchmark problem of a circular disk of radius r ¼ 0:25 m positioned slightly above a rectangle at a distance of 1:5h was considered as shown in Fig. 18. When the gap between disk and 47

L.C. Ngo et al.

Finite Elements in Analysis and Design 162 (2019) 34–50

component varies linearly by moving top (u ¼ U) and bottom (u ¼  U) walls as sketched in Fig. 22. This problem was tested by using four different grids (G1: 120  80, G2: 180  120, G3: 240  160, and G4: 300  200). The shear flow is defined by the Reynolds number (Re) and capillary number (Ca) [14]. In the present study, Re of 10 and Ca of 0.025 were imposed. Under these physical conditions, drops are moderately deformed. The unit density and viscosity ratio were used, and no-slip boundary conditions were imposed on the walls. The time step was chosen as half of the critical time step Δt c computed by Eq. (11). Fig. 23 shows the evolution of interface and velocity field with singularity treatment. The effect of singularity treatment on mass error was examined. The mass/area error (%) is defined as follows:

rectangle is small enough such that the smoothing zones of interfaces are in contact, the kink along the dashed line in Fig. 18 affects the curvature estimation. The size of the computational domain is 1:5 m  1:5 m, and the surface of the rectangle is positioned at y ¼ 0:75 m. An unstructured grid with 100 nodes along the x and y directions was used. Fig. 19 compares curvature fields calculated with SLM and LOLEX and without singularity treatment. The figure indicates that the approach without singularity treatment provides inaccurate curvature along the kink whereas SLM and LOLEX methods provide smooth curvature fields (see Fig. 19b–c). In order to investigate the effect of grid size on the estimated curvature along the interface with and without singularity treatment, we test the static problem by using four different unstructured grids with 50, 100, 200 and 400 nodes in each direction of x and y. The curvature error (ε) is defined as follows:

ε¼

κext  κ  100 ð%Þ κ ext

E¼

(16)

AðtÞ  Ai  100 Ai

(17)

where Ai and AðtÞ represent the initial area and the area of drop at time t, respectively. Mass error histories with singularity treatment were nearly similar to those without singularity treatments until the instant marked with a circle in Fig. 24a, at which smoothing zones of two interfaces are within the distance of one element. However, mass errors with singular treatment were smaller than those without the treatment for all grid resolutions after collision. Fig. 24 shows that mass error also depends on grid resolution and the maximum mass error during the entire simulation was linearly proportional to h (the size of a finite element) with the present LOLEX method (Fig. 25). Fig. 26 compares the interface evolutions and the curvature calculations without singularity treatment, and with the SLM and LOLEX methods. It shows that the interface evolutions are quite similar until t ¼ 2:2 s. However, curvature field without the treatment is quite different from those with the treatment near the merging instant. Consequently, the method without the treatment provides a completely different pattern of droplets coalescence at t ¼ 3:2 s. The incorrect curvature estimation without singularity treatment led to earlier/unphysical coalescence of the interfaces.

where κext and κ are the exact and average values of estimated curvatures along a portion of the circle which is near the straight line and within 0:6  x  0:9 where singularity may occur. Fig. 20 plots error (ε) versus grid size. The curvature obtained by LOLEX method is nearly constant and equal to exact curvature. On the other hand, the curvature by SLM tends to approach to the exact one when grid resolution increases because SLM uses the geometrical weighted average scheme to compute the curvature at singular nodes and does not employ the treatment of dependent nodes of LOLEX method. The curvature obtained without singularity treatment deviates far from the exact one due to singularity and cannot be corrected by increasing the grid resolution. 4.3. Circles with a straight interface We examined a complicated static case with the interfaces of three circles and a straight-line, where two circles and a straight-line are joined. The size of the computational domain is 1:5 m  1:5 m and the same grid resolution as that of Section 4.2 is employed. The straight line is positioned at y ¼ 0:75 m and the radius of circles is r ¼ 0:2 m. The circles from left to right are centered at ð0:282; 0:93Þ, ð0:65; 0:972Þ and ð0:075; 0:972Þ. Comparison of curvature fields is shown in Fig. 21 for LOLEX, SLM, and without singularity treatment. The figure shows that the method without singularity treatment produces unphysical curvatures at several regions; the bottom of the second and third circles and in the middle of second and third circles. In the case of SLM, curvature distribution is erroneous at the bottom of the second circle where the singular nodes exist between the second circle and the straight line, which are joined as one body. On the other hand, the present LOLEX method produces negative curvatures only where they should be. However, in regions where the two circles join, all the methods fails to correctly estimate the curvature field. In this case, the SLM method cannot be applicable since there exists only one body. The LOLEX method can treat the problem in regions such as the bottom part of the middle circle, but fails to handle problem at the junctions of the two circles. This is due to the local domain building from a singular node in those regions contains only one body, and the normal and curvature are computed as similar to the ordinal method.

4.5. Partial coalescence of water drop into pool Charles and Mason [35] found the phenomenon of “partial coalescence cascade” in the experiment of a drop falling onto a fluid-fluid interface. They claimed that the effect of capillary pinch-off was the reason of partial coalescence. The phenomenon is difficult to simulate because it is necessary to resolve both very small time and length scales [36]. In this section, two immiscible liquids are considered with a droplet of one liquid (water) being placed above the pool of the other liquid. In the experimental work by Chen et al. [37], the droplet was made to rest on the pool and merging occurred after some time. The liquid pool surrounding the droplet is a mix of 20% polybutene in decane. The initial droplet diameter D is 1:1 mm. To be consistent with existing studies of this partial coalescence [17,36], we employed an axisymmetric geometry as shown in Fig. 27. The rectangular computational domain has H ¼ 7D and R ¼ 3D. The depth of water H1 ¼ 3D, and we initially placed the drop very close to the flat interface (with the gap of 2h, where h ¼ D=70). The dimensionless numbers governing the flow for this problem are the Bond number (Bo) and the Ohnesorge (Oh) number defined as follows:

ρ2 gr 2 μ1 ffi : ; Oh ¼ pffiffiffiffiffiffiffiffiffiffiffi σ 2ρ1 σ r

4.4. Drop collision in shear flow

Bo ¼

Merging of two drops in shear flow was considered. This case was chosen to investigate the accuracy of curvature calculation in dynamic simulation with and without singularity. Two drops of radius R ¼ 0:5 m are initially placed at a distance of d ¼ 5R as shown in Fig. 22. The physical domain size is 12R  8R. The distance s between centerline and drop center is 0:84R, that is, 0:42 m. The initial velocity profile of the x-

Here, the value of the Bond number (Bo) is 9:58  102 and the Ohnesorge number (Oh) is 5:53  103 . The physical parameters of the two fluid are: ρ1 ¼ 1; 000 kg ⋅ m3 , μ1 ¼ 0:001 kg ⋅ m1 ⋅ s1 and ρ1 ¼ 760 kg ⋅ m3 , μ1 ¼ 0:002 kg ⋅ m1 ⋅1 , respectively. The time step was chosen as ⋅t  0:5 ⋅ t c ¼ 2 ⋅ 106 s. Fig. 28 shows the effect of singularity 48

(18)

L.C. Ngo et al.

Finite Elements in Analysis and Design 162 (2019) 34–50

treatment on pressure and curvature fields. Without singularity treatment, the curvature field near the gap was erroneous whereas the present method of singularity treatment provided a smoother curvature field. Fig. 28b shows the effect of erroneous curvature estimation on pressure field. Unlike the pressure field with SLM method, the pressure field without singularity treatment appears unphysical because of the erroneous estimation of surface tension. Because SLM and LOLEX produced almost identical results, only results of SLM were presented. Fig. 29 and Fig. 30 compare the experimental visualization and the present numerical results of interface evolution. Initially, while draining the drop into the liquid pool, a capillary wave (shown by an arrow in the experimental pictures) is generated. The drop then evolves into a cylindrical body of liquid because of the upward momentum of the capillary wave. Finally, pinch-off occurs with a daughter droplet suspended above the deformed liquid pool at the end of simulation. Results from SLM and LOLEX methods were both in a good qualitative agreement with the experimental data. In the present simulations, the inability of SLM to accurately estimate the curvature after coalescence (Fig. 31) affected the solution a little. Ervik et al. [17] reported that some of their simulation using SLM crashed after coalescence. Lastly, we compared the mass error histories from numerical simulations by SLM and LOLEX. The mass error (loss/gain) is defined as follows.

ε¼

MðtÞ  Mi  100 Mdi

Acknowledgement This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIP) (No.NRF2017R1A2A2A05001177). Many comments from Professor Sanghun Choi were greatly appreciated. Appendix A. Supplementary data Supplementary data to this article can be found online at https:// doi.org/10.1016/j.finel.2019.05.005. References [1] S.O. Unverdi, G. Tryggvason, A front-tracking method for viscous, incompressible, multi-fluid flows, J. Comput. Phys. 100 (1992) 25–37. [2] G. Ryskin, L.G. Leal, Numerical solution of free-boundary problems in fluid mechanics. Part 1. The finite-difference technique, J. Fluid Mech. 148 (1984) 1. [3] C.W. Hirt, B.D. Nichols, Volume of fluid (VOF) method for the dynamics of free boundaries, J. Comput. Phys. 39 (1981) 201–225. [4] S. Osher, J.A. Sethian, Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations, J. Comput. Phys. 79 (1988) 12–49. [5] M. Sussman, E.G. Puckett, A Coupled level set and volume-of-fluid method for computing 3D and axisymmetric incompressible two-phase flows, J. Comput. Phys. 162 (2000) 301–337. [6] S. Shin, D. Juric, Modeling three-dimensional multiphase flow using a level contour reconstruction method for front tracking without connectivity, J. Comput. Phys. 180 (2002) 427–470. [7] S.J. Cummins, M.M. Francois, D.B. Kothe, Estimating curvature from volume fractions, Comput. Struct. 83 (2005) 425–434. [8] M. Sussman, E. Fatemi, P. Smereka, S. Osher, An improved level set method for incompressible two-phase flows, Comput. Fluids 27 (1998) 663–680. [9] S. Osher, R.P. Fedkiw, Level set methods: an Overview and some recent results, J. Comput. Phys. 169 (2001) 463–502. [10] P. Smereka, Semi-implicit level set methods for curvature and surface diffusion motion 19 (2003) 439–456. [11] P. Macklin, J. Lowengrub, Evolving interfaces via gradients of geometry-dependent interior Poisson problems: application to tumor growth, J. Comput. Phys. 203 (2005) 191–220. [12] P. Macklin, J. Lowengrub, An improved geometry-aware curvature discretization for level set methods: application to tumor growth, J. Comput. Phys. 215 (2006) 392–401. [13] P. MacKlin, J.S. Lowengrub, A new ghost cell/level set method for moving boundary problems: application to tumor growth, J. Sci. Comput. 35 (2008) 266–299. [14] K.Y. Lervåg, B. Müller, S.T. Munkejord, Calculation of the interface curvature and normal vector with the level-set method, Comput. Fluids 84 (2013) 218–230. [15] B. Scholtes, M. Shakoor, A. Settefrati, P.-O. Bouchard, N. Bozzolo, M. Bernacki, New finite element developments for the full field modeling of microstructural evolutions using the level-set method, Comput. Mater. Sci. 109 (2015) 388–398. [16] D. Salac, W. Lu, A local eemi-implicit level-set method for interface motion, J. Sci. Comput. 35 (2008) 330–349. [17] Å. Ervik, K.Y. Lervåg, S.T. Munkejord, A robust method for calculating interface curvature and normal vectors using an extracted local level set, J. Comput. Phys. 257 (2014) 259–277. [18] L.C. Ngo, H.G. Choi, Efficient direct re-initialization approach of a level set method for unstructured meshes, Comput. Fluids 154 (2017) 167–183. [19] M. Shakoor, B. Scholtes, P.-O. Bouchard, M. Bernacki, An efficient and parallel level set reinitialization method – application to micromechanics and microstructural evolutions, Appl. Math. Model. 39 (2015) 7291–7302. [20] M.H. Cho, H.G. Choi, S.H. Choi, J.Y. Yoo, A Q2Q1 finite element/level-set method for simulating two-phase flows with surface tension, Int. J. Numer. Methods Fluids 70 (2012) 468–492. [21] J.U. Brackbill, D.B. Kothe, C. Zemach, A continuum method for modeling surface tension, J. Comput. Phys. 100 (1992) 335–354. [22] E. Bӓnsch, P. Morin, R.H. Nochetto, A finite element method for surface diffusion: the parametric case, J. Comput. Phys. 203 (2005) 321–343. [23] G.C. Buscaglia, R.F. Ausas, Variational formulations for surface tension, capillary and wetting, Comput. Methods Appl. Mech. Eng. 200 (2011) 3011–3025. [24] D.P. Mu~ noz, J. Bruchon, S. Drapier, F. Valdivieso, A finite element-based level set method for fluid-elastic solid interaction with surface tension, Int. J. Numer. Methods Eng. 93 (2013) 919–941. [25] J. Bruchon, Y. Liu, N. Moulin, Finite element setting for fluid flow simulations with natural enforcement of the triple junction equilibrium, Comput. Fluids 171 (2018) 103–121. [26] S. Choi, M.H. Cho, H.G. Choi, J.Y. Yoo, A Q2Q1 integrated finite element method with the semi-implicit consistent CSF for solving incompressible two-phase flows with surface tension effect, Int. J. Numer. Methods Fluids 81 (2016) 284–308. [27] S. Zahedi, M. Kronbichler, G. Kreiss, Spurious currents in finite element based level set methods for two-phase flow, Int. J. Numer. Methods Fluids 69 (2012) 1433–1456.

(19)

where Mdi , Mi and MðtÞ denote the initial mass of the drop and the total initial mass of water and the total mass of water with time, respectively. Fig. 32 presents the mass error histories from numerical simulations with SLM and LOLEX for the three different unstructured meshes: G1 (50 nodes/D), G2 (60 nodes/D), and G3 (70 nodes/D). As the mesh resolution increased, smoother evolutions of mass errors were provided. Both SLM and LOLEX seem to conserve mass well for both the base and finer grids. Mass error is normalized with the initial mass of the drop instead of the total initial mass of water. 4.5.1. Limitation of the present simulation The film drainage theory [38] reported that the process of droplet merging consisted of three stages: attachment, drainage of trapped film of the continuous phase, and rupture and coalescence. In this study, an event at a very small length scale was not considered because the computational cost to resolve the whole microscopic process of the film drainage is prohibitive as a result of the adoption of an extremely fine mesh. 5. Conclusion We developed algorithms of two singularity treatments (SLM and LOLEX) for unstructured mesh by constructing a separate level set field for each body. An iterative method for constructing a local level set field on unstructured mesh was proposed for the present LOLEX method. Three static cases with given analytical curvature fields were tested. Results showed that both the present SLM and LOLEX methods provided more accurate curvature fields than the approach without singularity treatment when singular nodes existed. We also confirmed that unlike SLM, the present LOLEX method could remove an erroneous curvature field of a singular region of two joined interfaces that were identified as parts of a single body. In dynamic test cases, we found that the present treatments of singularity not only accurately estimated the normal vector and curvature fields but also improved the mass conservation. Simulations of partial coalescence of the drop into a pool agreed well with the existing experimental result. While the present SLM and LOLEX methods provided consistent results with the existing experimental results, the present LOLEX provided more accurate curvature fields than the SLM after coalescence. 49

L.C. Ngo et al.

Finite Elements in Analysis and Design 162 (2019) 34–50

[28] S. Hysing, A new implicit surface tension implementation for interfacial flows, Int. J. Numer. Methods Fluids 51 (2006) 659–672. [29] S. Gallot, D. Hulin, J. Lafontaine, Differential Manifolds (2004) 1–49. [30] M. Raessi, M. Bussmann, J. Mostaghimi, A semi-implicit finite volume implementation of the CSF method for treating surface tension in interfacial flows, Int. J. Numer. Methods Fluids 59 (2009) 1093–1110. [31] J. Donea, A Taylor–Galerkin method for convective transport problems, Int. J. Numer. Methods Eng. 20 (1984) 101–119. [32] A.J. Baker, J.W. Kim, A Taylor weak-statement algorithm for hyperbolic conservation laws, Int. J. Numer. Methods Fluids 7 (1987) 489–520. [33] M.H. Cho, H.G. Choi, J.Y. Yoo, A direct reinitialization approach of level-set/ splitting finite element method for simulating incompressible two-phase flows, Int. J. Numer. Methods Fluids 67 (2011) 1637–1654.

[34] D. Adalsteinsson, J. Sethian, The fast construction of extension velocities in level set methods, J. Comput. Phys. 148 (1999) 2–22. [35] G. Charles, S. Mason, The mechanism of partial coalescence of liquid drops at liquid/liquid interfaces, J. Colloid Sci. 15 (1960) 105–122. [36] P. Yue, C. Zhou, J.J. Feng, A computational study of the coalescence between a drop and an interface in Newtonian and viscoelastic fluids, Phys. Fluids 18 (2006). [37] X. Chen, S. Mandre, J.J. Feng, Partial coalescence between a drop and a liquidliquid interface, Phys. Fluids 18 (2006) 16–20, https://doi.org/10.1063/ 1.2201470. [38] A.K. Chesters, The modelling of coalescence processes in fluid-liquid dispersions : a review of current understanding, Trans IChemE 69 (1991) 259–270.

50