A three-dimensional volume-of-fluid method for reconstructing and advecting three-material interfaces forming contact lines

Journal of Computational Physics 307 (2016) 550–573

Contents lists available at ScienceDirect

Journal of Computational Physics www.elsevier.com/locate/jcp

A three-dimensional volume-of-fluid method for reconstructing and advecting three-material interfaces forming contact lines Ashish Pathak, Mehdi Raessi ∗ Department of Mechanical Engineering, University of Massachusetts Dartmouth, North Dartmouth, MA 02747, United States

a r t i c l e

i n f o

Article history: Received 19 May 2015 Received in revised form 15 November 2015 Accepted 18 November 2015 Available online 11 December 2015 Keywords: Contact line Contact angle Three-material Volume-of-fluid PLIC

a b s t r a c t We introduce a piecewise-linear, volume-of-fluid method for reconstructing and advecting three-dimensional interfaces and contact lines formed by three materials. The new method employs a set of geometric constructs that can be used in conjunction with any volumetracking scheme. In this work, we used the mass-conserving scheme of Youngs to handle two-material cells, perform interface reconstruction in three-material cells, and resolve the contact line. The only information required by the method is the available volume fraction field. Although the proposed method is order dependent and requires a priori information on material ordering, it is suitable for typical contact line applications, where the material representing the contact surface is always known. Following the reconstruction of the contact surface, to compute the interface orientation in a three-material cell, the proposed method minimizes an error function that is based on volume fraction distribution around that cell. As an option, the minimization procedure also allows the user to impose a contact angle. Performance of the proposed method is assessed via both static and advection test cases. The tests show that the new method preserves the accuracy and mass-conserving property of the Youngs method in volume-tracking three materials. © 2015 Elsevier Inc. All rights reserved.

1. Introduction In multi-phase flow simulations, the volume-of-fluid (VOF) method is one of the most commonly used interface capturing techniques due to (a) its mass conserving property, and (b) ease in handling large interfacial deformations. In the VOF method, an indicator function ψ defined as,



ψ(x) =

1, 0,

x ∈ Material 1 x ∈ / Material 1

(1)

is used to represent Material 1 shown in Fig. 1. In the discretized form, the VOF function is the fraction of cell volume V occupied by Material 1,

=

1



ψ dV

V V

*

Corresponding author. E-mail addresses: [email protected] (A. Pathak), [email protected] (M. Raessi).

http://dx.doi.org/10.1016/j.jcp.2015.11.062 0021-9991/© 2015 Elsevier Inc. All rights reserved.

(2)

A. Pathak, M. Raessi / Journal of Computational Physics 307 (2016) 550–573

551

Fig. 1. (a) Interface between two materials. (b) PLIC approximation as linear (planar in 3D) segments of the interface shown in (a).

Fig. 2. (a) Three-material configuration where the interface (I 23 ) between materials 2 and 3 is truncated by Material 1 interface (I 1 ), resulting in a triple point T (contact line in 3D). (b) PLIC representation of interfaces around the triple point bounded by S in (a).

A common approach in the VOF method is to use Piecewise-Linear-Interface-Calculation (PLIC) to reconstruct the material interfaces. Fig. 1 shows the interface approximated as linear (planar in 3D) segments. The numbers shown in Fig. 1(b) correspond to the VOF function F . There exist various volume-tracking methods that rely on PLIC. An example is the Youngs method [1], which has been used in this work. In the Youngs method, to reconstruct the interface, the orientation of each ∇ψ segment is computed by using nˆ ψ = |∇ψ| [2, p. 37], while its location is determined by matching exactly Material 1’s volume fraction. Now, consider the three-material configuration shown in Fig. 2. We denote by Material 1 (or M1), the material that shares a continuous interface with the other two materials, M2 and M3. For example, in contact line problems, M1 represents the material that forms the contact surface; e.g., a solid surface on which a droplet resides. The interface between materials M2 and M3 intersects with the Material 1 surface resulting in a triple point in 2D and a contact line in 3D (Fig. 2). We have developed a method in the PLIC-VOF context to reconstruct three-material configurations in 3D, capture contact lines and track volumes of three materials. A similar approach based on error minimization was proposed by Choi and Bussmann [3] to resolve three-material cells and the triple point in 2D. The method was later used by Ghasemi et al. [4] to resolve triple points in their study of the interaction between two-phase fluid flows and a moving rigid solid body. Choi and Bussmann [3] noted that their procedure could be extended to 3D conceptually, but was “significantly more complicated to implement” owing to “much more complicated geometry” [3, p. 1008]. We have presented in this paper a complete set of 3D geometrical constructs needed for the three-material reconstruction method. The proposed geometric constructs, detailed in Appendices and the Supporting Information document, are based on convex 3D polyhedra and are quite general. They are applicable to both structured and unstructured grids. The method presented here is not a direct 3D extension of Choi and

552

A. Pathak, M. Raessi / Journal of Computational Physics 307 (2016) 550–573

Bussmann’s [3]. It is different from their method both in the definition of the minimization function and in the geometric constructs. We will now describe previous efforts toward resolving interfaces in multi-material cells in the context of VOF. “Onion skin method” [5,6] can be identified as the simplest approach belonging to this category. In this approach, each interface in a multi-material cell is assumed to separate two materials and not to become truncated due to the presence of another material. The scheme is, therefore, unable to resolve contact lines, and is suitable only for layered configurations. In layered configurations, material ordering becomes important, which was addressed by Mosso and Clancy [7] and Benson [8]. They determined the dynamic ordering according to the relative positions of the material centroids “fitted” along a line. More recently, Niem et al. [9] performed material ordering following a strategy based on intersection checks. Sijoy and Chaturvedi [10] provided a nice comparison of different material ordering schemes in the VOF context. They also proposed a new method by combining two earlier schemes. Caboussat et al. [11] formulated a 2D approach that located the triple point by solving a constrained minimization problem inside three-material cells. They switched to the onion-skin method of Youngs [5,6] when the triple point was estimated to lie outside the cell’s bounds. Performance of the method in resolving the Y-shaped interfaces seemed to be better than in resolving the T-shaped ones. The method is material order independent i.e., a priori information on the material ordering is not required. Schofield et al. [12] proposed a 2D method that estimated the centroid location of each material in multi-material cells using a linear approximation of the indicator function f . Power diagrams were then used to locate the interfaces between each pair of materials by matching exactly their volume fractions inside the cell. Interfaces were smoothed to achieve second-order accuracy. Schofield et al. [12] assessed various PLIC reconstruction techniques and noted that LVIRA [13] failed to yield correct interface orientation in multi-material cells because the error function calculation did not acknowledge the presence of an additional (third) material behind the extended interface line (plane in 3D). As will be presented later, this idea serves as the motivation for our three-material reconstruction algorithm: if the presence of the additional material is indeed taken into account while calculating the error function, LVIRA can be augmented to resolve correctly the three-material cells. Braeunig et al. [14] developed a method called Natural Interface Positioning (NIP), which was further improved in [15]. Here, the interfaces were assumed to be mono-dimensional, but multi-dimensionality was taken into account by considering interface normals. “Condensates” were created by joining contiguous material zones, and layers were advected in a Lagrangian manner. After the interface positions were determined, they were mapped back to the Eulerian grid. Normal ∇f vectors were calculated using the gradient of the indicator function f (nˆ = |∇ f | ). Blais et al. [16] extended the NIP method to handle more than two materials. Onion skin approach [5] was adopted, where normals were calculated using a material priority list provided by the user. To avoid interface intersections inside mixed cells, they assumed the same orientation for all interfaces. Automatic material ordering was achieved by the approximate material centroid locations available from the previous time-step. The Moment-of-Fluid (MoF) method, as presented by Ahn and Shashkov [17], is another approach to reconstructing and tracking multi-material interfaces. Ahn and Shashkov [17] provided a nice comparison of three reconstruction techniques in 3D: (1) Gradient based interface reconstruction (GRAD), (2) LVIRA and (3) MoF. MoF method uses the additional information on centroid position of each material to perform reconstruction in a multi-material cell. The interface orientation is determined by minimizing an error function, which is constructed as the difference between the prescribed (target) and the predicted material centroid locations. Matching exactly each material’s volume fraction acts as the constraint to the minimization problem. Nested dissection approach makes the method capable of automatic material ordering. Dyadechko and Shashkov [18] suggested a scheme called “Lagrangian remap” that calculated and advected cellwise material centroids making MoF suitable for fluid flow simulations in an Eulerian grid. The “Lagrangian remap” was assessed against 2D twomaterial test cases in [19]. More recent adaptations of MoF have been mostly in the Arbitrary-Eulerian–Lagrangian (ALE) framework [20–22], where the computational volumes followed the material motion during the Lagrangian stage. These later adaptations of MoF were tested against 2D advection problems that included triple point configurations. The visualization community has also been involved in multi-material reconstruction. Bonnell et al. [23] used a dual mesh, the vertices of which were associated with the tuples constructed from the material volume fractions. The dual mesh cells were then intersected with Voronoi cells in barycentric space. The intersection points were mapped back to the physical space to determine interface positions. Another method, based on volume fraction averaging, was proposed by Meredith [24]. Bilinear interpolations using the averaged nodal values were carried out to construct volume fraction interpolations for each material. Looking at each point inside the cell for the material whose interpolated volume fraction dominates, one can reconstruct multiple interfaces. The method was order independent. Anderson et al. [25] proposed a method for interface reconstruction based on energy minimization using simulated annealing. A mixed cell was discretized into a number of subcells. Labeling on the subcells was allowed to change in order to achieve a minimum energy state. Of all the schemes discussed so far, only multi-material reconstructions presented in [17] and [23–25] were performed in 3D. The methods presented in [23–25] did not preserve volumes exactly, however, the volume error in [25] was the lowest and could be improved by discretizing a mixed cell into more number of subcells. More subcells would, however, amount to slower simulations. In contrast to MoF [17], the 3D reconstruction algorithm presented in this paper relies only on the available volume fraction information to calculate locations and orientations of the interfaces present in a three-material cell. We do not advect additional quantities like cellwise material centroids. The material volume fraction advection, carried out

A. Pathak, M. Raessi / Journal of Computational Physics 307 (2016) 550–573

553

with the “operator split” scheme [26], remains the same as in the two-material configuration; the three-material algorithm is only invoked during the reconstruction stage. While the focus of this paper is on the volume-of-fluid method, we briefly mention some of the approaches in the context of the level set method for resolving multi-material configurations and contact lines. To resolve the contact line formed by liquid-gas-solid phases, [27] and [28] extend the liquid-gas interface into the solid domain. Then, an additional level set equation, which imposes a contact angle, is solved in a small region around the contact line. An elegant approach was proposed by Saye and Sethian [29,30], where an unsigned level set function is advanced in time. Then, at the new time,  level sets are used to reconstruct “Voronoi interface” using the Eikonal reinitialization algorithm. The Voronoi interface represents the actual interface separating multiple phases. The three-material interface reconstruction method presented in this paper is order dependent and requires a priori information on material ordering. However, it is suitable for the contact line applications because the material representing the contact surface is always known. Although inspired by LVIRA, our algorithm does not rely on the iterative scheme of LVIRA (which could become computationally expensive as noted by [17]) for reconstruction in two-material cells. Instead, ∇f the fast method of Youngs [1] is used there, where the interface orientation is calculated by nˆ f = |∇ f | (see Ref. [2, p. 37] for details). However, for accurate interface reconstruction in two-material cells located near the contact line (e.g., cell r in Fig. 2(b)), for which the volume fraction field is truncated by the third material, the truncated interface is extended into the third material and a virtual volume fraction field is generated and used in the above formula. The structure of this paper is as follows: We describe our algorithm in Section 2 providing the mathematical formulation of the problem. Section 3 assesses the performance of the proposed reconstruction using both static and advection test cases. Finally, our work is summarized in Section 4. The 3D geometric constructs have been presented in the Appendices. The pseudo-codes corresponding to the geometric constructs are included in the Supporting Information document. 2. Resolving three material cells Our description henceforth will be in a 3D context. The figures used to explain our algorithm are mostly 2D for clarity and convenient representation. Consider a three-material configuration as depicted in Fig. 2(a). We first denote by M1 the material whose interface is untuncated. To be succinct, we denote by I 1 the untruncated interface of M1, and by I 23 the interface between M2 and M3 as shown in Fig. 2(a). In contact line problems, the triple point, denoted by T in Fig. 2(a), will lie on I 1 . Note that in 3D computational cells, both I 1 and I 23 are represented by 3D planar polygons. We also denote the color functions representing M1 and M2 by ψ and f , respectively, where



ψ(x) =  f (x) =

1, 0,

x ∈ M1 x ∈ / M1

(3)

1, 0,

x ∈ M2 x ∈ / M2

(4)

In the discretized form, ψ and f assume a value between 0 and 1 inside the cells containing portions of I 1 and I 23 , respectively. The volume fraction of M3 can then be given as α = 1 − ψ − f . I 1 can be easily reconstructed with the method ∇ψ of Youngs [1] by using the orientation calculated from nˆ ψ = |∇ψ| , which was designed for two-phase applications, owing to I 1 ’s continuity in the neighborhood of three-material cells. The orientation of the interface I 23 in three-material cells is estimated using a minimization procedure inspired from LVIRA [13]. Our algorithm would require the knowledge of material M1, which we reconstruct first, a priori and is order dependent. 2.1. Reconstructing I 23 in three-material cells Main focus of the present work is the geometrical reconstruction of the interface I 23 in three-material cells. A precise geometrical reconstruction would allow for an accurate transport of mass (and momentum) in three-material cells and their neighborhood. As shown in Fig. 3, a three-material cell can fall in one of the following three categories: (a) interface I 1 lying completely “below” the plane of I 23 , (b) the two interfaces intersecting and (c) interface I 1 lying “above” the plane of I 23 . The algorithm presented in this paper can reconstruct I 23 in all three scenarios. We have introduced the terms “above” and “below” a plane to facilitate the ensuing discussion. Consider Fig. 3(a) where x f lies on the plane of I 23 and x is an arbitrary point in space. We calculate the quantity d = (x − x f ) · nˆ f , where nˆ f is normal to the plane I 23 . d > 0, d < 0 and d = 0 mean, respectively, that the point x lies below, above or on the plane of I 23 . Now consider the three-material cell shown in Fig. 3(b). Here, the volume occupied by M1 and M2 are polyhedrons denoted by V1 and V2, respectively. The orientation of I 23 , namely nˆ f , is calculated by minimizing an error function g (ˆn f ). Details on the construction of g (ˆn f ) will be presented later in this section. After determining the actual orientation nˆ f , the interface I 23 is located following the method described below. M1 interface I 1 is first reconstructed by matching the volume ∇ψ fraction ψ with a plane oriented as nˆ ψ = |∇ψ| (method of Youngs [1]). The plane of I 1 is then used to “clip” the polyhedron V1 , shown in red in Fig. 3(b), from the computational cell (BCEFG is clipped from ACEF i.e. V1 = Vcell − V1). The interface I 23 will, therefore, reside inside polyhedron V1 . We find the location of I 23 by translating it with the orientation nˆ f inside polyhedron V1 , until the volume V2 (BCDH in Fig. 3(b)) of M2 matches exactly with the known (target) M2 volume. For

554

A. Pathak, M. Raessi / Journal of Computational Physics 307 (2016) 550–573

Fig. 3. Possible configurations in a three-material cell: (a) The interface I 1 lies completely “below” the plane of I 23 . (b) Interfaces I 1 and I 23 intersecting. (c) Interface I 1 lies “above” the plane of I 23 .

Fig. 4. I 23 is translated along line segment v 1 v 2 inside polyhedron V1 (shown in red). 1 , 2 , . . . , 5 show various locations of I 23 during the translation. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

this purpose, we first find the greatest extent by which I 23 could be translated inside polyhedron V1 . Since V1 is a convex polyhedron, the extremities of the translation will lie on the polyhedron’s vertices. We compute a quantity dt = nˆ f · x v for each vertex of the polyhedron V1 , where x v is the vertex’s coordinate. The vertices with the minimum and maximum values of dt are labeled v 1 and v 2 having coordinates x v 1 and x v 2 , respectively. The plane I 23 can now be translated on the line segment v 1 v 2 to match the known (target) M2 volume. This has been shown in Fig. 4, where 1 , 2 , . . . , 5 show various locations of the interface I 23 as it is translated along the line segment v 1 v 2 in the polyhedron V1 shown in red. Location 1 corresponds to V2 = 0, and 5 corresponds to V2 = volume of polyhedron V1 . To pinpoint the location of I 23 on the line segment, we find the root of the following equation:

f p (x f ) = f

(5)

where, x f = μx v 1 + (1 − μ)x v 2 is a point on the line segment v 1 v 2 as μ ∈ [0, 1], and f is the target M2 volume fraction in the cell. f p is the predicted volume fraction corresponding to I 23 location x f . The root of Eq. (5) can be found using bisection, secant or Brent’s algorithm [31, Ch. 9]. We use the Brent’s algorithm, which is robust, fast and ensures convergence. Given the orientation nˆ f , locating the plane I 23 in a three-material cell, by matching exactly the M2 volume fraction, will also match the M3 volume fraction exactly. Since this scheme relies on the location of polyhedron vertices, the above explanation is valid for any other configuration than shown in Fig. 3. Next, we focus on the construction of the error function g (ˆn f ), which is minimized to find the actual orientation of I 23 , nˆ f , in a three-material cell. Consider a three-material cell c as shown in Fig. 5(a). We look in a 3 × 3 × 3 stencil around c to determine the cell with the largest value of the quantity (1 − ψ) · f · α , where, α = 1 − f − ψ is the M3 volume fraction. In the example shown in Fig. 5(a), this cell is b. For an arbitrary orientation nˆ f , I 23 is located inside b by matching the M2 volume as described in the previous paragraph. Once located, the plane I 23 is extended in the 3 × 3 × 3 stencil around cell b. In each cell N j , located in the above stencil (e.g., cell c), the interface I 1 , if it exists, is first reconstructed. Extended

A. Pathak, M. Raessi / Journal of Computational Physics 307 (2016) 550–573

555

Fig. 5. (a) Neighborhood of a three-material cell c. The interface I 23 is located by matching exactly M2 volume (V2) in cell b. I 23 is then extended in cells around cell b to construct g (ˆn f ). (b) Magnified view of cell c. Interface I 23 extended from cell b intersects with the clipped polyhedron V1 (shown in red) to obtain volume V2. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

interface I 23 is then intersected with the clipped polyhedron V1 (shown in red in Fig. 5(b)) to calculate the M2 volume, p V2, which lies “below” I 23 . We define f N as the predicted volume fraction corresponding to the volume V2 in cell N j . j

p

f N = V2/ V cell

(6)

j

where, V cell is the cell volume. Denoting the known (target) volume fraction in each cell N j by f N j , we can define a function g (ˆn f ) as,

g (ˆn f ) =

 p ( f N j − f N j )2

(7)

Nj

As mentioned before, nˆ f associated with the minimum value of g (ˆn f ) is the true I 23 orientation in the three-material cell c. Since nˆ f can be described in 3D by polar and azimuthal angles as shown in Fig. 6, Eq. (7) represents a two-parameter function. To minimize g (ˆn f ), we use the Powell’s method [31, Ch. 10], which was chosen because the derivative of g (ˆn f ) function is not readily available. In the above discussion, we have chosen cell b to locate the plane I 23 and extend it to the 3 × 3 × 3 neighborhood of cell b. Initially we chose the three-material cell c for this purpose. We had found large errors in the calculated nˆ f in the three-material cells that were almost fully occupied by M1 (i.e., ψ → 1). In such three-material cells, due to small “opening”, namely small (1 − ψ ) values, small errors in VOF advection resulted in large errors in locating the plane I 23 . This prompted us to choose a neighboring cell b with the largest value of (1 − ψ) · f · α . We now describe the construction of nˆ f as a function of polar and azimuthal angles. Fig. 6(a) shows a typical contact line formed in a three-material computational cell, where the plane I 1 with a normal vector nˆ ψ meets I 23 with a normal vector nˆ f . Denoting the contact angle between the planes I 1 and I 23 by φ , it is evident that the angle between the normal vectors nˆ ψ and nˆ f is also equal to φ . In general, nˆ f is allowed several orientations, the locus of which results in a cone whose base is a circle C shown in Fig. 6(b) lying on the plane I 1 . Although φ is shown by an angle in Fig. 6, the  acute 

→ −

current discussion applies also to obtuse contact angles. Since, nˆ f is a unit vector, the radius AB of circle C is sin φ and

− →

OA = cos φ nˆ ψ . We seek two perpendicular unit vectors on the plane I 1 : aˆ and bˆ shown in Fig. 6(b). aˆ can be any vector on the plane I 1 and can easily be extracted from the 3D reconstruction scheme [1], which is used to locate the plane I 1 . Then

− → a perpendicular vector bˆ can be calculated as aˆ × nˆ ψ . AB can now be represented in terms of aˆ and bˆ as:

− →









→ → − − AB = AB cos θ aˆ + AB sin θ bˆ

(8)

556

A. Pathak, M. Raessi / Journal of Computational Physics 307 (2016) 550–573

Fig. 6. (a) Contact angle φ shown between M1 interface I 1 and plane I 23 . (b) For contact angle φ between nˆ f and nˆ ψ , the locus of nˆ f is a cone whose base is a circle C lying on the plane I 1 .

−  →

where AB = sin φ . Now, we can write:

− →

nˆ f = OB

− → − → = OA + AB = cos φ nˆ ψ + sin φ cos θ aˆ + sin φ sin θ bˆ

(9)

Representing nˆ f as a function of φ and θ , as shown in Eq. (9), makes g (ˆn f ) a two-parameter function for a given nˆ ψ . A typical g (ˆn f ) surface profile, taken from the test case that will be presented in Section 3.2.1, is plotted against θ and φ in Fig. 7(a). For contact line applications, the contact angle φ is either computed by empirical relations [32,33], or available from experimental data [34]. Using the additional information provided by the contact angle φ , reduces nˆ f , and hence g (ˆn f ), from a two-parameter variable to a single-parameter variable. An example for the single-parameter case is shown in Fig. 7(b), where the imposed contact angle is 90◦ , i.e., φ = π /2, and g (ˆn f ) is plotted as a function of θ . The curve shown in Fig. 7(b) is a slice made at φ = π /2 of the surface profile shown in Fig. 7(a). In all test cases presented in this paper, we have considered the most general two-parameter formulation of the function g (ˆn f ). The minimization has been performed by employing the multi-dimensional Powell’s method [31, Ch. 10]. If one wishes to use the one-parameter formulation of g (ˆn f ), the minimization can be performed using one-dimensional routines such as the golden section search method or Brent’s method [31, Ch. 10]. In Fig. 7(a), we can see a local minimum near the crest. If the initial guess is not appropriate, the minima calculators are susceptible to capturing this local minimum instead of the global one. To address this issue, we keep the initial guess provided to Powell’s method sufficiently close to the global minimum. To do this, we compute g (ˆn f ) at [N × N] regularly spaced points in (φ, θ ) space, where φ ∈ [0, π ] and θ ∈ [0, 2π ]. We have found N = 10 to be sufficient. The point corresponding to the minimum g (ˆn f ) among [ N × N ] points is provided as the initial guess to the Powell’s method. Similar approach can be applied to the one-parameter formulation of g (ˆn f ) to obtain the initial bracketing needed by the 1D minimization routines. 2.2. Additional contact line configurations In a 3D computational domain, the contact line (intersection of interfaces I 23 and I 1 ) would generally lie in threematerial cells. Resolving the interfaces in three-material cells would, therefore, be tantamount to capturing the contact line (if it exists). Although rare, there can be special configurations when a contact line does not fall exactly inside three-material cells, but on cell faces instead. Fig. 8 shows all configurations (including three-material cells) in which a contact line could exist. These additional configurations need to be identified and resolved in order to capture the contact line. As shown in Fig. 8(b), if the contact line lies on the face of a two-material cell c (containing M2 and M3), the I 23 orientation nˆ f can be calculated by constructing and minimizing the function g (ˆn f ) using the steps described in the previous section. Information aˆ and bˆ on M1 interface I 1 can be retrieved from one of the neighboring cells containing M1 (cell b in Fig. 8(b)). For the cases shown in Figs. 8(c) and 8(d), the interface I 23 is coincident with one of the cell faces, so there is no need to construct a minimization function there. 2.3. Extending I 23 interface into Material 1 In two material cells, using the color function f to compute the normal vector to the interface I 23 , for example via nˆ f =

∇f |∇ f | [2, p. 37], requires sufficient stencil, which may not be available when the f -field is truncated by the third material, i.e.,

A. Pathak, M. Raessi / Journal of Computational Physics 307 (2016) 550–573

557

Fig. 7. (a) Surface plot of g (ˆn f ) against angles φ and θ . (b) g (ˆn f ) against θ for a constant φ = π /2. In both plots, a local minimum is visible near the crest. The example is taken from the test case presented in Section 3.2.1.

material M1 represented by ψ . In such cases, the f -field is discontinuous and ∇ f calculation may be, therefore, inaccurate. As suggested by Ghasemi et al. [4], a “virtual” f -field can be constructed by extending the I 23 interface, and effectively the f -field, to material M1 as described in the example shown in Fig. 9. The continuous virtual f -field is then used to compute the orientation of the interface I 23 . We adopt the same approach in the present work and extend it to 3D. After locating the contact line in a three-material cell, or in one of the special configurations as mentioned in the last section, the plane I 23 in that three-material cell is ∇f extended into M1 in a 5 × 5 × 5 stencil. This stencil is large enough for calculating the orientation nˆ f = |∇ f | in two-material cells containing M2 and M3, which are located in the vicinity of the contact line, by following the procedure explained next. The extended M2 volume fraction is stored in a variable called f virt. , which is different from the known (actual) volume fraction f and is shown in Fig. 9. The f virt. field can now be used to compute the orientation of I 23 in the two-material ∇ f virt.

cells via nˆ f =  virt.  . Present algorithm thus not only calculates the interface orientation in three-material cells, but also ∇f ensures reasonably accurate calculation of the interface orientation in those two-material cells that lie in the vicinity of the contact line. 2.4. Advection in three-phase cells Following the reconstruction of the interfaces I 1 and I 23 , we can choose from a variety of advection schemes. We adopt the operator split advection scheme [1,26,35] to transport the material volumes. In two-material cells, we use the analytical formulae provided by Youngs [1] to calculate the flux volumes. The advection in three-material cells will require a discussion here. Consider a three-material cell as shown in Fig. 10 where the planes I 1 and I 23 have been already located.

558

A. Pathak, M. Raessi / Journal of Computational Physics 307 (2016) 550–573

Fig. 8. Contact line can lie in one of the four configurations: (a) in a three material cell, (b) on the cell face of a two-material (M2–M3) cell, (c) I 23 coincides with a cell face of an M2 filled cell, (d) I 23 coincides with a cell face of an M1–M2 filled cell.

Fig. 9. Interface I 23 extended into M1 in a 5 × 5 × 5 neighborhood around the contact line T. Extended interface has been shown by the dashed line. The M1 volume lying below the dashed line plus the M2 volume (shown in blue) make f virt. . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Considering advection in the positive x-direction, the total volume advected out (flux volume) would reside in the rectangle (cuboid in 3D) BCIH. With the geometric constructs presented in the Appendices and the Supporting Information document, it is possible now to compute the flux volumes of M1, M2 and M3 represented by polyhedrons BCDF, FEG and GEDIH, respectively. Material volumes advected in time-step t are then used to solve geometrically the scalar transport equations of f and ψ ,

∂f + u · ∇ f = 0 ∂t ∂ψ + u · ∇ψ = 0 ∂t

(10) (11)

Eqs. (10) and (11) would implicitly satisfy, on a discrete level, the transport equation for M3, the volume fraction (color function) of which is given by α = 1 − f − ψ .

A. Pathak, M. Raessi / Journal of Computational Physics 307 (2016) 550–573

559

Fig. 10. Flux volumes due to advection by the velocity component u through the right face of a three-material cell. The advected volume (volume flux) of M1, M2 and M3 are represented by polyhedrons BCDF, FEG and GEDIH, respectively.

3. Results In this section, we will assess the performance of the new method in resolving three-material configurations. We start with static test cases which report the accuracy with which the new methodology captures the interface orientations in three-material cells. In the static test cases, the exact volume fractions are initialized. Then, the new methodology is used to compute the interface normal vectors. The results are compared with the exact interface orientation for error assessment. Advection test cases are then presented which demonstrate the performance of the new methodology when working in conjunction with an advection scheme [1,26,35]. We are not imposing a contact angle in any of the test cases presented in this section. The error function g (ˆn f ) is formulated and minimized as a two-parameter function, which is the most general case. In the figures presented here, yellow–brown is used to represent material M1, while M2 is represented by blue–green. M3 is left transparent for the sake of contact line’s visibility. 3.1. Static problems We define here three error norms, l∞ , l1 and l2 to evaluate the accuracy of the overall algorithm in estimating the ˆ Note that nˆ f = ∇ f (see Ref. [2, p. 37] for details) is used to calculate I 23 orienI 23 orientation nˆ f = nx ˆi + n y ˆj + n z k. |∇ f | tation in two-material cells while the new reconstruction algorithm, presented in section 2, estimates I 23 orientation in three-material cells.

 

 

 

 

 

 

. . . e j = ncalc − nexact − nexact − nexact  + ncalc  + ncalc  x y z x y z

l∞ = max e j

(13)

j

l1 =

1  N

(12)

ej

(14)

j

   1 e 2j l2 =  N

(15)

j

Here, e j in Eq. (12) adds absolute differences between calculated and exact normal vector components in each direction. The error norms’ calculations are performed over both two-material and three-material computational cells containing the interface I 23 . N represents the number of such cells belonging to either category. 3.1.1. Reconstruction of a planar interface intersecting a curved surface Problem setup is a planar interface I 23 intersecting a hemisphere of material M1 at 90◦ contact angle as shown in Fig. 11. The computational domain size is 0.625 × 0.625 × 0.3125. The radius of the hemisphere is 5/32 and 20 cells have been used to resolve it. The hemisphere’s center is at (0.3125, 0.3125, 0.0). The exact volume fractions corresponding to the planar interface I 23 , shown in Fig. 11,√are initialized. The planar interface passes through the center of the hemisphere and √ has an exact orientation of nˆ f = (1/ 2, 0, 1/ 2), which will be compared against the orientations computed by the new algorithm here. Table 1 shows the accuracy of the overall algorithm in computing the plane I 23 orientation nˆ f . The new methodology ∇f

calculates nˆ f exactly in three-material cells. Computing nˆ f in two-material cells relies on the ∇ f calculation (nˆ f = |∇ f | ),

560

A. Pathak, M. Raessi / Journal of Computational Physics 307 (2016) 550–573

Fig. 11. A planar interface I 23 intersecting a hemisphere of material M1 at a contact angle of 90◦ . Material M2 is shown in transparent blue for visibility of the hemisphere. Material M3 is left completely transparent.

Table 1 Error norms associated with nˆ f calculation by (1) the new methodology in three-material cells, and by (2) ∇ f method in two-material cells. 1st column reports the number of cells belonging to each category. # Cells

l∞

l1

l2

2-material cells 2722

1.675 × 10−14

3.875 × 10−16

7.065 × 10−16

3-material cells 66

3.345 × 10−14

5.181 × 10−15

9.237 × 10−15

Fig. 12. Spherical cap of material M2 sitting on an M1 sphere. The interfaces I 23 and I 1 meet at a contact angle of 120◦ .

and is accurate to the machine precision for the present I 23 orientation.1 Furthermore, in the two-material cells near the contact line, the f virt. field, extended into M1, was used in the nˆ f calculation. Having obtained exact nˆ f in all two-material cells, as shown in Table 1, indicates that the f virt. extension was also done accurately. First column in Table 1 provides the number of the two-material and three-material cells where the interface I 23 lies. Because the number of two-material cells is significantly more than that of three-material cells, using ∇ f method to calculate nˆ f in two-material cells saves computational time appreciably. 3.1.2. Reconstruction of a non-planar interface intersecting a curved surface We now demonstrate the performance of the proposed algorithm in reconstructing a non-planar interface I 23 . Fig. 12 shows a spherical cap made of material M2 sitting on top of a sphere, which is made of M1. The contact line where the interfaces I 23 and I 1 meet is a circle. The two spheres have the same radius R and the distance between their centers is d. From the mid-plane slice of Fig. 12, the contact angle φ between the interfaces I 23 and I 1 is then:

1

It should be noted that ∇ f returns nˆ f with machine precision in this particular interface orientation because the discrete

∂f ∂f and computed here ∂x ∂z ∇f

have the same magnitude due to the 45◦ orientation of the plane and its symmetry with respect to the xz plane. For a general interface orientation, will not compute “exact” normals [13].

A. Pathak, M. Raessi / Journal of Computational Physics 307 (2016) 550–573

561

Table 2 Error in nˆ f calculated by ∇ f method in two-material cells and the new methodology in three-material cells.

x

# Cells

2-material cells 1/64

l∞ 2.506 × 10−1

1636

1/128

6660

1.305 × 10−1

1/256

26 992

1.364 × 10−1

1/512

108 488

1.387 × 10−1

64

4.089 × 10−1

3-material cells 1/64

Order

1/128

116

1.763 × 10−1

1/256

228

9.603 × 10−2

1/512

436

4.338 × 10−2

l1

Order

4.027 × 10−2 0.94

−0.064 −0.024

3.742 × 10−2 3.734 × 10−2 3.712 × 10−2

0.88 1.15

1.391 × 10−1 8.086 × 10−2 3.315 × 10−2

Order

5.258 × 10−2 0.11 0.003 0.009

3.181 × 10−1 1.21

l2

4.768 × 10−2 4.805 × 10−2 4.818 × 10−2

0.14

−0.011 −0.004

3.250 × 10−1 1.19 0.78 1.29

1.410 × 10−1 8.185 × 10−2 3.366 × 10−2

1.20 0.78 1.28

Fig. 13. (a) Spherical cap made of material M2 sits on an M1 sphere. (b) Interface I 23 is extended inside M1 in a narrow band around the contact line. The already present M2 (in blue) plus the extended material (in purple) make the f virt. field. For the images shown, 80 cells were used to resolve the diameter. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

φ = 2 sin−1 ( √

d 2R

)

(16)

We choose d = 3R, R = 5/32, which yields a contact angle φ = 120◦ . Table 2 shows the accuracy with which I 23 orientation nˆ f is calculated by ∇ f method [2] in two-material cells and by the new methodology in three-material cells at different resolutions. The second column reports the number of cells belonging to either category. The ∇ f method used to compute nˆ f is known to be zeroth-order accurate [36,37]. Same order of convergence can be seen in Table 2 for nˆ f calculated in two-material cells. Table 2 also shows that the new methodology yields the interface I 23 orientation with 1st order accuracy in three-material cells. Fig. 13(b) shows the extension of I 23 inside M1 in a narrow band around the contact line. The extension was performed using the scheme described in Section 2.3. The purple patch plus the M2 volume (shown in blue) make the f virt. field, which is used to calculate nˆ f in two-material cells via ∇ f virt. . The location of l∞ in two-material cells (see Table 2) is well away from the contact line at all resolutions except at x = 1/64. It means that the l∞ error is due to the error in ∇ f calculation and not due to f virt. extension. At the resolution x = 1/64, the location of l∞ error in two-material cells is in the vicinity of the contact line because large errors are associated with nˆ f calculation in three-material cells and hence in f virt. extension. 3.2. Advection problems We now assess the performance of the three-material reconstruction algorithm when used in concert with an advection scheme. For the results presented here we have used operator split volume-of-fluid advection scheme [1,26,35]. To assess the accuracy of advection, we estimate the error in the volume fraction distribution of material M1, represented by ψ , and M2, represented by f , at the end of advection. The formula for L1 error “E c ” has been borrowed from [3] and is given by:

562

A. Pathak, M. Raessi / Journal of Computational Physics 307 (2016) 550–573

Fig. 14. Advection of T-shaped interfaces formed by three materials.

    n   n   f i , j ,k − f ie, j ,k  + αi , j ,k − αie, j ,k 

Ec =

i , j ,k

i , j ,k

i , j ,k

f ie, j ,k

+



i , j ,k

αie, j,k

(17)

where, α = 1 − f − ψ . Superscript n denotes the volume-fraction at the end of the simulation and superscript e denotes the corresponding exact volume fraction. Also reported in this section are errors in the normal vector nˆ f at the end of advection. The definitions of nˆ f error norms are the same as given in Eqs. (13)–(15). In general, the volume fraction errors at the end of advection are expected to converge with first-order accuracy similar to the method of Youngs as reported in [13,36,37]. The normal vector nˆ f computed from ∇ f , therefore, is expected to show zeroth-order convergence [36,37]. 3.2.1. Advection of T-shaped interfaces We begin with a test case where the interfaces I 23 and I 1 are both planar and form a T-shape as shown in Fig. 14(a). The computational domain size is 0.625 × 0.625 × 0.625, which is resolved using a uniform 80 × 80 × 80 grid. M1 interface I 1 has a normal vector nˆ ψ = (0, 0, −1) and nˆ f of I 23 is (−1, 0, 0), corresponding to 90◦ contact angle. The contact line at t = 0 passes through the point (20.5 x, 0, 20.5 x), where x is the uniform finite volume cell size in each direction. A constant velocity of (1/2, 0, 1) is prescribed in the entire computational domain. Time step is set as t = x/10. The simulation was performed for a total time t f = 400 t. Fig. 14 shows snapshots of the T-shaped interfaces at different times. Table 3 reports error E c in volume fractions at the end of the simulation, along with the numerator (N r ) and denominator (D r ) in the expression for E c in Eq. (17). V f α is the total exact volume of M2 and M3 at the end of the simulation, which is known analytically and could also be computed by multiplying D r by cell volume x3 . Table 4 shows errors in normal vector nˆ f at the end of advection. Since the planar interfaces are captured exactly by the overall algorithm, a grid convergence study for this test case was not performed. Note that since the I 23 orientation remains (−1, 0, 0) during the simulation, nˆ f is computed exactly using ∇ f method inside two-material cells. 3.2.2. Advection of Y -shaped interfaces In this test case, M1 interface I 1 meets the interface I 23 at a contact angle of 45◦ as shown in Fig. 15(a). The computational domain size, cell size, time-step, prescribed velocity and total simulation time t f are the same as the previous test case (Section 3.2.1). Fig. 15 shows snapshots of the simulation at different times. Tables 3 and 4 report errors in volume

A. Pathak, M. Raessi / Journal of Computational Physics 307 (2016) 550–573

563

Table 3 VOF advection errors in the advection of T- and Y -shaped interfaces. Simulation

Nr

Dr

V f α (= D r x3 )

Ec

T-advection Y -advection

2.364 × 10−9 2.671 × 10−9

1.248 × 105 2.592 × 105

5.951 × 10−2 1.236 × 10−1

1.895 × 10−14 1.030 × 10−14

Table 4 nˆ f error norms in the advection of T- and Y -shaped interfaces. Simulation

l∞

l1

l2

2-material cells T-advection Y -advection

6.968 × 10−17 7.432 × 10−14

1.271 × 10−18 8.504 × 10−15

5.624 × 10−18 1.867 × 10−14

3-material cells T-advection Y -advection

2.833 × 10−16 3.900 × 10−14

2.833 × 10−16 3.724 × 10−14

2.833 × 10−16 3.725 × 10−14

Fig. 15. Advection of Y -shaped interfaces formed by three materials.

fractions and in nˆ f , respectively, at the end of the advection. Again, because the planar interfaces are captured exactly by the overall algorithm, a grid convergence study for this test case was not conducted. Note that since the I 23 orientation remains 45◦ in the x–z plane during the simulation, nˆ f is computed exactly by ∇ f method inside two-material cells. It is worth mentioning here that for both T- and Y -shaped interface advections, the time-step t was chosen such that the contact line would lie on a cell-face many times during the simulation. There would then be a situation when the contact line would exist even in the absence of three-phase cells. The algorithm formulation (as mentioned in Section 2.2) is capable of resolving contact line in such scenarios. It is, therefore, able to preserve the accuracy with which it captures planar interfaces to the end of the simulation. 3.2.3. Advection of a curved contact line Initial problem setup for this test case is same as the static test case presented in Section 3.1.1. A rotational velocity, ω = (0, 1, 0), is prescribed about an axis that passes through the center of the M1 hemisphere located at (0.3125, 0.3125, 0) .

564

A. Pathak, M. Raessi / Journal of Computational Physics 307 (2016) 550–573

Fig. 16. Advection of a curved contact line formed by three materials. Angular velocity of the planar interface is (0, 1, 0).

Table 5 VOF errors in the advection of a curved contact line shown in Fig. 16. Last column presents the order of convergence of error E c defined in Eq. (17).

x

Nr

Dr

V f α (= D r x3 )

Ec

1/64

1.225 × 102

2.991 × 104

1.141 × 10−1

4.096 × 10−3

1/128

4.632 × 102

2.392 × 105

1.141 × 10−1

1.936 × 10−3

1/256

1.807 × 103

1.914 × 106

1.141 × 10−1

9.441 × 10−4

1/512

7.073 × 103

1.531 × 107

1.141 × 10−1

4.620 × 10−4

Order 1.08 1.04 1.03

Table 6 nˆ f error norms in the advection of a curved contact line shown in Fig. 16. No convergence is shown by the error norms.

x

l∞

2-material cells 1/64

7.103 × 10−2

1/128

6.467 × 10−2

1/256

7.875 × 10−2

1/512

9.199 × 10−2

3-material cells 1/64

6.824 × 10−2

1/128

6.351 × 10−2

1/256

7.984 × 10−2

1/512

1.044 × 10−1

Order

l1

Order

4.500 × 10−3 0.14

−0.28 −0.22

2.619 × 10−3 2.049 × 10−3 1.478 × 10−3

−0.33 −0.39

2.759 × 10−2 2.377 × 10−2 3.186 × 10−2

Order

1.115 × 10−2 0.78 0.35 0.47

3.057 × 10−2 0.10

l2

6.533 × 10−3 5.704 × 10−3 4.823 × 10−3

0.77 0.20 0.24

3.676 × 10−2 0.15 0.22

−0.42

3.259 × 10−2 3.050 × 10−2 4.070 × 10−2

0.17 0.10

−0.42

The simulation was run for a total time t f = π /2 with a time-step t = π x/(25R ), where R = 5/32 is the radius of the M1 hemisphere. Fig. 16 shows snapshots of the simulation at different times. Table 5 reports the error in volume fractions E c at the end of the simulation at different grid resolutions. In this test case, E c shows first-order convergence similar to the method of Youngs [1] because the plane I 23 did not keep the same orientation during the simulation as in the last two test cases. Table 6 reports errors in nˆ f at the end of the simulation. As expected, no convergence is shown by the error norms either in two-material cells or in three-material cells in this case. This is consistent with the previously reported results in the literature on the convergence of nˆ f computed by using ∇ f [36,37]. Note that the order of accuracy of nˆ f in the three-material cells reduces to that of two-material cells because the calculation relies on the advected f in two-material cells.

A. Pathak, M. Raessi / Journal of Computational Physics 307 (2016) 550–573

565

Fig. 17. Snapshots of rotation of a two-material sphere in the presence of a third material M1. M2 and M3 are shown in blue and transparent purple respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Table 7 VOF errors in the rotation of a two-material sphere presented in Fig. 17. First-order convergence is shown by the error E c . V f α (= D r x3 )

Ec

3

1.598 × 10−2

2.951 × 10−2

4.098 × 102

3.352 × 104

1.598 × 10−2

1.223 × 10−2

1/256

1.818 × 103

2.681 × 105

1.598 × 10−2

6.781 × 10−3

1/512

8.059 × 103

2.145 × 106

1.598 × 10−2

3.757 × 10−3

x

Nr

1/64

1.236 × 10

1/128

Dr 2

4.189 × 10

Order 1.27 0.85 0.85

3.2.4. Rotation of a two-material sphere Next, we present a 3D extension of the test case used by Choi and Bussmann [3], where a two-material sphere performs a full counter-clockwise rotation inside another material. We initialize a sphere of radius 5/32. As shown in Fig. 17(a), the hemisphere in blue is made of M2 and that in transparent purple is made of M3. The sphere is surrounded by M1, which has been made completely transparent in Fig. 17 for the sphere’s visibility sake. Initially, the interface I 23 is a plane passing through the center of the sphere with the normal vector orientation nˆ f set as (cos θ1 cos θ2 , sin θ1 , cos θ1 sin θ2 ), where θ1 = π /6 and θ2 = −π /3. At the contact line, the interfaces I 1 and I 23 meet at a contact angle of 90◦ . The domain size is  = (0, −1, 0) is specified in the 1 × 0.5 × 1 and the sphere’s center is initially located at (0.75, 0.25, 0.5). Angular velocity ω computational domain about an axis passing through the center of the domain, (0.5, 0.25, 0.5). Time step t is 32π x/325 and the simulation was run for a total time t f = 2π . Fig. 17 shows snapshots of the two-material sphere at different times. Table 7 reports the volume fraction error E c at the end of the simulation at different mesh resolutions. First order convergence is observed in the errors. V f α is the total volume of M2 and M3 and, in the present case, is equal of the volume of the sphere. Table 8 reports the errors in nˆ f at the end of the simulation. Again, as expected, no convergence is shown by the error norms. The negative convergence seen in some cases here is due to the VOF advection errors, which are more pronounced in this test compared to others; see the volume fraction errors, E c , in Table 7. The VOF advection errors can emanate from two sources: (1) inaccuracy associated with the calculation of the interface normal, which negatively impacts the interface reconstruction and volume flux calculations, and (2) directional VOF advection schemes. We will test the former in Section 3.3, where the exact interface normal vectors are used in the two-material cells and the impact on the three-material interface reconstruction and on E c is examined. Testing the latter source is beyond the scope of this work, but we believe multidimensional advection schemes, such as that proposed in [38], can reduce the volume fraction errors, which in turn, further improves the accuracy of the three-material interface reconstruction. 3.2.5. Advection of a spherical cap along the surface of a sphere Finally, we consider a test case that involves the advection of a non-planar I 23 interface. Similar to the setup presented in Section 3.1.2, a spherical cap made of M2 resides on an M1 sphere. The radii of the two spheres are the same and

566

A. Pathak, M. Raessi / Journal of Computational Physics 307 (2016) 550–573

Table 8 nˆ f errors in the rotation of a two-material sphere presented in Fig. 17. No convergence shown by the error norms.

x

l∞

2-material cells 1/64

2.827 × 10−1

1/128

2.950 × 10−1

1/256

3.538 × 10−1

1/512

5.277 × 10−1

3-material cells 1/64

4.000 × 10−1

1/128

3.441 × 10−1

1/256

5.320 × 10−1

1/512

5.748 × 10−1

Order

−0.061 −0.26 −0.58

l1

Order

5.728 × 10−2 4.216 × 10−2 4.343 × 10−2 4.574 × 10−2

−0.63 −0.11

7.178 × 10−2 9.823 × 10−2 1.544 × 10−1

Order

7.100 × 10−2 0.44

−0.043 −0.074

1.016 × 10−1 0.22

l2

5.403 × 10−2 5.625 × 10−2 6.612 × 10−2

0.39

−0.058 −0.23

1.218 × 10−1 0.50

−0.45 −0.65

1.043 × 10−1 1.362 × 10−1 2.044 × 10−1

0.22

−0.38 −0.59

Fig. 18. Advection of a spherical cap along the surface of a sphere in the presence of a third material. The circular contact line traverses 120◦ .

√

are equal to R = 5/32. The distance between the centers of the spheres √ is 3R. The center of the M1 sphere is at (1/2, 15/64, 15/64) and that of M2 spherical cap is (17/64, 15/64, (15 + 5 3)/64). The line segment joining the two centers is perpendicular to the y-axis and makes a 30◦ angle with the negative x-axis at t = 0; see Fig. 18(a). The domain size is  = (0, 1, 0) is prescribed about M1 sphere’s center. The time step t is 1.0 × 0.46875 × 0.78125 and a rotational velocity ω set as 2π x/(75R ) and the simulation was run for t f = 2π /3. The contact line in this case is a circle which traverses an angle of 2π /3 on M1 sphere’s surface. Fig. 18 shows snapshots of the simulated results at different times. Table 9 reports the VOF error E c at the end of the simulation at different mesh resolutions. First order convergence is seen in E c errors. We have also presented in Table 10 nˆ f errors at the end of the simulation. No definite convergence is seen in the nˆ f errors.

A. Pathak, M. Raessi / Journal of Computational Physics 307 (2016) 550–573

567

Table 9 VOF errors in the advection of a spherical cap along the surface of another sphere, presented in Fig. 18.

x

Nr

Dr

V fα

Ec

1/64

1.134 × 102

9.181 × 104

3.502 × 10−1

1.235 × 10−3

1/128

4.501 × 102

7.345 × 105

3.502 × 10−1

6.128 × 10−4

1/256

1.889 × 103

5.876 × 106

3.502 × 10−1

3.215 × 10−4

1/512

7.664 × 103

4.701 × 107

3.502 × 10−1

1.630 × 10−4

Order 1.01 0.93 0.98

Table 10 nˆ f errors in the advection of a spherical cap along the surface of another sphere, presented in Fig. 18.

x

l∞

2-material cells 1/64

3.408 × 10−1

1/128

2.157 × 10−1

1/256

1.911 × 10−1

1/512

2.273 × 10−1

3-material cells 1/64

5.123 × 10−1

1/128

3.262 × 10−1

1/256

2.554 × 10−1

1/512

2.296 × 10−1

Order

l1

Order

l2

6.009 × 10−2 0.65 0.17

−0.25

7.510 × 10−2 0.23

5.138 × 10−2

6.209 × 10−2

0.16

4.607 × 10−2

5.710 × 10−2

0.068

4.394 × 10−2

5.522 × 10−2

3.627 × 10−1 0.65 0.35 0.15

Order

0.27 0.12 0.048

3.721 × 10−1 0.93

1.905 × 10−1

1.991 × 10−1

0.65

1.214 × 10−1

1.325 × 10−1

0.48

8.725 × 10−2

1.018 × 10−1

0.90 0.59 0.38

Table 11 Volume fraction error, E c , and error norms of nˆ f in three-material cells in a special case of the test presented in Section 3.2.4 where the orientation of I 1 as well as the orientation of I 23 in the two-material cells are exactly prescribed.

x

Ec

1/128

6.49 × 10−3

Order

1/256

3.57 × 10−3

l∞

Order

2.02 × 10−1

(0.86)

2.35 × 10−1

l1

Order

3.34 × 10−2

(−0.22)

4.69 × 10−2

l2

Order

4.62 × 10−2

(−0.49)

6.54 × 10−2

(−0.50)

3.3. Relation between interface reconstruction accuracy in two- and three-material cells ∇ψ In the test cases presented thus far, we have reconstructed the interface I 1 using nˆ ψ as |∇ψ| and I 23 in two-material ∇f

cells using nˆ f as |∇ f | following an approach suggested by [2]. If more accurate methods of calculating nˆ ψ and nˆ f , such as height functions [39–41], are employed, the accuracy of the interface reconstruction in three-material cells can be improved too. To test this argument, we repeated the test case presented in Section 3.2.4 at resolutions x = 1/128 and x = 1/256. However, this time, the I 1 orientation nˆ ψ and the I 23 orientation nˆ f in two-material cells were prescribed using the exact values throughout the simulation. The interface reconstruction in three-material cells was still carried out using the error minimization method presented in this paper. At the end of simulations, the error E c associated with volume fractions and the errors associated with nˆ f in three-material cells were evaluated and are reported in Table 11. Comparing the errors presented here with those reported in Tables 7 and 8, it is observed that the magnitude of the errors was halved when the exact normal vectors were used. That observation indicates that using more accurate methods for advecting two-material interfaces would also benefit the accuracy of the new method in reconstructing and advecting three-material interfaces. 3.4. Computational time of three-material interface reconstruction Although the proposed three-material interface reconstruction may seem computationally expensive, it does not become a “computational bottleneck” when the simulation resolution is increased. That is because the ratio of the three-material to two-material cells decreases as the resolution is increased. To demonstrate this, we report in Table 12 the computational

568

A. Pathak, M. Raessi / Journal of Computational Physics 307 (2016) 550–573

Table 12 Computational time associated with VOF advection in the test case presented in Section 3.2.5. At each resolution, the total time for VOF advection is presented in the 2nd column, and the time taken by the three-material interface reconstruction is shown in the 3rd column. The 4th column gives the percentage of the total time that is taken by the three-material interface reconstruction. The number of three- and two-material cells are reported in the 5th and 6th columns, respectively.

x

Total time (sec.)

3-mat. recons. time (sec.)

%

No. of 3-mat. cells

No. of 2-mat. cells

1/64 1/128 1/256 1/512

3.43 9.72 49.14 517.0

3.16 6.69 12.71 27.40

92 69 26 5.3

70 126 254 502

1676 6806 27 496 110 408

time associated with the VOF advection for the test case presented in Section 3.2.5. Table 12 presents the total time2 for the VOF advection (2nd column) and the time spent on the three-material interface reconstruction (3rd column) for one time cycle at each resolution. It is seen from the 4th column that while the three-material interface reconstruction takes 92% of the total computational time at the lowest resolution x = 1/64, it takes only 5.3% of the total time as the resolution is increased to x = 1/512. That is because as the resolution increases, the number of three-material cells (5th column) increases at a smaller rate compared to the two-material cells (6th column). The table shows that the number of threematerial cells is almost doubled with each increase in the resolution, while the number of two-material cells is quadrupled. 4. Conclusion We have introduced a novel VOF-based method to perform 3D interface reconstruction and advection in three-material cells. The method relies on a two-parameter minimization of an error function. The minimization problem formulation also allows the user to impose a contact angle at the contact line formed by three materials. If a contact angle is imposed, the minimization reduces to a single-parameter problem. In addition to three-material cells, the above methodology can be used in two-material cells for interface reconstruction; however, that would increase the computational cost significantly. Therefore, in two-material cells, it is more efficient to use the gradient of volume fraction, ∇ f , to calculate the interface orientation and then employ the method of Youngs [1] to reconstruct the interface. The proposed three-material reconstruction is general and can also be used in an unstructured tetrahedral grid in conjunction with the analytical two-material reconstruction developed by Lv et al. [42]. Although the volume advection presented in this paper was performed using the operator split scheme [26], the proposed reconstruction is also suitable for the unsplit advection scheme [26]. Performance of the new algorithm was assessed via both static and advection test cases. The algorithm is able to reconstruct planes exactly and capture the orientation of non-planar interfaces with first-order accuracy. For advection test cases, the accuracy of the new methodology in reconstructing three-material interfaces is the same as the accuracy of the method used for reconstructing two-material interfaces, which was the Youngs’ method [1] here. If the accuracy of the twomaterial interface reconstruction is increased by, e.g., using more accurate methods for calculating interface orientation such as height function [39–41], then the accuracy of the three-material interface reconstruction and advection is also improved. A complete set of geometric constructs are presented in the Appendices. They are based on convex 3D polyhedra making the proposed reconstruction suitable for both structured and unstructured grids. Acknowledgements The research support from the National Science Foundation under the CBET Grant Nos. 1236462 and 1336232 is gratefully acknowledged. The simulations presented here were performed on the UMass Dartmouth HPC cluster, which was purchased through partial support from the NSF (MRI Grant No. CNS-0959382) and AFOSR (DURIP Grant No. FA9550-10-1-0354). AP gratefully acknowledges the Fellowship from the Office of the Associate Provost for Graduate Studies at UMass Dartmouth. Appendix A Here we present the geometric functions necessary to perform the new three-material reconstruction algorithm. The definitions of the data structures and the pseudo-codes are available in the Supporting Information document. A.1. Data structures Data structures are important for conciseness and enhanced comprehensiveness of the geometric functions. We define four data structures: (1) st_point, (2) st_face, (3) st_polyhedron, and (4) Quat. Structure st_point is used to represent both Cartesian coordinates of a 3D point and a vector. Structures st_polyhedron and st_face store the

2 For measuring time, the simulations were performed on a single processor of a compute node with two Intel Xeon (quad-core) E5620 2.4 GHz processors, 24 GB 1596 MHz RAM.

A. Pathak, M. Raessi / Journal of Computational Physics 307 (2016) 550–573

569

Fig. 19. (a) A cubic polyhedron representing a computational cell. (b) Cubic polyhedron when cut by a plane shown in red. The volume lying below the plane is shown in blue. Point xc and vector aˆ (in IJ direction) lie on the plane of the polygon IJKLMN.

information pertaining to a polyhedron and a polyhedron face, respectively. Quat represents a quaternion, which will be used to rotate a 3D polygon. A.2. 3D geometric functions We now describe the functions constituting the 3D geometrical toolbox. We will follow an order of precedence while describing the functions: the functions that are called are described before the calling functions. This order would make it easier for the reader to comprehend expeditiously the overall implementation. Since only the volume-fraction information is needed, all calculations are performed on a cubic cell’s local coordinate system. In the local coordinate system, cell span is [0, 1] × [0, 1] × [0, 1] instead of actual cube vertices’ coordinates. The functions are designed in a way to keep the computational costs to a minimum. Interface reconstruction in two-material cells Template: void area_3D (double v m , st_point nˆ m , st_point* xc , st_point* aˆ ). The function accepts as inputs volume fraction v m of a material m and a unit normal vector nˆ m . nˆ m is the orientation of the plane representing material m interface in a given computational cell. In the current context, v m can be f or ψ and nˆ m can be nˆ f or nˆ ψ . The function performs geometrical reconstruction of a 3D polygon resulting from the intersection of the material m interface plane and the computational cell, and yields as outputs: (1) centroid of the 3D polygon, xc , and (2) a unit vector aˆ lying in plane of the polygon. Here, we used the 3D reconstruction scheme of Youngs [1], designed for two-phase flows, and modified it to extract xc and aˆ . For succinctness, we do not provide a pseudo-code for this particular function. The reader is referred to [1] for a detailed description of the 3D two-material reconstruction scheme. Fig. 19(a) represents a cubic cell, which is cut in Fig. 19(b) by a plane, shown in red, having an orientation nˆ m . The material m lying below the plane is shown in blue, and has a volume v m . In the context of unstructured tetrahedral grids, the reader is referred to the analytical two-material reconstruction scheme developed by Lv et al. [42]. Constructing a simple polyhedron Template: void const_sim_polyhedron (st_polyhedron* cube_poly). The function constructs the polyhedron associated with a unit cubic cell. Individual faces (st_face) are constructed by storing the vertices in the counter-clockwise (CCW) direction. Detailed construction is provided in the pseudo-code available in the Supporting Information document. Constructing face of a polyhedron Template: void const_face (st_point* cell_vert, int n_verts, st_point* xc , st_point* nˆ m , st_ polyhedron* ptr, int m). The function constructs a new polygon resulting from the intersection of a plane (described by point xc and normal vector nˆ m ) and a 3D planar polygon cell_vert. Fig. 20(a) shows the new polygon constructed as a result of the intersection. To avoid repeated vertices, we consider a line segment between consecutive vertices [i, j) (i included, j not) of the polygon cell_vert and look for its intersection with the plane (xc , nˆ m ). Fig. 20(b) depicts that the point D, which would otherwise be stored twice due to the plane’s intersections with CD and ED, is stored only once following this technique. The newly constructed polygon is stored as the mth face of the polyhedron ptr. The vertices of the constructed face are in CCW direction when looked from outside the polyhedron. 3D rotation of a plane: Sorting vertices of a polyhedron face is computationally expensive, and we have so far avoided it by storing vertices in the CCW direction. However, in the function const_com_polyhedron(), which is to follow,

570

A. Pathak, M. Raessi / Journal of Computational Physics 307 (2016) 550–573

Fig. 20. (a) Polygon HDEFG is formed due to the intersection of a polygon ABCDEF with a plane (xc , nˆ m ). (b) Each line segment [i, j) (i included j not) is considered for the intersection with the plane (xc , nˆ m ). As a result of this technique, the point D is stored only once, due to the intersection of CD with the plane.

ˆ OF is along the z-axis of the coordinate system. The polygon rotates about axis ˆ = nˆ × kˆ to lie on an x– y Fig. 21. A 3D polygon having an orientation n. plane.

constructing the last face of a polyhedron would require sorting vertices in the CCW direction. For this purpose, we rotate a 3D polygon to lie on an x– y plane, perform sorting on the 2D polygon, then rotate back the sorted polygon to the original 3D space. A standard method of rotation in 3D space involves the Euler’s angles, which can be very cumbersome [43]. Quaternions offer a convenient and efficient mechanism of 3D rotation. A detailed description on the concept and operations of quaternions can be found in [43,44]. We will briefly present the mechanism to rotate a general polygon in 3D. Consider a 3D polygon P (can be face of a polyhedron) ABCDE, shown in Fig. 21, with a normal vector nˆ pointing outwards. Our goal is to rotate it to lie on an x– y plane. Amount of rotation is calculated as θ = cos−1 (ˆn · kˆ ) = cos−1 n z , ˆ We then compute the axis ˆ (x ˆi +  y ˆj + z kˆ ) about which rotation will be performed: ˆ = nˆ × k. ˆ where nˆ = nx ˆi + n y ˆj + n z k. The quaternions that would rotate each point on the plane of the polygon P to an x– y plane are defined as [44]:

q = cos(θ/2) + (x ˆi +  y ˆj + z kˆ ) sin(θ/2)

(18)

= cos(θ/2) − (x ˆi +  y ˆj + z kˆ ) sin(θ/2)

(19)

−1

q

The data structure Quat that stores information on a quaternion is presented in the Supporting Information document. Any point ( p x ˆi + p y ˆj + p z kˆ ), on the plane of the polygon P , is first converted to a quaternion p (function point2Quat() in the Supporting Information document) and then rotated to an x– y plane by a quaternion conjugation:

p = qpq−1 = (qp)q−1 = q(pq−1 )

(20)

It is convenient to perform the operation (20) by making use of the associative property of quaternion multiplication instead of constructing a rotation matrix. We first compute the product qp resulting in a new quaternion, which is again multi-

A. Pathak, M. Raessi / Journal of Computational Physics 307 (2016) 550–573

571

Fig. 22. Vertices of a polygon are sorted in CCW direction based on polar angle θ of each vertex. Polar angle calculation is performed in a coordinate system x – y  lying in the plane of the polygon. Centroid, xc , of the polygon serves as the origin O of the coordinate system.

plied with q−1 to yield (qp)q−1 , thus completing operation (20). The function for quaternion multiplication qt_mul() is provided in the Supporting Information document. Point ( p x ˆi + p y ˆj + p z kˆ ), extracted from quaternion p (using the function Quat2point() in the Supporting Information document), will now lie on an x– y plane. The operation (20) is repeated on all vertices of the polygon P , which will rotate them to lie on an x– y plane resulting in a new polygon P  . We now employ a sorting function sort_ccw() (see the Supporting Information document), which uses insertion sort [31] based on polar angle, to sort vertices of the polygon P  in the CCW direction. Polar angle θ , which is shown in Fig. 22, is calculated for each vertex of the polygon using the atan2() function offered by standard math libraries. While sorting, we consider only the x and y components, and disregard the z component. Afterwards, the sorted vertices of the polygon P  could be rotated back to the original plane by performing an inverse quaternion conjugation:

p = q−1 p q

(21)

The last operation (21) can be avoided if instead of storing vertex coordinates during the sort operation on P  , we store the vertex indices with respect to the unsorted polygon. Since the relative position of the vertices remains the same during polygon rotation, the indices could be used to fetch the vertices and create a sorted set from the original polygon P . This approach has been implemented for computational efficiency. complex polyhedron Template: void const_com_polyhedron (st_polyhedron* sim_ptr, st_polyhedron* com_ptr, st_point* xc , st_point* nˆ m , int srt). The function constructs a polyhedron com_ptr when a plane L, described by orientation nˆ m and point xc lying on it, cuts a polyhedron sim_ptr. sim_ptr is a general polyhedron, and vertices of its each face are sorted in the CCW direction. Fig. 19(b) shows the cutting plane L (shown in red) dividing the given polyhedron sim_ptr (the cube) into two volumes. The volume lying below the plane constitutes the new polyhedron com_ptr, whose information is stored. The plane L is checked for intersection with each face of the polyhedron sim_ptr. If intersection occurs, using const_face() (see the Supporting Information documented), a new face lying below the plane L is computed, which then constitutes a new face of the polyhedron com_ptr (see Fig. 20(a)). If there is no intersection and a face of sim_ptr lies completely below the plane L, the face directly becomes a new face of the polyhedron com_ptr. Last face of com_ptr lies on L and is represented by the polygon

Constructing

IJKLMN shown in Fig. 19(b). To construct the last face, we select all the previously constructed faces of the polyhedron com_ptr that were created as a result of the intersection between the plane L and the polyhedron sim_ptr’s faces. We then connect all the vertices created due to the intersections to construct the last face (IJKLMN). To avoid repeated vertices at connected faces, we use a scheme which can be found in the corresponding pseudo-code provided in the Supporting Information document. This scheme exploits the CCW storage of the function const_face(). The last face IJKLMN of the polyhedron com_ ptr is still unsorted. If new polyhedrons are to be constructed from com_ptr, e.g. by cutting it with another plane, the last face would need sorting in CCW direction (function sort_ccw() described in the Supporting Information document). Otherwise, sorting is not required at this stage. The input variable srt is set to 1 or 0 corresponding to the two scenarios, respectively. Because sorting is a computationally expensive procedure, we have kept its usage to a minimum. const_com_polyhedron()’s mechanism is similar to the clipping and capping algorithm of Stephenson and Christiansen [45]. All steps before the construction of the last face correspond to the clipping stage. Construction of the last face constitutes the capping stage. The method presented in Ref. [45] is different from our approach in two regards: (a) we avoid repeated vertices in a polyhedron face during its construction (see Fig. 20(b)), as opposed to hash tables used in [45], and (b) we construct the last face using a sort-based method described in the last paragraph, as opposed to using connectivity information described in [45]. Since we work with small number of faces, performing an insertion based sort is “less complicated to implement and more efficient” as argued by Maric et al. [46].

572

A. Pathak, M. Raessi / Journal of Computational Physics 307 (2016) 550–573

polygon area calculation Template: double area3D_polygon (st_point* vert, st_point* nˆ m , int n_vert). Given n_vert number of sorted (in CCW direction) vertices (st_point* vert) of a 3D polygon, the func-

3D

tion computes the polygon’s area. It uses the fast approach of Sunday [47], which calculates the area of the 3D polygon’s projection on one of x– y, y–z or z–x planes. The projected area is then scaled back using the appropriate normal vector component of the 3D polygon. Calculating polyhedron volume Template: double calc_poly_vol (st_polyhedron* ptr). The function calculates the volume of a polyhedron represented by the pointer ptr. The approach is to decompose the polyhedron into disjoint pyramids, calculate the volumes of individual pyramids, and then add the individual volumes to obtain the polyhedron volume. Appendix B. Supplementary material Supplementary material related to this article can be found online at http://dx.doi.org/10.1016/j.jcp.2015.11.062. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35]

D.L. Youngs, An interface tracking method for a 3D Eulerian hydrodynamics code, Tech. Rep. 44/92/35, AWRE, 1984. M. Bussmann, A three dimensional model of an impacting droplet, Ph.D. thesis, University of Toronto, 2000. B. Choi, M. Bussmann, A piecewise linear approach to volume tracking a triple point, Int. J. Numer. Methods Fluids 53 (2007) 1005–1018. A. Ghasemi, A. Pathak, M. Raessi, Computational simulation of the interactions between moving rigid bodies and incompressible two-fluid flows, Comput. Fluids 94 (2014) 1–13. D.L. Youngs, Time-Dependent Multi-Material Flow with Large Fluid Distortion, Academic, New York, 1982. D.J. Benson, Volume of fluid interface reconstruction methods for multi-material problems, Appl. Mech. Rev. 55 (2002) 151–165. S. Mosso, S. Clancy, A geometrically derived priority system for Young’s interface reconstruction, Tech. Rep. LA-CP-95-0081, Los Alamos National Laboratory, Los Alamos, NM, 1995. D.J. Benson, Eulerian finite element methods for the micromechanics of heterogeneous materials: dynamic prioritization of material interfaces, Comput. Methods Appl. Mech. Eng. 151 (3–4) (1998) 343–360, containing papers presented at the Symposium on Advances in Computational Mechanics. D. de Niem, E. Kührt, U. Motschmann, A volume-of-fluid method for simulation of compressible axisymmetric multi-material flow, Comput. Phys. Commun. 176 (3) (2007) 170–190. C.D. Sijoy, S. Chaturvedi, Volume-of-fluid algorithm with different modified dynamic material ordering methods and their comparisons, J. Comput. Phys. 229 (10) (2010) 3848–3863. A. Caboussat, M.M. Francois, R. Glowinski, D.B. Kothe, J.M. Sicilian, A numerical method for interface reconstruction of triple points within a volume tracking algorithm, Math. Comput. Model. 48 (2008) 1957–1971. S.P. Schofield, R.V. Garimella, M.M. Francois, R. Loubère, A second-order accurate material-order-independent interface reconstruction technique for multi-material flow simulations, J. Comput. Phys. 228 (2009) 731–745. J.E. Pilliod Jr., E.G. Puckett, Second-order accurate volume-of-fluid algorithms for tracking material interfaces, J. Comput. Phys. 199 (2) (2004) 465–502. J.-P. Braeunig, B. Desjardins, J.-M. Ghidaglia, A totally Eulerian finite volume solver for multi-material fluid flows, Eur. J. Mech. B, Fluids 28 (4) (2009) 475–485. R. Loubere, J.-P. Braeunig, J.-M. Ghidaglia, A totally Eulerian finite volume solver for multi-material fluid flows: Enhanced Natural Interface Positioning (ENIP), Eur. J. Mech. B, Fluids 31 (2012) 1–11. B. Blais, J.-P. Braeunig, D. Chauveheid, J.-M. Ghidaglia, R. Loubère, Dealing with more than two materials in the FVCF–ENIP method, Eur. J. Mech. B, Fluids 42 (2013) 1–9. H.T. Ahn, M. Shashkov, Multi-material interface reconstruction on generalized polyhedral meshes, J. Comput. Phys. 226 (2007) 2096–2132. V. Dyadechko, M. Shashkov, Reconstruction of multi-material interfaces from moment data, J. Comput. Phys. 227 (11) (2008) 5361–5384. H.T. Ahn, M. Shashkov, Adaptive moment-of-fluid method, J. Comput. Phys. 228 (8) (2009) 2792–2821. S. Galera, J. Breil, P.-H. Maire, A 2D unstructured multi-material Cell-Centered Arbitrary Lagrangian–Eulerian (CCALE) scheme using MOF interface reconstruction, in: 10th ICFD Conference Series on Numerical Methods for Fluid Dynamics, Comput. Fluids 46 (1) (2011) 237–244. J. Breil, T. Harribey, P.-H. Maire, M. Shashkov, A multi-material ReALE method with MOF interface reconstruction, Comput. Fluids 83 (2013) 115–125. M. Kucharik, M. Shashkov, Conservative multi-material remap for staggered multi-material arbitrary Lagrangian–Eulerian methods, J. Comput. Phys. 258 (2014) 268–304. K.S. Bonnell, M.A. Duchaineau, D.R. Schikore, B. Hamann, K.I. Joy, Material interface reconstruction, IEEE Trans. Vis. Comput. Graph. 9 (4) (2003) 500–511. J. Meredith, Material interface reconstruction in visit, Tech. Rep. UCRL-CONF-209351, Lawrence Livermore National Laboratory, Livermore, CA, United States, 2005. J.C. Anderson, C. Garth, M.A. Duchaineau, K. Joy, Discrete multi-material interface reconstruction for volume fraction data, in: Proc. of Eurographics/IEEE-VGTC Symposium on Visualization, Computer Graphics Forum 27 (3) (2008). W.J. Rider, D.B. Kothe, Reconstructing volume tracking, J. Comput. Phys. 141 (2) (1998) 112–152. G. Son, A level set method for incompressible two-fluid flows with immersed solid boundaries, Numer. Heat Transf., Part B, Fundam. 47 (5) (2005) 473–489. M. Fujita, O. Koike, Y. Yamaguchi, Computation of capillary interactions among many particles at free surface, Appl. Phys. Express 6 (3) (2013) 036501. R.I. Saye, J.A. Sethian, The Voronoi implicit interface method for computing multiphase physics, Proc. Natl. Acad. Sci. 108 (49) (2011) 19498–19503. R. Saye, J.A. Sethian, Analysis and applications of the Voronoi implicit interface method, J. Comput. Phys. 231 (18) (2012) 6051–6085. W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in C, 2nd edn., The Art of Scientific Computing, Cambridge University Press, New York, NY, USA, 1992. ´ C. Tropea, Dynamic contact angle of spreading droplets: experiments and simulations, Phys. Fluids 17 (6) Š. Šikalo, H.-D. Wilhelm, I. Roisman, S. Jakirlic, (2005) 062103. Y. Sui, P.D. Spelt, An efficient computational model for macroscale simulations of moving contact lines, J. Comput. Phys. 242 (2013) 37–52. K. Yokoi, D. Vadillo, J. Hinch, I. Hutchings, Numerical studies of the influence of the dynamic contact angle on a droplet impacting on a dry surface, Phys. Fluids 21 (7) (2009) 072102. M. Bussmann, J. Mostaghimi, S. Chandra, On a three-dimensional volume tracking model of droplet impact, Phys. Fluids 11 (6) (1999) 1406–1417.

A. Pathak, M. Raessi / Journal of Computational Physics 307 (2016) 550–573

573

[36] M. Raessi, J. Mostaghimi, M. Bussmann, Advecting normals vectors, a new method for calculating interface normals and curvatures when modeling two-phase flows, J. Comput. Phys. 226 (2007) 774–797. [37] S.J. Cummins, M.M. Francois, D.B. Kothe, Estimating curvature from volume fractions, Comput. Struct. 83 (6) (2005) 425–434. [38] D.J. Harvie, D.F. Fletcher, A new volume of fluid advection algorithm: the stream scheme, J. Comput. Phys. 162 (1) (2000) 1–32. [39] S. Afkhami, M. Bussmann, Height functions for applying contact angles to 3D VOF simulations, Int. J. Numer. Methods Fluids 61 (8) (2009) 827–847. [40] M. Sussman, A second order coupled level set and volume-of-fluid method for computing growth and collapse of vapor bubbles, J. Comput. Phys. 187 (1) (2003) 110–136. [41] J. Helmsen, P. Colella, E.G. Puckett, Non-convex profile evolution in two dimensions using volume of fluids, Tech. Rep. LBNL-40693, Lawrence Berkeley National Laboratory, 1997. [42] X. Lv, Q. Zou, Y. Zhao, D. Reeve, A novel coupled level set and volume of fluid method for sharp interface capturing on 3D tetrahedral grids, J. Comput. Phys. 229 (7) (2010) 2573–2604. [43] K. Shoemake, Animating Rotation with Quaternion Curves, in: SIGGRAPH ’85 Proceedings of 12th Annual Conference on Computer Graphics and Interactive Techniques, vol. 9, ACM, 1985. [44] J.B. Kuipers, Quaternions and Rotation Sequences, vol. 66, Princeton University Press, Princeton, 1999. [45] M.B. Stephenson, H.N. Christiansen, A polyhedron clipping and capping algorithm and a display system for three dimensional finite element models, SIGGRAPH Comput. Graph. 9 (3) (1975) 1–16. [46] T. Maric, H. Marschall, D. Bothe, voFoam-A geometrical volume of fluid algorithm on arbitrary unstructured meshes with local dynamic adaptive mesh refinement using OpenFOAM, arXiv preprint, arXiv:1305.3417. [47] D. Sunday, Fast polygon area and Newell normal computation, J. Graph. Tools 7 (2) (2002) 9–13.