A coupled volume-of-fluid and level set method based on general curvilinear grids with accurate surface tension calculation

Journal of Computational Physics 396 (2019) 799–818

Contents lists available at ScienceDirect

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

A coupled volume-of-fluid and level set method based on general curvilinear grids with accurate surface tension calculation Zhizhu Cao a , Dongliang Sun b,∗ , Jinjia Wei a,c , Bo Yu b,∗ , Jingfa Li b a b c

School of Chemical Engineering and Technology, Xi’an Jiaotong University, Xi’an 710049, China School of Mechanical Engineering, Beijing Institute of Petrochemical Technology, Beijing 102617, China State Key Laboratory of Multiphase Flow in Power Engineering, Xi’an Jiaotong University, Xi’an, 710049, China

a r t i c l e

i n f o

Article history: Received 11 January 2019 Received in revised form 5 July 2019 Accepted 8 July 2019 Available online 11 July 2019 Keywords: Two-phase flow Level set Volume-of-fluid Interface normal Curvilinear grid Sharp interface

a b s t r a c t In this paper, a coupled volume-of-fluid (VOF) and level set (LS) method on general curvilinear grids, combining the advantages of VOF method and LS method is developed for interfacial flow simulations. In this method, an iterative geometric operation is introduced to calculate the level set function on physical domain and computational domain. In particular, a direct construction of the distance function is successfully implemented on the physical domain and thus the numerical error induced by contravariant transformation from the computational domain could be avoided. The level set based continuum surface tension (CSF) model on curvilinear grids is derived to model surface tension accurately, and then the model as well as the proposed coupled volume tracking method is incorporated into an incompressible Navier–Stokes solver. Finally, four benchmark problems including time-reversed single vortex flow, static droplet, dam break and single bubble rising problems are tested on irregular domains to validate the accuracy and robustness of the proposed method to deal with complex domains. The results show good agreement with previous experimental data or numerical results, good characteristics of mass conservation are also observed, indicating that the present method can accurately simulate incompressible two-phase flows. © 2019 Elsevier Inc. All rights reserved.

1. Introduction Unsteady two-phase flows in irregular domain receive special attention in computational fluid dynamics, because it involves the challenges not only in dealing with moving surfaces but also the mesh generation in complex domain. Most researches mainly focus on unstructured meshes to deal with complex domain [1]; an alternative is using the body-fittedcoordinate (BFC) technology. The main idea of this technology is that an arbitrary geometry of interest indicated as the physical domain, is mapped onto a rectangular region called computational domain in the BFC system. Accordingly, the solution in the physical domain can be transformed from a curvilinear grid system to the orthogonal computational domain [2]. Since the general curvilinear grid is one kind of structured grids, it possesses advantages of simplicity in data structure, fast speeds in grid generation, easy implementation of high order numerical scheme and good convergence prop-

*

Corresponding authors. E-mail addresses: [email protected] (D. Sun), [email protected] (B. Yu).

https://doi.org/10.1016/j.jcp.2019.07.016 0021-9991/© 2019 Elsevier Inc. All rights reserved.

800

Z. Cao et al. / Journal of Computational Physics 396 (2019) 799–818

erty. In the following, literature reviews on moving surfaces tracking methods such as volume-of-fluid and level set on general curvilinear grid are discussed. It is well known that VOF method has an excellent mass conservation property [3]. However, for the reason that the volume fraction function is a discontinuous function, it is difficult to obtain the accurate normal and curvature. Kothe et al. [4] developed a VOF method on the general curvilinear grid where the interface was reconstructed on the physical domain. Nikseresht et al. [5] proposed a 3D VOF method with Lagrangian propagation in general curvilinear coordinates and applied the method to solve viscous incompressible free surface flows around air-cushion vehicles (ACV). López and Hernandez [6] presented an analytical and geometrical tool for 3D VOF methods on a general grid system. Sarthou et al. [7] put forward a global methodology coupled with an immersed boundary method or with phase advection methods such as VOF-PLIC (piecewise linear interface construction), Front-tracking approaches to deal with fictitious domains of all kinds on curvilinear grids. They found that any operations can be performed much faster on orthogonal computational grids than on curvilinear grids after transformation. It is note that recently Nguyen and Park [8,9] presented a 3D algebraic VOF method on the general curvilinear grid for complex two-phase incompressible flows. In their method, an interface-sharpening equation is derived by combining the artificial compression and anti-diffusion terms. There is no need to reconstruct interfaces geometrically, however, compared to the geometric VOF method, mass losses due to numerical diffusion in algebraic VOF method are still obvious since the reported mass loss error is as high as 2% (see Fig. 4 in [8]). In the LS method [10], the moving surface is implicitly captured by the zero level set. The curvature can be computed accurately by using the level set function, which is a smooth and continuous function. Yue et al. [11] presented a LS method that was coupled with incompressible Navier–Stokes equations in a generalized curvilinear coordinate system for the study of free surface flows. Different numerical schemes for advection of the level set function are compared within a generalized curvilinear framework, including the third order quadratic upwind interpolation for convective kinematics (QUICK) scheme, and the second and third order essentially non-oscillatory (ENO) schemes. Similar to other grids systems (such as the Cartesian grid, unstructured grid), the major drawback of LS method on curvilinear grid system is its poor mass conservation property, i.e., loss/gain of mass introduced by numerical dissipation. To overcome this issue, Wang et al. [12] developed a mass conservative level set method based on a curvilinear coordinate system, and applied it to simulate the spilling breaking wave problem, they found that there are significant impacts of temporal and spatial schemes on computational results. To combine the advantages of the volume-of-fluid method with the level set method, Sussman and Puckett [13] developed a coupled level set/volume-of-fluid (CLSVOF) method for capturing the free surface in two-phase flow problems. Subsequently, several improved CLSVOF methods [14–16] were presented. It is also noted that most of the researches focused on the extension of the CLSVOF method to unstructured grids for dealing with complex domain [17,18]. However, it is hard to find any related literature about the implementation of the CLSVOF on curvilinear grids. Although CLSVOF method is more accurate than either the VOF or LS method alone, it is more complicated because both the volume fraction and level set advection equations need to be solved in this method. Recently, Sun and Tao [19] presented a coupled volume-of-fluid and level set (VOSET) method, which is simpler than the CLSVOF method. In the method, only the volume fraction advection equation needs to be solved and the LS function is computed by a simple iterative geometric operation. The main idea of this operation is that the shortest distance from any point to the interface can be computed directly in a geometric manner. After that, the VOSET method was extended to the 3D structured grids by Ling et al. [20] and to unstructured grids by the present authors Cao et al. [21–23]. It is worth mentioning that Wang et al. [24] developed a second-order volume-of-fluid method (VOF) on general structured grids by using a reconstructed distance function (RDF) technique. In their method, the distance function on the computational domain is constructed geometrically following the procedure given in [25]. The distance on the physical domain is obtained by contravariant transformation from the calculated distance on the computational domain. It is noted that the possible error could be induced due to the contravariant transformation. Moreover, several geometric configurations for different cases requiring special treatment in the transformation seem very complicated. To the best of our knowledge, there is no report on direct construction of the distance function on the physical domain as “it is very difficult since only the reconstructed interface on the computational domain is known [24]”. In this paper, we present a coupled volume-of-fluid and level set method on general curvilinear grids for simulating the two-phase flows in complex irregular domains. In this method, an iterative geometric operation is extended to curvilinear grids for calculating the level set function φ near interfaces, which can be applied to compute the accurate curvature κ and smooth the discontinuous physical quantities near interfaces. It is worth mentioning that the level set function on the physical domain is computed by direct construction of the distance function rather than contravariant transformation as in [24]. Afterwards a PLIC method is developed to capture gas-liquid interfaces on curvilinear grids, which can conserve the mass and overcome the disadvantage of non-conservation of mass in LS method. The structure of the paper is as organized as follows. The solution procedure of the VOSET method on general curvilinear grids is first discussed in Section 2. The proposed volume tracking method is validated and then incorporated into an incompressible flow solver in Section 2 as well, in which the level set based CSF model on curvilinear grids is also derived and implemented for dealing with the surface tension dominated problems. Subsequently, several classical tests i.e. static droplet, dam break and single bubble rising problems, are conducted to verify the proposed method in Section 3. Some conclusions are made in Section 4.

Z. Cao et al. / Journal of Computational Physics 396 (2019) 799–818

801

Fig. 1. Grids on the physical domain and the computational domain.

2. VOSET on general curvilinear grids 2.1. Level set function It is noted that the results could be significantly affected by the initial volume fraction field. To eliminate this possible negative influence and calculate the initial volume fraction accurately, an algorithm based on quad-tree grids that can subdivide an interface cell is implemented in our code work. Specifically, take a round interface for example, if all vertexes of a quadrilateral cell on curvilinear coordinates are outside the interface, then the volume fraction α = 0; if inside, then α = 1; otherwise, the cell is an interface cell with 0 < α < 1. For an interface cell, its volume fraction can be computed by using divisions of the cell as follows. first, let us insert and connect the middle point of each edge of a quadrangle, then the cell turns into four new cells with approximately equal areas. After subdivisions with the generation level NL, the parent interface cell can be divided into 4NL new cells. It is easy to obtain the number Q of quadrilateral cells inside the interface by judging whether their locations of center points are inside the interface. Finally, α = Q /4NL . It is found that this initialization method is simple, robust and high accurate for quadrilateral meshes. In the paper, seven refinement levels are used in the initialization in each interface cell. In PLIC method, the interface in a cell is represented by a line segment and is determined by the volume fraction function α together with a normal vector. The normal vector is calculated either by the volume fraction or by the level set function. It is worth mentioning that the level set function φ p on the physical domain is not a level set function on the computational domain (see Fig. 1). Considering the computational domain is regular and orthogonal, the level set function φc on the computational domain can be obtained following the iterative geometric operation given by Sun and Tao [19]. A brief review on the implementation process of the iterative geometric operation is presented as follows: 1. Solve the interface normal vector by using volume fraction α and reconstruct interface with the analytical PLIC-VOF approach. 2. Set the initial value of φ in the whole computational domain.



φi0, j =

−max( L , W ), max( L , W ),

αi, j ≥ 0.5 αi, j < 0.5

(1)

where the subscript (i , j) refers to the cell index, L and W represent the length and width of the domain, respectively. 3. Mark the cells near interfaces. This step can save computational cost significantly. The marked cells is recommended to be within 3h away from center points of the interfacial cells, which is wide enough to compute the interface normal and curvature. 4. Calculate the level set function φ in the marked region. This step obtains the shortest distance between the center point of each cell in the marked region and the interface geometrically. 5. Reconstruct the interface again with the more accurately interface normal vector by using φ , repeat steps 2-4 until the iteration times reach the preset valve. Unlike Wang et al.’s contravariant transformation method [24], a direct construction of the distance function geometrically on the physical domain is proposed to obtain the level set function. For 2D cases, the interface in PLIC is a line segment with its starting point and ending point, for 3D cases, the interface is a polygon with several corners. When the interface on the computational domain is reconstructed by using φc and α , the interface polygon corners coordinates are obtained accordingly. For 2D cases, the interface corners coordinates on physical domain can be calculated based on linear interpolations as follows

  (ξ − ξ1 ) (η − η1 ) P ξ,η = ( P ξ ,η + P ξ2 ,η2 − P ξ2 ,η1 − P ξ1 ,η2 ) + P ξ2 ,η1 − P ξ1 ,η1 (ξ2 − ξ1 ) (η2 − η1 ) 1 1 (η − η1 ) + ( P ξ ,η − P ξ1 ,η1 ) + P ξ1 ,η1 (η2 − η1 ) 1 2

(2)

where P ξ,η is the coordinates of point (x, y ) on the physical domain corresponding to the coordinates of point (ξ, η) on the computational domain.

802

Z. Cao et al. / Journal of Computational Physics 396 (2019) 799–818

Fig. 2. The minimum distance from cell central point P (i , j ) to an interface EF in a cell.

When the interface is constructed with its corners coordinates, the level set function φ p in the flagged region on the physical domain can be computed. To be specific, by comparing all of the minimum distances from cell central point P (i , j ) to any interface within the 6 h × 6 h stencil around the cell (i , j ), the shortest distance di , j at cell (i , j ) is obtained. Fig. 2 gives three possible cases of the minimum distance from cell central point P (i , j ) to an interface EF in a cell. If θ1 > 90◦ , PE segment length is the minimum distance; if θ2 > 90◦ , PF is the minimum distance; else if θ1 < 90◦ and θ2 < 90◦ , the vertical line PM is the minimum distance. Then, the level set function φ P can be calculated by

⎧ ⎨ −di , j , αi , j > 0.5 αi, j = 0.5 φ i , j = 0, ⎩ di , j , αi, j < 0.5

(3)

Finally, the more accurate normal vector to the interface can be calculated by using the above level set function on physical domain

 = ∇φ/|∇φ| n

(4)

The above distance function gradient is discretized following Wang et al.’s treatment [24] in which the derivative of distance function in each coordinate direction is approximated using the forward, backward, and central difference schemes, respectively. Take the ξ direction for example,

(nξ )1 = (φi , j − φi −1, j )/ξ, (nξ )2 = (φi +1, j − φi , j )/ξ,

(5)

(nξ )3 = (φi +1, j − φi −1, j )/2ξ. For convenience, ξ = η = 1 in this study. Among the 9 combinations for 2D cases, the one making ||∇φ| − 1| the minimum will be chosen to calculate the normal vector in Eq. (4). In summary, the flow chart of the geometric operation on computational and physical domains is given in Fig. 3. To verify the present method of calculating the level set function, a circular interface with radius of 0.15 is located at (0.5, 0.75) in a unit circular domain. The Laplace equation is used for body-fitted coordinate grids generation. Fig. 4 depicts the signed distance function near the circular interface (Fig. 4a) and the corresponding error contour (Fig. 4b) with three different mesh densities (h = [1/25, 1/50, 1/100]) calculated by the present method. In Fig. 4a the red thick lines denote the interfaces and the black thin lines represent the contour lines of the signed distance function. It can be seen that the accurate level set function can be obtain by the proposed iterative geometrical operation for curvilinear grids. The error contours in Fig. 4b show that the accuracy of constructed distance function is positively correlated with gird orthogonality. The error L 1 is defined by

L1 =

M 1 

M

|dcal − dexact |

(6)

i =1

where M refers to the number of the interfacial cells along the circumference of the circle, dcal and dexact represent the calculated distance and the exact distance solutions, respectively. The convergence of error L 1 on the level set function calculation for the circular interface with respect to increasing number of iterations N is shown in Table 1. It is found that the error obviously decreases from N = 1 to N = 2, while there is no change or very little change if N is greater than 2. 2.2. Volume fraction advection In the VOF method, the interface in a grid cell is determined by volume fraction α , which is a value between 0 and 1. Specifically, if the cell is fulfilled by the target fluid, then α is equal to 1, if the cell is fulfilled by the non-target fluid, then

Z. Cao et al. / Journal of Computational Physics 396 (2019) 799–818

803

Fig. 3. Flow chart of the geometric operation on computational and physical domains. Table 1 The convergence of error L 1 on the level set function calculation for a circular interface (radius 0.15) at different mesh densities with respect to increasing N. Grid

Total cells

Flagged cells

VOSET on curvilinear grid N =1

N =2

N =3

25 × 25 50 × 50 100 × 100 200 × 200 300 × 300

625 2500 10000 40000 90000

32 68 134 294 446

8.51 × 10−4 1.93 × 10−4 5.21 × 10−5 1.61 × 10−5 1.03 × 10−5

8.32 × 10−4 1.84 × 10−5 4.82 × 10−5 1.10 × 10−5 5.11 × 10−6

8.32 × 10−4 1.84 × 10−5 4.82 × 10−5 1.09 × 10−5 5.12 × 10−6

α is equal to 0, otherwise, the cell is an interface cell and α represents the volume ratio of target fluid. It can be seen that the volume fraction on physical and computational domains is the same. For incompressible flows, the advection equation of volume fraction α can be written as

∂α + u · ∇ α = 0 ∂t

(7)

Eq. (7) can be transformed to the following expression

∂u j ∂ α ∂(α u j ) + −α =0 ∂t ∂xj ∂xj

(8)

The governing equation (Eq. (8)) in the predefined Cartesian coordinate system (x, y ) is transformed into an equation defined in the general curvilinear coordinate system (ξ, η) (see Fig. 1) as given below

∂( αJU ) ∂( αJV ) ∂( UJ ) ∂( VJ ) ∂α +( ) = α( ) + + ∂t ∂ξ ∂η ∂ξ ∂η

(9)

where J is Jacobian determinant, U and V are the contravariant velocities, J = xξ y η − xη y ξ , U = u y η − vxη , V = vxξ − u y ξ . The right term is zero if the continuum equation is satisfied, and then the discretized equation of Eq. (9) is given by

(α

n +1 ∗ ) i, j

=α

n i, j

+

1

ξ η



F in−1/2, j J i −1/2, j

−

F in+1/2, j J i +1/2, j

+

F in, j −1/2 J i , j −1/2

−

F in, j +1/2 J i , j +1/2



(10)

where F in−1/2, j , F in+1/2, j , F in, j −1/2 , F in, j +1/2 is the area of the volume fraction flux polygon passing through west, east, south and north faces of cell (i , j ) during time step t at time level n, respectively. The volume fraction flux polygon can be computed after the interface is reconstructed by PLIC algorithm, in which 16 possible cases for the interface shape on

804

Z. Cao et al. / Journal of Computational Physics 396 (2019) 799–818

Fig. 4. The level set function near the circular interface and the corresponding error contour with different mesh densities calculated by the iterative geometrical operation. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)

rectangular grids must be considered [19]. Finally, the interface is propagated by updating the volume fraction value in the whole computational domain. However, because discrete velocity divergences are not exactly zero due to the numerical error, the right term in Eq. (9) is considered as a divergence correction to decrease errors induced by the mass residual

(α

n +1 ) i, j

= (α

n +1 ∗ ) i, j

t + ξ η



U in−1/2, j J i −1/2, j

η −

U in+1/2, j J i +1/2, j

η +

V in, j −1/2 J i , j −1/2

ξ −

V in, j +1/2 J i , j +1/2

ξ

(11)

Z. Cao et al. / Journal of Computational Physics 396 (2019) 799–818

805

2.3. Governing equations and their solutions 2.3.1. Governing equations The governing equations on the computational domain for transient, incompressible, Newtonian, two-phase flows including the continuum and momentum equations with surface tension and gravity, can be expressed as

1 ∂U

1 ∂V

+

=0

(12)

∂ u ˆ c ˆ v ˆ +E =E +S ∂t

(13)

J ∂ξ

J ∂η

Eˆ c , Eˆ v , Sˆ in the momentum equations represent the convective term, diffusive term and source term, respectively. Diffusive term:

Eˆ v =

1

∇ · μ(φ)∇ u      1 1 ∂ μ(φ)  1 ∂ μ(φ)  = α (u )ξ − β(u )η + γ (u )η − β(u )ξ ρ (φ) J ∂ξ J J ∂η J

ρ (φ)

(14)

where α = x2η + y 2η , γ = x2ξ + y 2ξ , β = xξ xη + y ξ y η . Convective term:

Eˆ c =

) 1 ∂(U u J

∂ξ

+

) 1 ∂( V u

(15)

∂η

J

Source term:

Sˆ = −

1

ρ (φ)

∇p +

1

ρ (φ)

+ F sv + g

1

ρ (φ)

Eˆ 0v

 ∇ p = ξx p ξ + ηx p η , ξ y p ξ + η y p η  )T Eˆ 0v = ∇ · μ(φ)(∇ u   ∂u ∂ ∂v ∂ ∂u ∂ ∂v ∂ = (μ(φ) ) + (μ(φ) ), (μ(φ) ) + (μ(φ) ) ∂x ∂x ∂y ∂x ∂x ∂y ∂y ∂y

(16)

(17)

If x = x(ξ, η), y = y (ξ, η), ξ = ξ(x, y ), η = η(x, y ), then

ξx =

1 J

yη ,

1

1

1

J

J

J

ηx = − y ξ , ξ y = − xη , η y = xξ

(18)

We have

   ∂  ∂ ∂u = μ(φ) μ(φ)(u ξ ξx + u η ηx ) ∂x ∂x ∂x     = μ(φ)(u ξ ξx + u η ηx ) ξ ξx + μ(φ)(u ξ ξx + u η ηx ) η ηx       yη ∂ yη yη yξ yξ ∂ yξ = μ(φ) μ(φ) − uξ − uη uξ − uη J ∂ξ J J J ∂η J J    ∂ ∂  ∂v = μ(φ) μ(φ)( v ξ ξx + v η ηx ) ∂y ∂x ∂y       xη ∂ yη yη yξ xξ ∂ yξ μ(φ) μ(φ) + vξ − vη vξ − vη =− J ∂ξ J J J ∂η J J           yη ∂ yη yη yξ yξ ∂ yξ ∂ ∂u ∂v ∂ + = − μ(φ) μ(φ) μ(φ) μ(φ) uξ − uη uξ − uη ∂x ∂x ∂y ∂x J ∂ξ J J J ∂η J J       xη ∂ yη yη yξ xξ ∂ yξ μ(φ) μ(φ) + − vξ − vη vξ − vη J ∂ξ J J J ∂η J J    ∂  ∂v ∂ = μ(φ) μ(φ)( v ξ ξ y + v η η y ) ∂y ∂y ∂y     = μ(φ)( v ξ ξ y + v η η y ) ξ ξ y + μ(φ)( v ξ ξ y + v η η y ) η η y       xη ∂ xη xη xξ xξ ∂ xξ =− μ(φ) − v ξ + v η + μ(φ) − v ξ + v η J ∂ξ J J J ∂η J J

(19)

(20)

(21)

(22)

806

Z. Cao et al. / Journal of Computational Physics 396 (2019) 799–818





 ∂  μ(φ)(u ξ ξ y + u η η y ) ∂x     = μ(φ)(u ξ ξ y + u η η y ) ξ ξx + μ(φ)(u ξ ξ y + u η η y ) η ηx       yη ∂ xη xη xξ yξ ∂ xξ = μ(φ) − u ξ + u η − μ(φ) − u ξ + u η J ∂ξ J J J ∂η J J     ∂ ∂u ∂v ∂ + μ(φ) μ(φ) ∂x ∂y ∂y ∂y       yη ∂ xη xη xξ yξ ∂ xξ μ(φ) − u ξ + u η − μ(φ) − u ξ + u η = J ∂ξ J J J ∂η J J       xη ∂ xη xη xξ xξ ∂ xξ μ(φ) − v ξ + v η + μ(φ) − v ξ + v η − J ∂ξ J J J ∂η J J ∂u ∂ μ(φ) ∂x ∂y

=

(23)

(24)

The density and viscosity in Eq. (13) can be computed by using the smoothed Heaviside function

ρ (φ) = (1 − H (φ)) ρ g + H (φ)ρl

(25)

μ(φ) = (1 − H (φ)) μ g + H (φ)μl

(26)

where the smoothed Heaviside function is written as

H (φ) =

⎧ ⎪ ⎨ 0, ⎪ ⎩



πφ 1 + ε + π sin( ε ) ,

1 2

φ

1

if φ < −ε if |φ| ≤ ε

(27)

if φ > ε

1,

It is worth pointing out that ε = 1.5h in this study, where h refers to the grid size on the physical domain. It is noted that the grid on physical domain is not uniform, and thus the gird size should vary with coordinates. In the present method, the longest distance from the central points of an interfacial cell to any central points of its nearby cells is chosen as the grid size. 2.3.2. The level set based CSF model on curvilinear grids The surface tension force F sv is computed by the LS function based on the CSF model [26,27] and its formula on the computational domain is derived as follows. The surface tension force on physical domain can be expressed as

F sv = σ κ (φ)δε (φ)∇φ

(28)

where φ represents the real level set function on physical domain, and function is obtained from

δε (φ) =

dH (φ) dφ



1 2ε

=



πφ



1 + cos( ε ) ,

0,

σ is surface tension coefficient. The smoothed delta

if |φ| ≤ ε

(29)

if |φ| > ε

The gradient of the LS function can be calculated by

∇φ =

φξ y η − φη y ξ  J

i+

φη xξ − φξ xη  J

j

(30)

κ (φ) formula derivation is listed as below     ∂ ∂φ 1 ∇φ ∂ ∂φ 1 + κ (φ) = ∇ · = |∇φ| ∂ x ∂ x |∇φ| ∂ y ∂ y |∇φ| ⎧ ⎧ ⎫ ⎫ ⎪ ⎪ ⎪ ⎪ ⎨ ⎨ ⎬ ⎬ φ ξ y η − φη y ξ φ ξ y η − φη y ξ 1 = (  yη −  yξ ) ⎪ J ⎪ ⎩ α (φξ )2 + γ (φη )2 − 2βφξ φη ⎪ ⎩ α (φξ )2 + γ (φη )2 − 2βφξ φη ⎪ ⎭ ⎭ ξ η ⎧ ⎧ ⎫ ⎫ ⎪ ⎪ ⎪ ⎪ ⎨ ⎨ ⎬ ⎬ φη xξ − φξ xη φη xξ − φξ xη 1 xξ −  xη ) + (  ⎪ J ⎪ ⎩ α (φξ )2 + γ (φη )2 − 2βφξ φη ⎪ ⎩ α (φξ )2 + γ (φη )2 − 2βφξ φη ⎪ ⎭ ⎭

Accordingly, the curvature

η

where

α = xη + y η , γ = xξ + y ξ , β = xξ xη + y ξ y η . 2

2

2

2

ξ

(31)

Z. Cao et al. / Journal of Computational Physics 396 (2019) 799–818

807

As shown in Fig. 2, φξ on east and west boundaries, φη on north and south boundaries are approximated by using the second order central difference scheme as follows

(φξ )e = (φ E − φ P )/ξ (φξ ) w = (φ P − φW )/ξ (φη )n = (φ N − φ P )/η (φη )s = (φ P − φ S )/η

(32)

In contrast, φη on east and west boundaries, φξ on north and south boundaries in Eq. (31) are obtained by using a distance-weighted linear interpolation from its neighboring points, respectively.

(φη )e = (1 − λ E →e )(φη ) E + λ E →e (φη ) P (φη ) w = (1 − λ P → w )(φη ) P + λ P → w (φη )W (φξ )n = (1 − λ N →n )(φξ ) N + λ N →n (φξ ) P (φξ )s = (1 − λ P →s )(φξ ) P + λ P →s (φξ ) S

(33)

where the distance-weighted coefficient is defined as

λ E →e = d E →e /(d E →e + d P →e ) λ P → w = d P → w /(d P → w + d W → w ) λ N →n = d N →n /(d N →n + d P →n ) λ P →s = d P →s /(d P →s + d S →s )

(34)

Similar to Eq. (32), the partial derivatives of level set function at central points in Eq. (33) are approximated by the central difference

(φη ) E = (φ N E − φ S E )/2η (φη ) P = (φ N − φ S )/2η (φη )W = (φ N W − φ S W )/2η (φξ ) N = (φ N E − φ N W )/2ξ (φξ ) P = (φ E − φW )/2ξ (φξ ) S = (φ S E − φ S W )/2ξ

(35)

To verify the accuracy of the present method on curvature calculation and make a comparison with other state-of-the-art volume tracking methods on curvilinear grids such as the second order VOF-RDF method [24], the curvature for a circular interface with a radius of 0.3 centered at (0.5, 0.207) is estimated with different grid resolutions. The grid is generated as follows

xi , j = (i − 1)/nx y i, j =

√

2( j − 1)/ny − (i − 1)/nx

(36)

where nx and ny denote grid numbers in i and j directions, respectively. The grid number changes from 32 × 32 to 512 × 512 with a factor of 2. The signed distance functions computed by the proposed method are shown in Fig. 5, it can be seen that smooth contour lines are obtained. For sake of comparison, the L 2 error norm for the curvature is defined by

  M 1  L2 =  (κc − κe )k2 M

(37)

k =1

where κc and κe are computed and exact curvatures in the interfacial cells, respectively, and M is the number of the interfacial cells. The L 2 error norm for the curvature is shown in Table 2, which includes the results computed by the present method and those by the second-order VOF-RDF method in [24]. It can be noted that the relative curvature error by the VOSET method, which varies from 1.94% to 4.95%, is much smaller than that varies from 8.1% to 12.9% of the VOF-RDF method. However, it is worth mentioning that both methods produce the curvature with no more than first order accuracy. This is because the curvatures in both methods are calculated by taking the second-order derivative of the constructed distance function that is second-order accuracy.

808

Z. Cao et al. / Journal of Computational Physics 396 (2019) 799–818

Fig. 5. The signed distance function for the curvatures calculation (red bold line represents the interface). Table 2 The convergence of error L 2 between the computed curvature using the level set function and exact curvature. Method

Second order VOF-RDF [24]

VOSET

Grid

Absolute error

Relative error

Absolute error

Relative error

32 × 32 64 × 64 128 × 128 256 × 256 512 × 512

– 0.375 0.337 0.281 0.431

– 11.2% 10.1% 8.4% 12.9%

0.107 0.0648 0.0943 0.106 0.165

3.21% 1.94% 2.82% 3.18% 4.95%

Fig. 6. Physical domain of the time reversed single vortex flow problem in a circle.

2.3.3. Discretization and solution of governing equations The governing equations are discretized using a collocated finite volume method based on the curvilinear coordinate system. The diffusion flux through cell faces is discretized by using a central difference scheme. The SIMPLE algorithm [28] is adopted to decouple the velocity and pressure (Eqs. (12)–(13)). The deferred-correction power law scheme is used to approximate the convective term. To assure the numerical stability, the time step t at time level n is restricted by the CFL condition and can be computed by



t = C F L × min

h



u n 

where the constantCFL is a value between zero and unity, in this paper, CFL = 0.1.

(38)

Z. Cao et al. / Journal of Computational Physics 396 (2019) 799–818

809

3. Results and discussion In order to verify the accuracy, efficiency, and capability of the present method, three numerical experiments, including the classical single vortex flow [29], static droplet [30], dam break [31] and single bubble rising problems [32] are conducted in this section. 3.1. Time-reversed single vortex flow The time-reversed single vortex flow problem is frequently used to validate the accuracy of a volume tracking method. In this problem, an initially circular fluid body with radius R = 0.15 m is usually positioned at (0.5 m, 0.75 m) in a unit square domain. However, to demonstrate the robustness of the proposed method in dealing with irregular domain, the test is performed in a unit circular domain as shown in Fig. 6. The test is performed in a unit circular domain as shown in Fig. 6.

Fig. 7. Single vortex flow test: (a) CLSVOF on the Cartesian grid (1282 ) [24]; (b-c) results on the curvilinear grid: (b) VOSET (1002 ); (c) VOF-MYC (1282 ) [24].

810

Z. Cao et al. / Journal of Computational Physics 396 (2019) 799–818

Fig. 8. Pressure fields and spurious velocity of a static round droplet after one time step (10−4 s) with grid size h = 1/40 calculated by (a) VOSET on Cartesian coordinates, (b) VOSET on Curvilinear coordinates.

If the boundary mesh density is given, the non-orthogonal grid within the circular domain can be generated by solving the discrete Laplace equations iteratively. A time-reversed single vortex flow field same as those in [24] is imposed and given as follows

u = − sin(2π y ) sin2 (π x) cos(π t / T ), v = sin(2π x) sin2 (π y ) cos(π t / T )

(39)

where T refers to the time period and is equal to 8 s in this study. It can be noted that the flow direction is clockwise until t = T /2 and then is reversed to the anticlockwise. The fluid body will return to its initial position after a time period under the effect of time-reversed velocity field. To facilitate the comparison, the reconstructed interfaces at t = 4 s and t = 8 s are shown in Fig. 7, which includes the proposed VOSET method on the curvilinear grid, the CLSVOF method on the Cartesian grid and VOF with a mixed Youngs’ and center column (height function) (MYC) methods on the curvilinear grid presented by Wang et al. [24] in a unit square domain. The grid number in the unit circle domain is 100 × 100, the equivalent grid size is 8.86 × 10−3 ; while the total grid number of the unit square domain is 128 × 128, and the corresponding grid size is 7.81 × 10−3 . From Fig. 7a, it can be seen that the interface reconstructed by the VOSET method on the curvilinear grid at t = 4 s is very close to those by the CLSVOF method on the Cartesian grid. In contrast, more filament in the tail region appeared in the VOF-MFC method even if it has higher mesh density compared to the present VOSET method. Fig. 7b gives the reconstructed interfaces when the fluid body returns to its initial position at t = 8 s, it can be found that the VOSET result is much better than those by VOF-MFC method, demonstrating the high accuracy of the present method even if for irregular domain. 3.2. Static droplet The static droplet in equilibrium without gravity is frequently used to validate the numerical methods on calculating surface tension. In the case, an initial quiescent droplet with radius R = 0.25 is place at the center of a unit square domain fulfilled with the second fluid phase. The densities of the droplet and its surrounding fluid are 1.0 kg m−3 , the viscosities are set to 1.0 Pa s. The coefficient of surface tension is 0.1 N m−1 . Free-slip boundary conditions are imposed on all the walls. In the absence of gravitational and external forces, the round droplet and the surrounding fluid will remain stationary and the exact pressure jump across the interface can be given as follows

( p )exact = σ / R

(40)

Fig. 8 shows the pressure fields and spurious velocity with grid size h = 1/40 by VOSET based on Cartesian coordinates and curvilinear coordinates, respectively. It is found that the present results obtained by curvilinear coordinates are very

Z. Cao et al. / Journal of Computational Physics 396 (2019) 799–818

811

Fig. 9. Comparison on pressure curves on the slice y = 0.5 between Cartesian coordinates and Curvilinear coordinates.

Fig. 10. Physical domain and initial configuration (a) and grid structure (b) for dam break problem.

close to that by Cartesian coordinates. The exact pressure jump calculated by Eq. (40) is 0.40 Pa. As shown in Fig. 9, the pressure jump in the slice y = 0.5 computed by VOSET on curvilinear coordinates almost coincides with that by VOSET on Cartesian coordinates, both results show very small relative errors (less than 0.025%) compared to exact pressure jump,

812

Z. Cao et al. / Journal of Computational Physics 396 (2019) 799–818

Fig. 11. Interface locations at different time of dam break problem.

which is much more accurate than those results by using different volume tracking methods including geometric VOF in FLUENT, level set in COMSOL on unstructured grids in [22]. Thus it can be concluded that the present VOSET method for curvilinear coordinates can calculate the surface tension accurately and result in the correct pressure jump. 3.3. Dam break problem The dam break is another classical test for two-phase flow numerical methods. The structural diagram and associated curvilinear grid of dam break problem are shown in Fig. 10a and Fig. 10b, respectively. An initial liquid column with L × 2L is located at the bottom-left corner of the irregular region, where L equals to 0.146 m. The fluid physical parameters are described in the following: The water density ρl = 1.0 × 103 kg m−3 , viscosity μl = 0.5 Pa s, air density ρ g = 1.0 kg m−3 , viscosity μ g = 0.5 × 10−3 Pa s, gravity g = 9.8 m s−2 . Slip boundary condition is set on all boundaries. The simulation time is 0.9 s, and the time step is calculated by Eq. (38). Fig. 11 shows the evolution of the interfaces at different times by using the proposed method. Fig. 12 depicts the curve of the dimensionless liquid front location X ∗ = x/ L with the dimensionless time t ∗ = t (2g / L )1/2 . It can be found that the present VOSET results match pretty well with the experimental data [31] and the numerical results of the phase field method [33], further verifying the accuracy of the present VOSET method coupled with Navier–Stokes solver. 3.4. Single gas bubble rising problem The single bubble rising problem has become a benchmark case for testing incompressible interfacial flow codes. Hysing et al. [32] proposed the quantitative reference solutions for the benchmark quantities of a 2D single bubble rising problem while undergoing shape deformation in a liquid column. It is noted that the reference data is obtained in a  = [0, 1] × [0, 2] rectangular domain, however, we design a irregular domain to check the proposed method on the curvilinear grid. The gird coordinates are calculated by

xi , j = x0 y i , j = 2 y 0 + a sin(4π x0 + c )

  where x0 = (i − 1)/nx, y 0 = ( j − b)/(ny − b), a = 0.08 exp −0.02(x20 + y 20 ) , b = (10 + ny )/11 and c =

(41) 3π j . 2(ny +1)

Z. Cao et al. / Journal of Computational Physics 396 (2019) 799–818

813

Fig. 12. The predicted liquid front location for the dam break problem on irregular region by the proposed method.

Fig. 13. Grid structure for bubble rising problem and the bubble initial position.

Fig. 13a shows the mesh with grid number 50 × 100, and Fig. 13b shows the initial bubble of diameter d = 0.5 which is centered at (0.5,0.5). It is worth mentioning that contours of volume fraction in Fig. 13b are not uniform because all the volume fractions are stored at the center points of non-uniform grid cells. In the method, the interface in a cell is sharp and presented by a line segment for 2D problems. Two cases with low density and viscosity ratios and high density and viscosity ratios are considered, respectively. The corresponding fluid physical parameters are listed in Table 3. Slip boundary condition is imposed on the vertical walls, while no-slip boundary condition is set on the top and bottom walls. The simulation time is 3 s. For the sake of comparison, the benchmark quantities are defined as

 α y J dξ dη y c =   α J dξ dη  α v y J dξ dη  vc =   α J dξ dη

(42)

(43)

where y c is the vertical component of the bubble centroid, and v c denotes the mean bubble rising velocity component.

814

Z. Cao et al. / Journal of Computational Physics 396 (2019) 799–818

Table 3 Physical parameters defining the rising bubble test cases. Case

ρ1

ρ2

μ1

μ2

g

σ

ρ1 /ρ2

μ1 /μ2

1 2

1000 1000

100 1

10 10

1 0.1

0.98 0.98

24.5 1.96

10 1000

10 100

Fig. 14. Bubble shape and the surface tension vector as well as the signed distance function near the interface at t = 1.5 s, 3 s for low density ratio bubble rising problem.

Fig. 14 presents the results for case 1 at t = 1.5 s and t = 3.0 s. From Fig. 14a it can be observed that the bubble shape deforms as time goes on. The signed distance function as well as the surface tension vector near the interface are shown in Fig. 14b. It is clear that accurate dynamic distance function can be obtained by the current VOSET method.

Z. Cao et al. / Journal of Computational Physics 396 (2019) 799–818

815

Fig. 15. Bubble rising problem for case 1. (a) Bubble centroid. (b) Rise velocity (v c ). (c) Bubble locations compared with Hysing et al. [32] at time t = 3 s.

The quantitative results such as bubble centroid y c , rising velocity v c are given in Fig. 15a and Fig. 15b, respectively. The present results show good agreement with the reference solutions in [32]. Furthermore, as shown in Fig. 15c, the predicted bubble shape and location at the final time t = 3 s with a grid number 50 × 100 match very well with those reported by Hysing et al.. Fig. 16 depicts the variation curve of mass errors with time in the single gas bubble rising problem for case 1, in which the mass error can be defined as follow

E r = | M (t ) − M (0)|/ M (0) (44)  − 3 where M (t ) =  α (x, t ) J dξ dη . It can be found that the relative mass error is less than 10 during the whole simulation

process, which turns to increase mainly when the bubble undergoes significant deformation. The case 2 involves a rising bubble with high density and viscosity ratios. Final bubble shapes at t = 3 s with different grid resolutions (grid number = [104 , 2 × 104 , 8 × 104 ]) are plotted in Fig. 17, it can be seen that both the shape and the location of the bubble on the coarse grid are approximately the same with those on the fine grid. It is worth mentioning that Hysing et al. [32] have done comprehensive researches on this case by using three different codes (two including TP2D, FreeLIFE representing Eulerian level set finite-element codes and one- MooNMD representing an arbitrary LagrangianEulerian moving grid approach). Their results show that although all codes obtained a similar shape for the main bulk of the bubble, no agreement with respect to the thin filamentary regions has been reached. Fig. 18a depicts how the vertical position of the center of mass at different times converges from coarse grid to fine grid. It is surprising that the present results accord with the reference data perfectly. The possible reason is that the small thin filamentary regions do not influence the overall averaged quantities to a significant degree. Fig. 18b depicts the curves of the mean rising velocity of the

816

Z. Cao et al. / Journal of Computational Physics 396 (2019) 799–818

Fig. 16. Variation curve of mass errors with time in the single gas bubble rising problem for case 1. Mass conservation error E r = | M (t ) − M (0)|/ M (0). Here, M (t ) =  α (x, t ) J dξ dη .

Fig. 17. Bubble shape at the final time t = 3 s with different grid resolutions for case 2.

bubble at different times from a coarse mesh to fine mesh. It can be seen that the first maximum rise velocity is 0.252 m/s occurring at t = 0.725 s. The present results show good agreement with the reference solution that the first maximum rise velocity has a magnitude of 0.25 ± 0.01 m/s and occurs around t = 0.73 ± 0.02 s. 4. Conclusions In this paper, we develop a coupled volume-of-fluid and level set method on general curvilinear grids for interfacial flow simulations. The proposed method not only performs good characteristics of mass conservation but also predicts the surface tension with good accuracy. The concluding remarks of this study are summarized as follows.

Z. Cao et al. / Journal of Computational Physics 396 (2019) 799–818

817

Fig. 18. Bubble rising problem for case 2. (a) Bubble centroid. (b) Rise velocity(v c ).

(1) An iterative geometric operation is extended to physical domain and computational domain for calculating the level set function φ near interfaces, which can be applied to compute the accurate curvature κ and smooth the discontinuous physical quantities near interfaces. A direct construction of the distance function rather than contravariant transformation as in [24] is proposed for computing the level set function on the physical domain. (2) The level set based CSF model on curvilinear grids is derived. It is of critical importance to accurately model the surface tension especially for problems in which surface tension is the dominant force. The level set based CSF model as well as the proposed volume tracking method is incorporated into an incompressible Navier–Stokes solver. (3) To verify the accuracy and robustness of the proposed method to deal with complex domains, several benchmark problems including time-reversed single vortex flow, static droplet, dam break and single bubble rising problems are tested on irregular domains. The results demonstrate good agreement with previous experimental data or numerical results, indicating that the present method can accurately simulate incompressible two-phase flows. Acknowledgements The study is supported by the National Natural Science Foundation of China (Nos. 51636006, 51776019), the International Cooperative Project (No. 51611130060), Jointly Projects of Beijing Natural Science Foundation and Beijing Municipal Education Commission (No. KZ201810017023), the Project of Construction of Innovative Teams and Teacher Career Devel-

818

Z. Cao et al. / Journal of Computational Physics 396 (2019) 799–818

opment for Universities and Colleges under Beijing Municipality (No. IDHT20170507), the Program of Great Wall Scholar (No. CIT&TCD20180313). References [1] G. Kozyrakis, A. Delis, N. Kampanis, A finite difference solver for incompressible Navier–Stokes flows in complex domains, Appl. Numer. Math. 115 (2017) 275–298. [2] J.J. Quirk, An alternative to unstructured grids for computing gas dynamic flows around arbitrarily complex two-dimensional bodies, Comput. Fluids 23 (1) (1994) 125–142. [3] B. Nichols, C. Hirt, R. Hotchkiss, Sola-vof: A Solution Algorithm for Transient Fluid Flow with Multiple Free Boundaries, Tech. rep., Los Alamos Scientific Lab, 1980. [4] D. Kothe, W. Rider, S. Mosso, J. Brock, J. Hochstein, Volume tracking of interfaces having surface tension in two and three dimensions, in: 34th Aerospace Sciences Meeting and Exhibit, 1996, p. 859. [5] A. Nikseresht, M. Alishahi, H. Emdad, Complete flow field computation around an acv (air-cushion vehicle) using 3d vof with lagrangian propagation in computational domain, Comput. Struct. 86 (7–8) (2008) 627–641. [6] J. López, J. Hernández, Analytical and geometrical tools for 3d volume of fluid methods in general grids, J. Comput. Phys. 227 (12) (2008) 5939–5948. [7] A. Sarthou, S. Vincent, J.-P. Caltagirone, A second-order curvilinear to cartesian transformation of immersed interfaces and boundaries. Application to fictitious domains and multiphase flows, Comput. Fluids 46 (1) (2011) 422–428. [8] V.-T. Nguyen, W.-G. Park, A volume-of-fluid (vof) interface-sharpening method for two-phase incompressible flows, Comput. Fluids 152 (2017) 104–119. [9] V.-T. Nguyen, W.-G. Park, Enhancement of Navier–Stokes solver based on an improved volume-of-fluid method for complex interfacial-flow simulations, Appl. Ocean Res. 72 (2018) 92–109. [10] S. Osher, J.A. Sethian, Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations, J. Comput. Phys. 79 (1) (1988) 12–49. [11] W. Yue, C.-L. Lin, V.C. Patel, Numerical simulation of unsteady multidimensional free surface motions by level set method, Int. J. Numer. Methods Fluids 42 (8) (2003) 853–884. [12] Z. Wang, Q. Zou, D. Reeve, Simulation of spilling breaking waves using a two phase flow cfd model, Comput. Fluids 38 (10) (2009) 1995–2005. [13] 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 (2) (2000) 301–337. [14] G. Son, N. Hur, A coupled level set and volume-of-fluid method for the buoyancy-driven motion of fluid particles, Numer. Heat Transf., Part B, Fundam. 42 (6) (2002) 523–542. [15] M. Haghshenas, J.A. Wilson, R. Kumar, Algebraic coupled level set-volume of fluid method for surface tension dominant two-phase flows, Int. J. Multiph. Flow 90 (2017) 13–28. [16] M. Skarysz, A. Garmory, M. Dianat, An iterative interface reconstruction method for PLIC in general convex grids as part of a coupled level set volume of fluid solver, J. Comput. Phys. 368 (2018) 254–276. [17] X. Yang, A.J. James, J. Lowengrub, X. Zheng, V. Cristini, An adaptive coupled level-set/volume-of-fluid interface capturing method for unstructured triangular grids, J. Comput. Phys. 217 (2) (2006) 364–394. [18] A. Ferrari, M. Magnini, J.R. Thome, A flexible coupled level set and volume of fluid (FLEXCLV) method to simulate microscale two-phase flow in non-uniform and unstructured meshes, Int. J. Multiph. Flow 91 (2017) 276–295. [19] D. Sun, W. Tao, A coupled volume-of-fluid and level set (VOSET) method for computing incompressible two-phase flows, Int. J. Heat Mass Transf. 53 (4) (2010) 645–655. [20] K. Ling, Z.-H. Li, D.-L. Sun, Y.-L. He, W.-Q. Tao, A three-dimensional volume of fluid & level set (VOSET) method for incompressible two-phase flow, Comput. Fluids 118 (2015) 293–304. [21] Z. Cao, D. Sun, B. Yu, J. Wei, A coupled volume-of-fluid and level set (VOSET) method based on remapping algorithm for unstructured triangular grids, Int. J. Heat Mass Transf. 111 (2017) 232–245. [22] Z. Cao, D. Sun, J. Wei, B. Yu, A coupled volume-of-fluid and level set method based on multi-dimensional advection for unstructured triangular meshes, Chem. Eng. Sci. 176 (2018) 560–579. [23] Z. Cao, D. Sun, B. Yu, J. Wei, A coupled volume of fluid and level set method based on analytic PLIC for unstructured quadrilateral grids, Numer. Heat Transf., Part B, Fundam. 73 (4) (2018) 189–205. [24] Z. Wang, J. Yang, F. Stern, A new volume-of-fluid method with a constructed distance function on general structured grids, J. Comput. Phys. 231 (9) (2012) 3703–3722. [25] Z. Wang, J. Yang, B. Koo, F. Stern, A coupled level set and volume-of-fluid method for sharp interface simulation of plunging breaking waves, Int. J. Multiph. Flow 35 (3) (2009) 227–246. [26] J. Brackbill, D.B. Kothe, C. Zemach, A continuum method for modeling surface tension, J. Comput. Phys. 100 (2) (1992) 335–354. [27] K. Yokoi, A practical numerical framework for free surface flows based on CLSVOF method, multi-moment methods and density-scaled CSF model: numerical simulations of droplet splashing, J. Comput. Phys. 232 (1) (2013) 252–271. [28] W. Shyy, Elements of pressure-based computational algorithms for complex fluid flow and heat transfer, in: Advances in Heat Transfer, vol. 24, Elsevier, 1994, pp. 191–275. [29] J.B. Bell, P. Colella, H.M. Glaz, A second-order projection method for the incompressible Navier–Stokes equations, J. Comput. Phys. 85 (2) (1989) 257–283. [30] T. Wang, H. Li, Y. Feng, D. Shi, A coupled volume-of-fluid and level set (VOSET) method on dynamically adaptive quadtree grids, Int. J. Heat Mass Transf. 67 (1) (2013) 70–73. [31] J.C. Martin, W.J. Moyce, W.G. Penney, A. Price, C. Thornhill, Part iv. An experimental study of the collapse of liquid columns on a rigid horizontal plane, Philos. Trans. R. Soc. Lond. A 244 (882) (1952) 312–324. [32] S.-R. Hysing, S. Turek, D. Kuzmin, N. Parolini, E. Burman, S. Ganesan, L. Tobiska, Quantitative benchmark computations of two-dimensional bubble dynamics, Int. J. Numer. Methods Fluids 60 (11) (2009) 1259–1288. [33] P.-H. Chiu, Y.-T. Lin, A conservative phase field method for solving incompressible two-phase flows, J. Comput. Phys. 230 (1) (2011) 185–204.