A coupled volume-of-fluid and level set (VOSET) method based on remapping algorithm for unstructured triangular grids

International Journal of Heat and Mass Transfer 111 (2017) 232–245

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

A coupled volume-of-fluid and level set (VOSET) method based on remapping algorithm for unstructured triangular grids Zhizhu Cao a, Dongliang Sun b,⇑, Bo Yu b,⇑, Jinjia Wei a a b

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

a r t i c l e

i n f o

Article history: Received 11 October 2016 Received in revised form 7 March 2017 Accepted 27 March 2017

Keywords: Two-phase flow VOSET method Remapping algorithm Unstructured triangular grid

a b s t r a c t We propose a coupled volume-of-fluid and level set (VOSET) method for interfacial flow simulations on unstructured triangular grids. In this method, the volume fraction advection is performed using a Lagrangian–Eulerian remapping algorithm, and the level set function is calculated by a simple iterative geometric operation. The present VOSET method can not only satisfy mass conservation but also predict the surface tension with good accuracy, thus combining the advantages and overcoming the disadvantages of volume-of-fluid (VOF) and level set (LS) methods. Finally, the present method is verified by well known Zalesak’s slotted-disk revolution, single vortex flow, dam break and single bubble rising problems. The results illustrate that the present method can accurately simulate incompressible two-phase flows for unstructured triangular grids. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction Capturing the moving interface is a critical issue in the simulation of two-phase flows. In the past few decades, there have been a variety of numerical methods proposed for solving interfacial flow problems, such as volume-of-fluid (VOF) [1–6], level set (LS) [7– 10], boundary integral [11] and front tracking methods [12,13]. The VOF method was originally developed by Hirt and Nichols [1] to track interfacial dynamic behaviors. In this method, a volume fraction is described with a and denotes whether a grid cell is occupied by the target or non-target phase. If a is 1, then the grid cell is filled with target fluid, if a is 0, it is filled with non-target fluid, and when a is between 0 and 1, the interface is located in the cell. The volume fraction advection is mainly performed by two algorithms. One is based on Eulerian framework [2–4], in which a is evolved by solving a standard advection equation. The other is under Lagrangian-Eulerian framework, where a is calculated by employing a remapping algorithm [5,6]. The VOF method has good characteristics of mass conservation. However, it is hard to calculate the surface tension force accurately, because a changes abruptly across the interface and is subject to a step(discontinuous) function. The accuracy of surface tension force solved by the

⇑ Corresponding authors. E-mail addresses: [email protected] (D. Sun), [email protected] (B. Yu). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2017.03.096 0017-9310/Ó 2017 Elsevier Ltd. All rights reserved.

CSF model will become worse with the grid refinement as the error norm for curvature grows linearly when increasing the mesh resolution [14]. The weaknesses of the VOF method in the surface force calculation have also been demonstrated by Gerlach et al. [15]. The LS method has been widely used since its development in 1988 by Osher and Sethian [7]. In this method, a smooth function !

/ðx ; tÞ, called the level set function, is used to represent the inter!

face as the set where /ðx ; tÞ ¼ 0. The level set function / is defined to be a signed distance function [8] !

!

!

!

j/ðx ; tÞj ¼ dðx ; tÞ ¼ minxC 2C ðj x x C jÞ

ð1Þ

!

where C is the interface, /ðx ; tÞ < 0 on one side of the interface and !

/ðx ; tÞ > 0 on the other. The evolution of / in the velocity field is given by [16]

@/ ! þ u r/ ¼ 0 @t

ð2Þ

!

where u is the velocity on the interface. However, the shape of level set obtained by Eq. (2) may get severely distorted and the level set may vanish over several time steps due to numerical errors, a reinitialization process is needed to retain / as a signed distance function after every time step. It can be achieved by solving the following transient problem to the steady state.




/s þ Sð/0 Þð1  jr/jÞ ¼ 0 /ðx! ; 0Þ ¼ /0

ð3Þ

233

Z. Cao et al. / International Journal of Heat and Mass Transfer 111 (2017) 232–245

where s is an artificial time, not equivalent to the actual time t, /s is the time derivative, S is the sign function. The level set function / will satisfy jr/j ¼ 1 as time progresses. The advantage of the LS method is that topological changes of the interfaces such as breaking and merging can be performed automatically, because the zero-level set always gives the location of the propagating interface, and there is no need to reconstruct the interface explicitly as with the VOF method. Moreover, since / is a smooth function, the surface tension force as well as the curvature could be calculated with high order of accuracy [8] and the discontinuous physical parameters such as density and viscosity calculated near interfaces could be perfectly smooth as well [14] by using /. However, for the lack of geometric conservative law, there is no guarantee that the employed reinitialization process can preserve flow volume in time. The phenomenon of fluid mass loss/gain may happen in each time step. As time advances, even a negligibly small error can be accumulated to a sufficiently large magnitude [17]. Therefore, the LS method produces more numerical errors than the VOF method, leading to loss /gain of mass, especially when the interfaces undergo severe stretching and tearing. It is worth noting that some conservative LS methods have been reported to overcome the gain/loss of mass. Olsson et al. [8] proposed a conservative LS method in which a Heaviside function was used as the level set function. Conservative second order TVD methods are adopted to solve the level set function. Moreover, an intermediate step that maintains the resolution of contact discontinuities after each advection step is added to make sure that the smooth profile of / is preserved and the interface keeps its thickness and shape in the reinitialization process. Yu et al. [9] developed an improved interface preserving LS method. A mass correction term is added to the reinitialization equation to conserve mass bounded by the interface. Besides, Moghadam et al. [10] presented a compact conservative LS method, in which the hyperbolic tangent function used as LS function is conservatively transported by a high-order compact finite difference scheme. To combine the advantages of VOF and LS methods, Sussman and Puckett [18] proposed a coupled level set and volume-offluid (CLSVOF) method for computing two-phase flows. Later, the CLSVOF method was extended to the adaptive triangular grids for resolving complex interface changes and interfaces of high curvature efficiently and accurately [19]. Chakraborty et al. [20] verified the accuracy of the CLSVOF method, and adopted the method to study the dynamics of the formation of drops into air. Recently, a coupled volume-of-fluid and level set (VOSET) method was developed by the present author Sun and Tao [14]. In the VOSET method, only the volume fraction advection equation needs to be solved and the LS function is calculated 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 way. In the CLSVOF method, two differential equations,i.e., the level set advection equation (Eq. (2)) and the reinitialization equation (Eq. (3)), are needed to solve in a iterative way for calculating the level set function. These iterative processes have to be done in the whole region and thus it’s time-consuming. Moreover, it seems to be complex because the accuracy of / is dependent on numerical schemes adopted in the discretization process. Therefore, the VOSET method is simpler compared with the CLSVOF method. Afterwards, the VOSET method was extended to the three-dimensional structured grids by Ling et al. [21] and the dynamically adaptive quadtree/octree grids by Wang et al. [22,23]. Ansari et al. [24] applied the VOSET method to study the interface topological changes of gas–liquid two-phase flows. The corresponding simulation results of bubble rising agree very well with the experimental data. It can be seen that in addition to achieving mass conservation, the method computes the surface tension with a good accuracy, thus verifying the

superiority of the method. The VOSET method was also adopted to simulate the film boiling problems by Guo et al. [25], the nucleate boiling problems by Ling et al. [26,27], the dynamics of a ferrofluid droplet under a uniform magnetic field by Shi et al. [28]. To the authors’ knowledge, all of the above researches about the VOSET method were based on Cartesian structured grids. In general, it is hard to fit in with complex irregular domains. In order to overcome this problem, the present work is aimed at extending the VOSET method to unstructured triangular grids for simulating the two-phase flows in complex irregular domains. The paper is organized as follows. The solution procedure of the VOSET method on unstructured triangular grids is first introduced. Then the present method is verified by well known Zalesak’s slotted-disk revolution, single vortex flow, dam break and single bubble rising problems. Finally, some conclusions are drawn. 2. VOSET on unstructured triangular grids 2.1. Iterative geometric operation In the VOSET method, the iterative geometric operation is used to reconstruct interface and calculate the LS function /, i.e. the signed distance function. The detailed procedure of the operation is presented as follows. Step 1: Reconstruct interface with the piecewise linear interface construction (PLIC) method. The interface reconstruction is based on the idea that a normal vector together with the volume fraction function a determines a unique linear interface cutting the cell. The normal vector to an interface is estimated by !

n ¼ ra=jraj

ð4Þ

Unlike the structured grids, the least-square method is used to calculate ra in Eq. (4) based on a stencil which identifies neighboring center points of the cell [29]. !

Based on values of a and n , an efficient analytic algorithm [30] is adopted to quickly determine the interface position on unstructured triangular grids. There are two possible cases for the interface position in the triangular cell, as shown in Fig. 1. Firstly, the three vertices of the cell are sorted, where vertex A is in the fluid and vertex C is outside the fluid. For vertex B, there are two possibilities, inside or outside fluid. The sorting procedure is implemented in the following manner. The three vertices of the cell are assumed as p1 ðx1 ; y1 Þ; p2 ðx2 ; y2 Þ; p3 ðx3 ; y3 Þ. If p1 satisfies the relation: !

!

!

!

ðp1 p2  n Þ  ðp1 p3  n Þ > 0

ð5Þ ! p1 p2

!

! p1 p3

!

then p1 is vertex A. In this case, if n <  n , then p2 is vertex B and p3 is vertex C. Otherwise, p3 is vertex B and p2 is vertex C. Thus, the three vertices A, B and C can be sorted in a specific order. Secondly, the coordinates of the two endpoints of the interface segment (E and F) are determined. As shown in Fig. 1, the endpoint F is located at the cell edge AC, and the endpoint E is either on the cell edge AB or on the edge BC. For judging the location of the end!

point E, a line BD perpendicular to the normal vector n is drawn. Then, the judgment rule can be expressed as

8 < E 2 AB if ða ¼ SSAEF Þ 6 ða ¼ SSABD Þ ABC ABC     : E 2 BC if a ¼ SAEF > a ¼ SABD SABC SABC

ð6Þ

For the case of a 6 a , the coordinates of E and F can be easily obtained by

234

Z. Cao et al. / International Journal of Heat and Mass Transfer 111 (2017) 232–245

Interface

(a) Case of α ≤ α *

(b) Case of α > α *

Fig. 1. Two cases for the interface position in the triangular cell.

Flagged region rffiffiffiffiffi

pffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! a ! AB and X F ¼ X A þ a  a AC a

XE ¼ XA þ

ð7Þ

For the case of a > a , the coordinates of E and F can also be calculated by

rffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! 1a ! CB and X F ¼ X C þ ð1  aÞ  ð1  a Þ CA XE ¼ XC þ  1a

ð8Þ

Step 2: Set the initial value for the signed distance function / in the whole computational domain.


/0i ¼

if ai P 0:5 if ai < 0:5

M M

ð9Þ

where M is the max geometrical size of computational domain and the subscript i denotes the grid cell index. This step can guarantee that the signed distance function calculated in the following step is included between M and M, i.e. M < / < M. Step 3: Flag the cells near the interfaces. As shown in Fig. 2, the width of flagged region is 6D, where D denotes the grid size. The flagged region is wide enough to compute the surface tension force and smooth the discontinuous physical parameters near the interfaces. This step can avoid calculating the signed distance function in the whole computational domain for the sake of computational saving. Step 4: Calculate the signed distance function / in the flagged region. In order to obtain /, the shortest distance to the interface should be first computed. As shown in Fig. 3, the shortest distance di at cell (i) can be obtained by comparing all of the minimum distances from cell central point P(i) to any interface within the 6D  6D stencil around the cell (i). Fig. 4 shows the method of calculating the minimum distance from cell central 

point P(i) to an interface EF in a cell. When h1 > 90 , PE is the 

minimum distance; when h2 > 90 , PF is the minimum dis



tance; when h1 < 90 and h2 < 90 , the vertical line PM is the minimum distance. Then, the signed distance function can be obtained by

8 > < di /i ¼ 0 > : di

ai > 0:5 ai ¼ 0:5 ai < 0:5

ð10Þ

Based on the above signed distance function, the more accurate normal vector to the interface can be calculated by !

n ¼ r/=jr/j

ð11Þ

Here, the least-square fit is used to obtain r/ in Eq. (11) based on a stencil which identifies neighboring center points of the cell (i) (see Fig. 8). Specifically, the value of / in the jth neighboring cell can be expressed by a Taylor series expansion for a linear approximation at the ith cell. !

!

/j ¼ /i þ ðx j  x i Þ  ðr/Þi j ¼ 1; 2; . . . ; n

Fig. 2. The flagged region near interfaces.

ð12Þ

The above equations can be also expressed in the following matrix form

1 1 0 /1  /i x1  xi y 1  y i C C B B Bx  x y  y C B/  / C ! i B 2 B 2 iC 2 iC @/=@x C C B B ¼B C B . .. .. C C @/=@y C B . B C B . . . C |fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl}i B A A @ @ ðr/Þi /n  /i xn  xi yn  yi |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl} 0

L

ð13Þ

Y

The over-determined system of linear equations can be solved by the least-square fit [31] and the gradient can be obtained by 1

ðr/Þi ¼ ðLT LÞ LT Y

ð14Þ

Then, return to Step 1 and reconstruct interface again by using the accurate normal vector. Repeat the geometric operation until the iteration times is equal to the pre-specified times N. By the iteration, we can obtain more accurate signed distance function, and reconstruct more accurate and smoother interface simultaneously. Since the changes of / are small during each time step, the signed distance function obtained after several iterations is very close to the convergence solution, there is no need to set the iteration times N as a great number to save computational cost. In numerical experiments for the unstructured triangular grids, we have found that if N is set as 2, a good balance is achieved between the precision requirements and low computational cost. Therefore, the iteration times N is recommended to be set as 2. Fig. 5 shows the signed distance function near the interfaces with different shapes, which is calculated by the proposed iterative geometric operation. In Fig. 5 the red1 lines denote the interfaces and the blue lines represent the contour lines of the signed distance function. From the figure, it can be found that the accurate signed distance function can be calculated by the proposed iterative geometrical operation for unstructured triangular grids. Fig. 6 compares the traditional PLIC method with the PLIC method based on the proposed iterative geometric operation. It is obvious that the interface reconstructed by our proposed method is more accurate and smoother than that obtained by the traditional method. The underlying reason is that the more accurate normal vector to the interface can be calculated by the proposed method. 1 For interpretation of color in Fig. 5, the reader is referred to the web version of this article.

Z. Cao et al. / International Journal of Heat and Mass Transfer 111 (2017) 232–245

Interface

235

Shortest distance

P(i)

Minimum distances from P(i) to any interface in cells within 6 ×6 stencil Fig. 3. The shortest distance determined by comparing all of the minimum distances.

(a) θ1>90 °

(b) θ2>90°

(c) θ1<90°and θ2<90°

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

Interface

Interface

Signed distance function

Signed distance function

(a) Circle interface

(b) Zalesak's slotted-disk interface

Fig. 5. The signed distance function near the interfaces with different shapes calculated by the iterative geometrical operation.

2.2. Lagrangian-Eulerian remapping algorithm In the Lagrangian-Eulerian remapping algorithm, the projection of the volume fraction is equivalent to solving the volume fraction advection equation given below.

@a ! þ u ra ¼ 0 @t

ð15Þ

The Lagrangian-Eulerian remapping algorithm [19] is briefly reviewed as follows.

Step 1: Obtain the Lagrangian grids. Each grid cell can be considered as a material element. Then each material element transports with the flow. Based on this idea, we can obtain the coordinates of vertices of each Lagrangian grid cell by using a second order Runge–Kutta scheme (see Fig. 7). Step 2: Reconstruct interface on the Lagrangian grids. It is firstly considered that the volume fraction in each Lagrangian grid cell at the n + 1 time level are equal to that in the corresponding Eulerian grid cell at the n time level.

236

Z. Cao et al. / International Journal of Heat and Mass Transfer 111 (2017) 232–245

0.7

0.7

Initial circle Reconstructed interface

0.69

0.68

0.68

y

y

0.69

0.67

0.67

0.66

0.66

0.65

0.65

0.64 0.35

Initial circle Reconstructed interface

0.36

0.37

0.38

x

0.39

0.4

0.41

0.64 0.35

(a) Traditional PLIC method

0.36

0.37

0.38

x

0.39

0.4

0.41

(b) PLIC method based on the iterative geometric operation

Fig. 6. Piecewise linear interface reconstruction on unstructured triangular grids. (Zoom of a specific region of a circle.)

geometric operation to reconstruct interface and calculate the LS function on the Lagrangian grids. Step 3: Remap the Lagrangian volume fraction back to the Eulerian grids. The new volume fraction on the original Eulerian grids is obtained by remapping the Lagrangian volume fraction back to the Eulerian grids. The remapping procedure is introduced in the following. According to the mass conversation law, the density of target fluid in each Lagrangian grid cell at the n + 1 time level can be calculated by

ðqnþ1 ÞL ¼ i

ðqni ÞE ðAni ÞE ðani ÞE ðAnþ1 ÞL ðanþ1 ÞL i i

¼

ðqni ÞE ðAni ÞE ðAnþ1 ÞL i

ð17Þ

where A denotes the grid cell area. Through the polygon-polygon clipping method [32], the target fluid in a Lagrangian grid cell (i) can be clipped and allocated to its surrounding Eulerian grid cells (j), and the clipped area can be ÞL and denoted as ðSi!j ÞL . Due to the difference between ðqnþ1 i ðqni ÞE , the clipped area needs to be transferred from Lagrangian to Eulerian grids using

Fig. 7. Schematic of the Lagrangian–Eulerian remapping.

ðSi!j ÞE ¼ ðSi!j ÞL

ðqnþ1 ÞL i ðqni ÞE

ð18Þ

Finally, the new volume fraction on the Eulerian grid cell (j) at the n + 1 time level can be obtained by summing up all of the clipped areas allocated to it. 2.3. Governing equations and their solutions 2.3.1. Governing equations The governing equations for transient, incompressible, Newtonian, two-phase flows include the continuum and momentum equations, which can be expressed as Fig. 8. A stencil used to determine the interface thickness e corresponding to an interface cell (i). nþ1 ÞL i

ða

n i ÞE

¼ ða

ð16Þ

where the subscripts L and E denote the Lagrangian and Eulerian grids, respectively. Then, we can adopt the proposed iterative

!

ru ¼0

ð19Þ

!
  ! ! ! T ! ! @u ! 1 þ ðu rÞ u ¼ rp þ r  le ð/Þ r u þðr u Þ þqe ð/Þ g þF sv @t qe ð/Þ

ð20Þ

237

Z. Cao et al. / International Journal of Heat and Mass Transfer 111 (2017) 232–245

(1) Physical parameters The density and viscosity in Eq. (20) can be calculated by

qe ð/Þ ¼ Hð/Þqg þ ð1  Hð/ÞÞql

ð21Þ

le ð/Þ ¼ Hð/Þlg þ ð1  Hð/ÞÞll

ð22Þ

where the subscripts g and l denote the gas and liquid phases, respectively. The smoothed Heaviside function is written as

8 0 if / < e > <

p/ / 1 1 Hð/Þ ¼ 2 1 þ e þ p sin e if j/j 6 e > : 1 if / > e

ð23Þ

(2) Surface tension force The surface tension force in Eq. (20) is computed by the LS function based on the density-scaled CSF model [33] and can be expressed as !

[36] is used to couple the velocity and pressure (Eqs. (19) and (20)). The deferred-correction central-difference scheme, which has second-order accuracy, is applied for convection term in the momentum equation (Eq. (20)). The second-order TVD Runge– Kutta scheme treats the time marching. The time step Dt at time tn is determined from the CFL condition, surface tension, gravity and viscous terms [31,34], and it can be calculated by

Dt ¼ kmin

h ; jjun jj



h jjgjj

1=2 ;

qn ðhÞ2 ; ln



ql þ qg 4pr

!

1=2 ðhÞ

3=2

ð28Þ

where the constant k is a value between zero and unity, h is the dimensional grid size, k = 0.05 and the maximum time step size is set to 104 for current unstructured VOSET method. 3. Results and discussion

where r is surface tension coefficient and jð/Þ ¼ r  ðr/=jr/jÞ. The smoothed delta function is obtained from

In this section, the present unstructured VOSET method is tested and verified with four numerical experiments, including the well known Zalesak’s slotted-disk revolution [37], single vortex flow [38], dam break and single bubble rising problems.

dscaling ð/Þ ¼ 2He ð/Þde ð/Þ e

3.1. Zalesak’s slotted-disk revolution problem

F sv ¼ rjð/Þdscaling ð/Þr/ e

de ð/Þ ¼

dHð/Þ ¼ d/

(

ð24Þ

ð25Þ



 1 1 þ cos pe/ if j/j 6 e 2e 0

if j/j > e

ð26Þ

and dscaling ð/Þ satisfies the following relation: e

Z

e

e

dscaling ð/Þd/ ¼ 1 e

ð27Þ

It is worth mentioning that the adopted model is equivalent to the volume fraction based density-scaled CSF model [34,35]. If the density scaling is used, the force distribution as well as the peak of force distribution is shifted to higher density region and improves the stability of surface force computations. (3) Interface thickness e The interface thickness e plays a key role to smooth the discontinuous physical parameters (Eqs. (21) and (22)) and calculate the accurate surface tension force (Eq. (24)). The value of e should be chosen with caution. A smaller e will reduce the smearing of density, viscosity near the interface, however, since the surface tension force acts only in the transition region, which is of finite extent 2⁄e around the interface, a too small e compared to the grid size might lead to zero surface tension force existing in some interface cells when center points of these cells are outside of the transition region, and thus can have a negative effect upon the accuracy of calculation of surface tension force, whereas a too large e can violate the locality of the surface tension force [15]. There is a numerical restriction on how small we can choose e to obtain accurate interface normal and curvature. Besides, the interface thickness e should be set as different values at different locations because the unstructured triangular grids are usually non-uniform. The value of e corresponding to an interface cell can be determined on the basis of a stencil as shown in Fig. 8. An interface cell (i) is surrounded by the adjacent cells (j = 1–12). e can be obtained by comparing all of the maximum distances dmax from the adjacent cell central points to the interface cell vertex. Then the longest distance is chosen as e. 2.3.2. Discretization and solution of governing equations The governing equations are discretized on a collocated grid system by the finite volume method (FVM). The PISO algorithm

The slotted-disk revolution problem presented by Zalesak [37] is frequently employed for verifications of the interface-tracking methods. In this paper, the simulation conditions are the same as those in Ningegowda’s simulation [39] to provide better contrast between the CLSVOF method and the present VOSET method. Specifically, as shown in Fig. 9, in a 1  1 computational domain, a slotted disk with radius of 0.15 and slot width of 0.05 and length of 0.25 is initially located at (0.5, 0.75), and is revolved around the domain center (0.5, 0.5) along the counterclockwise direction. A spatially varying velocity field is given as the following

u ¼ ðy þ 0:5Þ=2;

v ¼ ðx  0:5Þ=2

ð29Þ

After a time period of 4p, the slotted disk completes one full rotation. We can compare the shape of the slotted disk with its initial geometry for verifying the accuracy of the VOSET method. As shown in Fig. 9, the shape obtained from the VOSET method on unstructured triangular grids is in good agreement with the CLSVOF result on structured grids [39]. The sharp corners are smoothed out due to the linear geometrical reconstruction of the interface. The difference between the numerical reconstruction and the initial geometry is visible only in the small regions near the sharp corners of the slotted disk. To be specific, the VOSET result for unstructured grids covers most of initial interfaces, and looks more accurate than the CLSVOF result at the right width of the rectangle. However, CLSVOF shows a more accurate result for corner regions of the slotted disk. Ling et al. [21] compared the VOSET for structured grids with the CLSVOF in 3D deformation test for an initial sphere, and their results showed that CLSVOF gave a better result. Of course, it should be admitted that CLSVOF is one of state-of-the-art methods and may be slightly better than VOSET, but VOSET is simpler than CLSVOF in which level set function and VOF function both need to be advected. 3.2. Single vortex flow problem To validate the present VOSET method, a single vortex flow problem is conducted. For the sake of contrast, we use the same test condition as the Ref. [14]. An initially circular fluid body with radius of 0.15 is located at (0.5, 0.75) in a 1  1 computational domain. The circular fluid body is evolved in a single vortex flow field, which can be expressed as

238

Z. Cao et al. / International Journal of Heat and Mass Transfer 111 (2017) 232–245

Initial geometry

One full rotation CCW

(a) CLSVOF method on structured grids [39] (grid size=1/100)

(b) VOSET method on unstructured triangular grids (grid size≈1/100)

Fig. 9. Rotation of Zalesak’s slotted disk simulated by different methods.

2

u ¼ sinð2pyÞsin ðpxÞ;

v ¼ sinð2pxÞsin2 ðpyÞ

ð30Þ

The single vortex flow field was first introduced by Bell et al. [38], and widely adopted to assess the integrity and capability of an interface-tracking method. The flow field contains a single vortex centered in the computational domain. The vortex spins fluid elements and stretches them into some thin filaments that spiral toward the vortex center. Thus each fluid element can undergo very significant deformation. As shown in Fig. 10, the evolution of the initially circular fluid body in the single vortex flow field was simulated by the VOSET method. It can be found that the fluid body is spun and stretched into a thinning filament, which agrees very well with those results obtained by other advanced methods [2,6,19]. In theory, the filament spirals toward the vortex center and gets thinner and thinner continuously. However, in the simulation process, the thin filament breaks into many small pieces beginning at its tail after long time evolution. The reason is that the thickness of the filament reaches the size of grid. To further verify the accuracy of the VOSET method, a timereversed single vortex flow field is also introduced and can be obtained by multiplying Eq. (30) by cos(pt/T) as follows. 2

u ¼ sinð2pyÞsin ðpxÞ cosðpt=TÞ; 2

¼ sinð2pxÞsin ðpyÞ cosðpt=TÞ

v ð31Þ

where T denotes the time period and is set as 8 s. The flow direction is initially clockwise and then is reversed to anticlockwise at time T/2 such that any fluid body returns to its initial position after a time period. Fig. 11 shows the interface locations, which are simulated by three different methods, after a time period in the timereversed single vortex flow field. As shown in this figure, the present VOSET method is more accurate than the VOF and LS methods. Fig. 11(a) and (b) show that there is a small piece separated from the main body in the simulation results of VOF method. Moreover, the VOSET method shows obvious advantage compared with the LS method even when the grid number of the later is ten times more than that of the former (see Fig. 11(c) and (d)), proving that the present method is of high accuracy for interface capturing. 3.3. Dam break problem The dam break problem is conducted to validate the present VOSET method by comparing the numerical results with the exper-

imental data [40]. Fig. 12 shows the structural diagram of dam break problem. An initial liquid column with width of L and height of 2L is located at the bottom-left corner of a 4L  4L square domain, where L equals to 0.146 m. The liquid density ql = 1.0  103kgm3, viscosity ll = 0.5 Pas, background gas density qg = 1.0 kgm3, viscosity lg = 0.5  103Pas, gravity g = 9.8 ms2 and surface tension coefficient r = 0.0755 Nm1. The dimensionless grid size is h = 4x/L = 1/50 and time step is determined by the time step criterion (Eq. (28)) with an upper limit of 4t⁄ = 1.16  103(Dt ¼ Dtð2g=LÞ1=2 ). Fig. 13 shows the dimensionless liquid front location X⁄ = x/L plotted along with the dimensionless time t⁄ = t(2g/L)1/2. By the comparison, it can be found that the present VOSET results agree very well with the experimental data [40] and the numerical results [14,41], verifying the accuracy of the present VOSET method. 3.4. Single gas bubble rising problem It is well known that transient behaviors of a rising gas bubble in liquid are categorized using four dimensionless numbers [42], including Morton number (M), Eotvos number (Eo), viscosity ratio (kl) and density ratio (kq).

M ¼ g l4l =ql r3 2

Eo ¼ gd ðql  qg Þ=r

ð32Þ

kl ¼ ll =lg kq ¼ ql =qg

where d denotes the initial bubble diameter. To assist with classifying simulations, other dimensionless parameters, such as dimensionless velocity U⁄, domain size X⁄, grid size Dx , time step Dt , Reynolds number (Re),viscosity l⁄, and density q⁄ are also defined as follows

U  ¼ u=U; X  ¼ x=d; Dx ¼ Dx=d; Dt  ¼ Dtðg=dÞ Re ¼ q1 ud=l1 ;



1=2

;



l ¼ l=l1 ; q ¼ q=q1

where U = (gd)1/2, d is the diameter of a circle bubble. The performance of handling interfacial problems involving large density ratio, such as, air/water interface, is an essential criterion to judge the robustness of a volume tracking method, because large gradients and material discontinuities exist near the interface which requires appropriate consideration [43]. In this

239

Z. Cao et al. / International Journal of Heat and Mass Transfer 111 (2017) 232–245

1

0.8

0.8

0.6

0.6

Y

Y

1

0.4

0.4

0.2

0.2

0

0

0.2

0.4

X

0.6

0.8

0

1

0

0.2

(a) t=0.5s

0.4

X

0.6

0.8

1

0.8

1

0.8

1

(b) t=1.0s

0.8

0.8

0.6

0.6

Y

1

Y

1

0.4

0.4

0.2

0.2

0

0

0.2

0.4

X

0.6

0.8

0

1

0

0.2

(c) t=2.0s

0.4

X

0.6

(d) t=4.0s

0.8

0.8

0.6

0.6

Y

1

Y

1

0.4

0.4

0.2

0.2

0

0

0.2

0.4

X

0.6

0.8

1

(e) t=6.0s

0

0

0.2

0.4

X

0.6

(f) t=8.0s

Fig. 10. Evolution of an initially circular fluid body in a single vortex flow simulated by the present VOSET method (grid size  1/100).

problem, a single gas bubble with the initial diameter of d (d = 0.01 m) is released from the position (2.5d, 2d) in a [0, 5d]  [0, 15d] liquid-filled domain. The free slip boundary conditions are applied on the surrounding walls. The density and viscosity ratios are all set as 1000, i.e. kq = kl = 1000. In this paper, two different cases are chosen for comparing the present unstructured VOSET method with the previous structured VOSET method [14].

These two cases are Eo = 10.0, M = 0.1 (case 1) and Eo = 100.0, M = 1000.0 (case 2), respectively. Two mesh sizes h = 1/20 and 1/40 corresponding to grids number 6.75  104 and 2.71  105 are adopted in the simulation, and the time step is calculated by Eq. (28). It should be pointed out that there are two large differences between the present unstructured VOSET method and the previous structured VOSET method. For the previous method, the

240

Z. Cao et al. / International Journal of Heat and Mass Transfer 111 (2017) 232–245

1

1

Initial geometry

Simulation result

0.9

0.8

0.7

0.7

Y

0.6

0.6

0.5

0.5

0.2

0.3

0.4

0.5

X

0.6

0.7

0.8

0.4

0.2

0.3

0.4

0.5

X

0.6

0.7

0.8

(b) VOF method (grid number=2.1×104)

(a) VOF method [6] 4

(grid number=0.35×10 ) 1

1

Simulation result

0.9

Initial geometry

0.7

0.7

Y

Y

0.8

Initial geometry

Simulation result

0.9

0.8

0.6

0.6

0.5

0.5

0.4

Initial geometry

Y

0.8

0.4

Simulation result

0.9

0.2

0.3

0.4

0.5

X

0.6

0.7

0.8

(c) LS method [31] (grid number=2.1×105)

0.4

0.2

0.3

0.4

0.5

X

0.6

0.7

0.8

(d) VOSET method (grid number=2.1×104)

Fig. 11. Interface locations after a time period in a time-reversed single vortex flow field.

volume fraction is evolved based on Eulerian framework and the surface tension force is calculated by the level set based CSF model. The present method adopts the Lagrangian-Eulerian remapping method to get the volume fraction and uses the level set based density-scared CSF model to calculate the surface tension force. In Fig. 14, the terminal bubble shapes and the flow fields computed by the present unstructured VOSET method are compared with the VOSET results on structured grids [14]. It can be found that both of the predicted simulation results agree well with each other, further verifying the accuracy of the present VOSET method. The rising velocity of the bubble is defined as the vertical component (vc) of the mean velocity:

vc ¼

R X R

av y dV X adV

ð33Þ

Fig. 15 shows the rising Reynolds number (Re = q1vcd/l1) of the single gas bubble with dimensionless time (t⁄ = t(g/d)1/2). The terminal velocities predicted by the current unstructured VOSET in case 1&2 are U⁄ = vc/(gd)1/2 = 0.4842 and 0.3680 respectively, whereas [14] reported terminal rising velocities of 0.1482 m/s and 0.1110 m/s,which are equivalent to U⁄ = 0.4734 and 0.3546. Hence, a small difference of 2.3% in case 1 and 3.8% in case 2 is

Fig. 12. Structural diagram of dam break problem.

Z. Cao et al. / International Journal of Heat and Mass Transfer 111 (2017) 232–245

Fig. 13. Comparison of the predicted liquid front location with the experimental data [40] and the numerical results [14,41].

observed between the present results and the structured VOSET results in [14]. To prove the mass conservative property of the proposed method, we investigate mass errors with time for the single gas bubble in these two cases (see Fig. 16), which shows that the maximum mass conservation error of the unstructured VOSET method is O(105). Since one of the main advantages for unstructured grids is that it can deal with complex domains, the following case is designed to

241

verify the present method’s adaptability to irregular regions. In this case, the dynamical behaviors of a single bubble through a converging channel are surveyed. The computational domain and initial configuration are shown in Fig. 17, a circular bubble of diameter d = 0.5 is initially placed in a converging channel of X = [0, 6d]  [0, 4d]. Both the bubble and its surrounding fluid are in the quiescent state at t = 0. No-slip boundary condition is imposed on all the walls. The physical parameters of the two fluids are listed in Table 1. The subscript 1 refers to the surrounding heavier fluid in X1 and the subscript 2 to the bubble in X2. In the case, the Eotvos number (Eo) equals to 90, and the Morton number (M) is 6.66  103. The time step size is determined by Eq. (28) with an upper bound of 4t⁄ = 4.43  104. To check whether our code solves the above problem correctly and also provide a way to make comparisons, the case is simulated by VOF in Fluent 6.3.26 with the same condition. Fig. 18 shows comparisons on the predicted bubble shapes at different times using the two different methods. The result shows that the shapes are quite similar with each other, and both bubbles tend to be a more non-convex shape and develop thin filaments at their tails after long time evolution. Moreover, there occurs minor difference in the rising velocities as the positions of two bubbles are approximately the same at different times. The temporal evolution of the Reynolds number and the center of mass of the bubble are shown in Fig. 19(a) and (b), respectively, it can be noted that how the Reynolds number converges as the mesh resolution increases from coarser grids (h = 1/20) to fine grids (h = 1/40). The maximum rise velocity (U⁄ = vc/(gd)1/2) has a magnitude of 0.4657 and occurs at dimensionless times 0.9695. When close to the converging wall,

(a) Case 1: Eo=10.0, M=0.1

(b) Case 2: Eo=100.0, M=1000.0 Fig. 14. Terminal bubble shapes and flow fields: (leftward panels) VOSET results on structured grids [14] (grid number = 3.00  104, dimensionless grid size h = 1/20); (rightward panels) present VOSET results on unstructured grids (grid number = 6.75  104, h = 1/20).

242

Z. Cao et al. / International Journal of Heat and Mass Transfer 111 (2017) 232–245

2d 5

2d

4

2d

3

Re

Ω1

4d

Γ

2 grid number: 6.75 104,h=1/20 1

3d

2d

5

grid number: 2.71 10 ,h=1/40

y x

0 0

2

4

6

8

10

12

14

t*

Ω2

d

d 6d

Fig. 17. Computational domain and initial configuration for the single bubble rising problem in a converging channel.

(a) Case 1: Eo=10.0, M=0.1

4. Conclusions 2.0

A coupled volume-of-fluid and level set (VOSET) method is proposed based on the unstructured triangular grids for interfacial flow simulations. The present unstructured VOSET method can not only guarantee mass conservation but also predict the surface tension with good accuracy, thus combining the advantages and overcoming the disadvantages of VOF and LS methods. To be specific, the main features of the VOSET method can be summarized as follows:

Re

1.5

1.0 grid number: 6.75 104,h=1/20

0.5

grid number: 2.71 105,h=1/40

0.0 0

2

4

6

8

10

12

14

t* (b) Case 2: Eo=100.0, M=1000.0

Er

Fig. 15. Curve of the rising Reynolds number (Re = q1ud/l1) of single gas bubble with dimensionless time (t⁄ = t (g/d)1/2).

1.0x10

-4

8.0x10

-5

6.0x10

-5

4.0x10

-5

2.0x10

-5

Case 1: Eo=10.0, M=0.1 Case 2: Eo=100.0, M=1000.0

0.0 0

2

4

6

8

10

12

14

t* Fig. 16. Variation curve of mass errors with time in the single gas bubble rising problem. Mass conservation error Er ¼ jMðtÞ  Mð0Þj=Mð0Þ. Here, R ! MðtÞ ¼ X ðaðx ; tÞÞdV.

the rising velocity of the bubble declines rapidly due to the effect of no-slip boundary condition.

(1) An iterative geometric operation is introduced based on unstructured triangular grids. Through the iterative geometric operation, we can, on one hand, reconstruct accurate and smooth interface, on the other hand, calculate accurate signed distance function, i.e. LS function. (2) The LS function is used to compute the accurate surface tension force by a density-scaled CSF model. Meanwhile, the discontinuous density and viscosity near the interfaces are also smoothed by the LS function. (3) The interface thickness plays a key role to calculate the accurate surface tension force and smooth the discontinuous physical parameters. The interface thickness should be set as different values at different locations because the unstructured triangular grids are usually non-uniform. We put forward a method to obtain the interface thickness on unstructured triangular grids. (4) The volume fraction advection is performed using a Lagrangian–Eulerian remapping algorithm, which can conserve the mass and overcome the mass-loss problem in the LS method. Finally, the present VOSET method is verified by four examples shown above. The results show that the present method can accurately simulate incompressible two-phase flows for unstructured triangular grids, proving its superiority. Thus it can be expected that the present VOSET method has promising potential to compute two-phase flow and heat transfer problems. In the future, we will extend the unstructured VOSET method to threedimensional simulations. The method for solving the distance function in three dimensions, the 3D analytic interface reconstruction and the 3D polyhedron clipping will present some challenges. Fortunately, the method computing the shortest distance from a point to a planar interface has been proposed by the present author Sun and his cooperators [21], the solution for the 3D analytic interface reconstruction on tetrahedral grids has been presented by Yang and James [30], and the Sutherland-Hodgeman algorithm

243

Z. Cao et al. / International Journal of Heat and Mass Transfer 111 (2017) 232–245 Table 1 Physical parameters defining the rising bubble test case.

q1

q2

l1

l2

g

r

q1/q2

l1/l2

1000

100

10

1

9.8

24.5

10

10

Time/methods

VOSET

VOF in Fluent

t*=0.89 (t=0.2)

t*=1.77 (t=0.4)

t*=2.66 (t=0.6)

t*=3.54 (t=0.8)

t*=4.43 (t=1)

Fig. 18. Bubble shapes for two dimensional rising bubble in a converging channel. Comparison between numerical predictions by VOSET and Fluent results (dimensionless grid size h = 1/20).

244

Z. Cao et al. / International Journal of Heat and Mass Transfer 111 (2017) 232–245

60 50

Re

40 30 20 4

grid number:1.58 10 ,h=1/20 4 grid number:6.21 10 ,h=1/40

10 0 0

1

2

3

4

t*

(a) 2.6 2.4 2.2

H*=yc/d

2.0 1.8 1.6 1.4 4

1.2

grid number:1.58 10 ,h=1/20 4 grid number:6.21 10 ,h=1/40

1.0 0.8

0

1

2

3

4

* t

(b) Fig. 19. The Reynolds number (Re = q1ud/l1) and center of mass of single gas bubble in a converging channel with dimensionless time (t⁄ = t (g/d)1/2). Here, R aydV . yc ¼ RX X

adV

adopted in this paper was also targeted to 3D polyhedron clipping [32]. The research regarding extensions to the case with heat transfer is also underway, and we will apply the present method to study nucleate pool boiling of FC-72 on micro-pin-finned surface under microgravity [44] in the future. Acknowledgements The study is supported by the National Natural Science Foundation of China (51325603, 51636006 and 51476054), the Program for New Century Excellent Talents in University (NCET-13-0792), the Beijing Talents Project and the Foundation of Shaanxi Key Laboratory of Energy Chemical Process Intensification (SXECPI201501). References [1] C.W. Hirt, D.B. Nichols, Volume-of-fluid (VOF) method for the dynamics of free boundaries, J. Comput. Phys. 39 (1981) 201–205. [2] J.E. Pilliod, E.G. Pucket, Second-order accurate volume-of-fluid algorithm for tracking material interfaces, J. Comput. Phys. 199 (2004) 465–502.

[3] J. López, J. Hernández, P. Gómez, et al., A volume of fluid method based on multidimensional advection and spline interface reconstruction, J. Comput. Phys. 195 (2004) 718–742. [4] C.B. Ivey, P. Moin, Conservative volume of fluid advection method on unstructured grids in three dimensions, center for turbulence research, Ann. Res. Briefs (2012) 179–192. [5] J. Dukowicz, J.R. Baumgardner, Incremental remapping as a transport/ advection algorithm, J. Comput. Phys. 160 (2000) 318–335. [6] K. Shahbazi, M. Paraschivoiu, J. Mostaghimi, Second order accurate volume tracking based on remapping for triangular meshes, J. Comput. Phys. 188 (2003) 100–122. [7] S. Osher, J.A. Sethian, Fronts propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulations, J. Comput. Phys. 79 (1988) 12–49. [8] E. Olsson, G. Kreiss, S. Zahedi, A conservative level set method for two phase flow, J. Comput. Phys. 210 (2005) 225–246. [9] C.H. Yu, Z.T. Ye, T.W.H. Sheu, et al., An improved interface preserving level set method for simulating three dimensional rising bubble, Int. J. Heat Mass Transfer 103 (2016) 753–772. [10] A.M. Moghadam, M. Shafieefar, R. Panahi, Development of a high-order level set method: compact conservative level set (CCLS), Comput. Fluids 129 (2016) 79–90. [11] T.Y. Hou, J.S. Lowengrub, M.J. Shelley, Boundary integral methods for multicomponent fluids and multiphase materials, J. Comput. Phys. 169 (2001) 302– 362. [12] S.O. Unverdi, G. Tryggvason, A front-tracking method for viscous, incompressible multi-fluid flows, J. Comput. Phys. 100 (1992) 25–37. [13] G. Tryggvasson, B. Bunner, A. Esmaeeli, A front tracking method for the computations of multiphase flow, J. Comput. Phys. 169 (2001) 708–759. [14] D.L. Sun, W.Q. Tao, A coupled volume-of-fluid and level set (VOSET) method for computing incompressible two-phase flows, Int. J. Heat Mass Transfer 53 (2010) 645–655. [15] D. Gerlach, G. Tomar, G. Biswas, et al., Comparison of volume-of-fluid methods for surface tension-dominant two-phase flows, Int. J. Heat Mass Transfer 49 (2006) 740–754. [16] S. Osher, R.P. Fedkiw, Level set methods: an overview and some recent results, J. Comput. Phys. 169 (2001) 463–502. [17] T.W.H. Sheu, C.H. Yu, P.H. Chiu, Development of level set method with good area preservation to predict interface in two-phase flows, Int. J. Numer. Meth. Fluids 67 (2011) 109–134. [18] M. Sussman, E.G. Puckett, A coupled level set and volume-of-fluid method for computing 3D and axisymmetric incompressible two-phase flows, J. Comput. Phys. 162 (2000) 301–337. [19] X.F. Yang, A.J. James, An adaptive coupled level-set/volume-of-fluid interface capturing method for unstructured triangular grids, J. Comput. Phys. 217 (2006) 364–394. [20] I. Chakraborty, M.R. Rubio, A. Sevilla, et al., Numerical simulation of axisymmetric drop formation using a coupled level set and volume of fluid method, Int. J. Multiphase Flow 84 (2016) 54–65. [21] K. Ling, Z.H. Li, D.L. Sun, et al., A three-dimensional volume of fluid & level set (VOSET) method for incompressible two-phase flow, Comput. Fluids 118 (2015) 293–304. [22] T. Wang, H.X. Li, Y.C. Feng, et al., A coupled volume-of-fluid and level set (VOSET) method on dynamically adaptive quadtree grids, Int. J. Heat Mass Transfer 67 (2013) 70–73. [23] T. Wang, H.X. Li, Y.F. Zhang, et al., Numerical simulation of bubble dynamics in a uniform electric field by the adaptive 3D-VOSET method, Numer. Heat Transfer, Part A 67 (2015) 1352–1369. [24] M.R. Ansari, R. Azadi, E. Salimi, Capturing of interface topological changes in two- phase gas–liquid flows using a coupled volume-of-fluid and level-set method (VOSET), Comput. Fluids 125 (2016) 82–100. [25] D.Z. Guo, D.L. Sun, Z.Y. Li, et al., Phase change heat transfer simulation for boiling bubbles arising from a vapor film by the VOSET method, Numer. Heat Transfer, Part A 59 (2011) 857–881. [26] K. Ling, Z.Y. Li, W.Q. Tao, A direct numerical simulation for nucleate boiling by the VOSET method, Numer. Heat Transfer, Part A 65 (2014) 949–971. [27] K. Ling, G. Son, D.L. Sun, et al., Three dimensional numerical simulation on bubble growth and merger in micro-channel boiling flow, Int. J. Therm. Sci. 98 (2015) 135–147. [28] D.X. Shi, Q.C. Bi, R.Q. Zhou, Numerical simulation of a falling ferrofluid droplet in a uniform magnetic field by the VOSET method, Numer. Heat Transfer, Part A 66 (2014) 144–164. [29] M. Huang, L.L. Wu, B. Chen, A piecewise linear interface-capturing volume-offluid method based on unstructured grids, Numer. Heat Transfer, Part B 61 (2012) 412–437. [30] X.F. Yang, A.J. James, Analytic relations for reconstructing piecewise linear interfaces in triangular and tetrahedral grids, J. Comput. Phys. 214 (2006) 41– 54. [31] N. Balcázar, L. Jofre, O. Lehmkuhl, et al., A finite-volume/level-set method for simulating two-phase flows on unstructured grids, Int. J. Multiphase Flow 64 (2014) 55–72. [32] I.E. Sutherland, G.W. Hodgeman, Reentrant polygon clipping, Commun. ACM 17 (1974) 32–42. [33] K. Yokoi, A practical numerical framework for free surface flows based on CLSVOF method, multi-moment methods and density-scaled CSF model:

Z. Cao et al. / International Journal of Heat and Mass Transfer 111 (2017) 232–245

[34] [35]

[36] [37] [38] [39]

numerical simulations of droplet splashing, J. Comput. Phys. 232 (2013) 252– 271. J.U. Brackbill, D.B. Kothe, C. Zemach, A continuum method for modeling surface tension, J. Comput. Phys. 100 (1992) 335–354. M.M. Francois, S.J. Cummins, E.D. Dendy, et al., A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework, J. Comput. Phys. 213 (2006) 141–173. G.J. Li, M. Huang, B. Chen, VOF method on unstructured grid using PISO algorithm, J. Eng. Thermophys. 34 (2013) 476–479. S.T. Zalesak, Fully multidimensional flux-corrected transport algorithms for fluids, J. Comput. Phys. 31 (1979) 335–362. J.B. Bell, P. Colella, H.M. Glaz, A second order projection method of the incompressible Navier-Stokes equations, J. Comput. Phys. 85 (1989) 257–283. B.M. Ningegowda, B. Premachandran, A coupled level set and volume of fluid method with multi-directional advection algorithms for two-phase flows with and without phase change, Int. J. Heat Mass Transfer 79 (2014) 532–550.

245

[40] J.C. Martin, W.J. Moyce, An experimental study of the collapse of liquid columns on a rigid horizontal plane, Philos. Trans. Roy. Soc. Lond, Ser. A, Math. Phys. Sci. 244 (1952) 312–324. [41] P.H. Chiu, Y.T. Lin, A conservative phase field method for solving incompressible two-phase flows, J. Comput. Phys. 230 (2011) 185–204. [42] J.R. Grace, Shapes and velocities of bubbles rising in infinite liquids, Trans. Inst. Chem. Eng. 51 (1973) 116–120. [43] J.M. Cubos-Ramírez, J. Ramírez, M. Salinas-Vázquez, et al., Efficient two-phase mass-conserving level set method for simulation of incompressible turbulent free surface flows with large density ratio, Comput. Fluids 136 (2016) 212– 227. [44] Y.F. Xue, J.F. Zhao, J.J. Wei, et al., Experimental study of nucleate pool boiling of FC-72 on micro-pin-finned surface under microgravity, Int. J. Heat Mass Transfer 63 (2013) 425–433.