On: 25 April 2014, At: 23:59 Publisher: Taylor & Francis

Informa Ltd Registered in England and Wales Registered Number: 1072954 Registered office: Mortimer House, 37-41 Mortimer Street, London W1T 3JH, UK

**Numerical Heat Transfer, Part B:**

**Fundamentals: An International Journal**

**of Computation and Methodology**

Publication details, including instructions for authors and subscription information:

http://www.tandfonline.com/loi/unhb20

**A VOF-Based Conservative Interpolation**

**Scheme for Interface Tracking (CISIT) of**

**Two-Fluid Flows**

Yeng-Yung Tsui a & Shi-Wen Lin a a

Department of Mechanical Engineering , National Chiao Tung University , Hsinchu , Taiwan , Republic of China

Published online: 15 Apr 2013.

**To cite this article: Yeng-Yung Tsui & Shi-Wen Lin (2013) A VOF-Based Conservative Interpolation**
Scheme for Interface Tracking (CISIT) of Two-Fluid Flows, Numerical Heat Transfer, Part B:
Fundamentals: An International Journal of Computation and Methodology, 63:4, 263-283, DOI:

10.1080/10407790.2013.756251

**To link to this article: **http://dx.doi.org/10.1080/10407790.2013.756251

PLEASE SCROLL DOWN FOR ARTICLE

Taylor & Francis makes every effort to ensure the accuracy of all the information (the “Content”) contained in the publications on our platform. However, Taylor & Francis, our agents, and our licensors make no representations or warranties whatsoever as to the accuracy, completeness, or suitability for any purpose of the Content. Any opinions and views expressed in this publication are the opinions and views of the authors, and are not the views of or endorsed by Taylor & Francis. The accuracy of the Content should not be relied upon and should be independently verified with primary sources of information. Taylor and Francis shall not be liable for any losses, actions, claims, proceedings, demands, costs, expenses, damages, and other liabilities whatsoever or howsoever caused arising directly or indirectly in connection with, in relation to or arising out of the use of the Content.

This article may be used for research, teaching, and private study purposes. Any substantial or systematic reproduction, redistribution, reselling, loan, sub-licensing, systematic supply, or distribution in any form to anyone is expressly forbidden. Terms & Conditions of access and use can be found at http://www.tandfonline.com/page/terms-and-conditions

### A VOF-BASED CONSERVATIVE INTERPOLATION

### SCHEME FOR INTERFACE TRACKING (CISIT) OF

### TWO-FLUID FLOWS

### Yeng-Yung Tsui and Shi-Wen Lin

Department of Mechanical Engineering, National Chiao Tung University, Hsinchu, Taiwan, Republic of China

A simple interface-reconstruction scheme (CISIT), based on the volume-of-ﬂuid (VOF) formulation and applicable to unstructured grids with arbitrary topology, is developed to track the interface in free-surface ﬂows. The interface is represented by the contour surface of VOF value f¼ 0.5, which is constructed simply by interpolation. The advancing of the interface is made by imposing mass conservation in the cells containing the interface. How-ever, problems arise when the interface moves across the cell face and into another cell. The cell originally occupied by the interface may be overﬁlled (f > 1) or underﬁlled (f < 1) as the interface advances (moving in the positive direction of the interface) to leave the cell. In contrast, the cell may be overdepleted (f < 0) or underdepleted (f > 0) as the interface retreats (moving in the negative direction of the interface) to move out of the cell. These situations are remedied via adjusting the transported volume of ﬂuid across the cell face by following the conservation law so that the VOF becomes 1 or 0 in these cells, depending on advancing or retreating of the interface. The resulting VOF distribution is uniform, either in 1 or 0, in the region outside the interface cells, and the interface occupies only one cell in its width. To smooth the ﬂuid properties across the interface region in the velocity calculation, an average smoothing technique is adopted. Application to a number of test cases, including model cases with known velocities and a number of real ﬂows, reveal that this method is accurate and robust.

INTRODUCTION

Multiﬂuid or multiphase ﬂow features an interface separating ﬂuids with differ-ent ﬂuid properties. Accurate prediction of this interfacial ﬂow presdiffer-ents a challenging task in numerical simulation because there exists a jump in ﬂuid density and viscosity across the interface. This jump phenomenon can be identiﬁed as a Heaviside (step) function in theory. Different from the high-speed compressible ﬂow, in which the ﬂow velocity, pressure, and temperature present a discontinuity across the shock

Received 3 October 2012; accepted 5 November 2012.

This work was supported by the National Science Council of the Rebublic of China under Contract Number 98-2221-E-009-130. The authors are indebted to Dr. J.-T. Yeh of Material and Chemical Research Laboratories, Industrial Technology Research Institute, for his encouraging and helpful discus-sions throughout this study.

Address correspondence to Yeng-Yung Tsui, Department of Mechanical Engineering, National Chiao Tung University, 1001 Ta-Hsueh Road, Hsinchu 300, Taiwan, Republic of China. E-mail: yytsui@mail.nctu.edu.tw

*Copyright # Taylor & Francis Group, LLC*
ISSN: 1040-7790 print=1521-0626 online
DOI: 10.1080/10407790.2013.756251

263

wave, these variables are continuous in interfacial ﬂows. Direct employment of shock-capturing schemes may not be appropriate for these ﬂows.

The geometric shape of the interface continues to change during ﬂuid ﬂow. The large surface distortion and surface breaking and merging rule out the use of Lagragian-grid methods in which grids are attached to the interface. An alternative is to adopt Eulerian grids (ﬁxed grids) with a method to track the movement of the interface. Tryggvason and co-workers [1, 2] developed a front-tracking method to fulﬁll this job. A ﬁxed grid is assigned to the velocity ﬁeld and a grid of lower dimen-sion to the ﬂuid front, or the interface. The interface grid is composed of triangular elements for a 3-D surface or line segments for a 2-D curve. The ﬂuid front moves in a Lagrangian manner, with the propagating velocity obtained by interpolation from the velocities on the Eulerian mesh. Since the interface grid may be greatly distorted in the front movement, re-meshing is required from time to time. An obvious short-coming is that the masses of the ﬂuids are not conserved during front propagation. In general, the coding of the method is rather complicated due to the interaction between the interface grid and the Eulerian grid, along with the requirement of re-meshing for the former.

Another method to track the interface movement is based on the use of a level function deﬁned as a signed distance away from the interface. The interface is repre-sented by the set of zero-level value points. Therefore, it was called the level set method by Osher and Sethian [3]. Different from the above front-tracking method, the advancing of the interface is made by solving the advection equation for the level set function U on an Eulerian grid as

### NOMENCLATURE

C volumetric ﬂuid ﬂux Cn Courant number Eo Eotvos number f volume fraction F ﬂux Fr Froude number gi gravitational acceleration Mo Morton number P pressure r gradient ratio Re Reynolds number S surface area ~ S S surface vector ~ T

T; Ti surface tension force

~ V

V ; Vj ﬂow velocity vector

w weighting factor
c(r) ﬂux limiter function
~_{d}_{d} _{distance vector}
Dt time-step size
Dv cell volume
m ﬂuid viscosity
q ﬂuid density

r surface tension coefﬁcient

sij viscous stress

/ Cartesian velocity components U level set function Subscripts

D downstream cell f cell face ﬁ ith cell face

g gas

l liquid

Ni ith neighboring cell P primary cell U upstream cell UU far upstream cell v vertex of the primary cell vi ith vertex

Superscripts

c convection d diffusion n new time step o old time step

w wetted area of the face

qU

qt þ ~VV rU ¼ 0 ð1Þ where ~VV is the velocity of the ﬂow. Hence, it is much easier to implement than the former one. A drawback is that it is prone to numerical errors in solving the hyper-bolic equation when the interface experiences severe stretching or tearing. Similar to the explicit tracking method described above, the major weakness of the method lies in the lack of mass conservation, especially when the interface is highly distorted. One of the sources that cause mass nonconservation is that the level set function will no longer be a distance function after some period of time [4]. Therefore, a redistan-cing or reinitialization becomes necessary to maintain it as a distance function after each time step. It has been shown that after the reinitialization process is done, the error of mass loss or gain is largely reduced [5, 6]. Another tactic to reduce the mass conservation problem is to reduce spatial discretization errors using high-order schemes together with adaptive grid technology [7].

Another class of methods using Eulerian grids is based on the volume-of-ﬂuid (VOF) function, which is deﬁned as the fraction of volume in the mesh cell occupied by one of the two ﬂuids. The VOF function f is either unity or zero in the cells con-taining a single ﬂuid and less than 1 and greater than 0 in the cells where the interface is located. Its evolution follows the same form given in Eq. (1).

qf

qtþ ~VV rf ¼ 0 ð2Þ Introducing the continuity equation, it can be cast into a conservation form.

qf

qtþ r ð~VV fÞ ¼ 0 ð3Þ Although eqs. (1) and (2) are of the same form, it is more difﬁcult to solve the equation for the VOF function than that for the level set function due to the former being a step function while the latter is continuous and smooth. A variety of schemes have been developed in the past, which, in general, can be classiﬁed into two cate-gories. One is to solve the VOF in a way similar to the shock capturing for super-sonic ﬂows. It is essential for this kind of method to limit numerical diffusion and dispersion which occur in the region near the jump. It is known that the ﬁrst-order upwind scheme is most stable but causes large diffusion, while the ﬁrst-order down-wind scheme has negative artiﬁcial viscosity. Although the downdown-wind scheme is unstable, the antidiffusive character tends to compress the jump region to make it sharp. Hirt and Nichols [8] were the ﬁrst to combine the upwind and downwind schemes in such a way as to minimize diffusion in interfacial ﬂow calculations. Rudman [9] adopted a method based on the ﬂux-corrected transport algorithm [10]. It consists of two stages. In the ﬁrst stage, the VOF value is calculated using the upwind scheme. This is followed by a correction stage in which the downwind scheme is used to reduce the diffusion generated in the ﬁrst stage. High-order schemes have also been introduced, including the HRIC of Muzaferija and Peric [11], CICSAM of Ubbink and Issa [12], STACS of Darwish and Moukalled [13],

and the latest FBICS of Tsui and Wu [14]. These schemes are based on the NVD formulation and satisfy the convection boundedness criteria [15]. They blend a high-resolution (bounded high-order) scheme with a compressive (bounded-downwind) scheme. The blending strategy is to switch from the high-resolution scheme to the compressive scheme in a continuous manner, depending on the angle between the interface and the grid lines. Another merit of these schemes is that they are developed with the aim of application to unstructured grids. The width of the resulting interface can be restricted to a region of several mesh sizes.

In another category of VOF method, the interface is approximately recon-structed using VOF information in local cells. The ﬂuid ﬂuxes crossing the faces of these cells are then determined with the help of the interface geometry in the cells. A simple way is to reconstruct the interface such that it is aligned with the grid line. In the SLIC algorithm of Noh and Woodward [16], this approximate interface is assumed to be vertical or horizontal, depending on the sweep direction of the direction-splitting solution method employed in this algorithm. A higher-order rep-resentation, introduced by Young as PLIC [17], is to approximate it by an oblique line segment. Comparing with the SLIC, the PLIC is more accurate, but more difﬁ-cult to implement. To construct the surface using the available volume fraction and its gradient (surface normal vector) in local cells, the use of iterative procedures is most common [18, 19]. It is further complicated by having a large number of possible interface conﬁgurations that need to be considered for determination of ﬂuid ﬂuxes through the cell faces. Different from the PLIC, in which the interface is ﬁtted inside the cell, the interface represented by a line segment is drawn across the face boundary in the FLAIR algorithm [20]. The linear equation of the approximate surface is determined by the volume fractions of the two cells on either side of the face. The ﬂuid ﬂux through this face during a time step can then easily be obtained by inte-gration. This algorithm was successfully applied to 2-D and axisymmetric ﬂow prob-lems [21, 22]. The major limitation of this technique is that it cannot be extended to be applicable to 3-D problems.

The above volume-tracking methods have become widely used, partly due to their capability to handle complex free-surface ﬂow problems in a relatively easy way. The accuracy of the interface-reconstruction VOF methods is highly dependent on the calculation of the interface normal and curvature from volume fractions, which is difﬁcult to compute accurately because of the discontinuity property of this function. In contrast, the level set function is continuous and smooth. The deriva-tives of this smooth function can be easily calculated with high-order accuracy. How-ever, as noted, the ﬂuid mass is not conserved. To take advantage of the merits of the VOF and level set methods, it is natural to combine the two to form a coupled method [23–27]. In these kinds of methods, the interface normal and curvature are evaluated from the level set function and the evolution of the surface is made using the VOF to enforce mass conservation.

The PLIC algorithm is one of the most popular methods used in many applica-tions. However, the rebuilding of the interface from the VOF is not straightforward even for rectangular grids. In addition, as required for the ﬂuxing across cell faces, it is necessary to consider 16 interface conﬁgurations for 2-D ﬂows and 64 conﬁgura-tions for 3-D ﬂows, though these numbers of conﬁguraconﬁgura-tions can be reduced after reg-ulating the cell geometry [28]. The extension of the method for unstructured meshes

with irregular cell geometry is not easy, partly due to the increased difﬁculty in reconstructing the interface and also due to the increased complexity in estimating the ﬂow ﬂuxes. To circumvent the latter problem, a Lagrangian-Eulerian advection algorithm was introduced by Mosso et al. [29]. In this method, the Eulerian grid cells containing the interface are moved in the Lagrangian sense. This is followed by reconstructing the interface in these cells because they are deformed during advec-tion. The last step is to remap the volume fraction from the Lagrangian cells to the original Eulerian cells. This method has been applied to triangular grids by Shahbazi et al. [30] and Yang et al. [26]. Instead of the mesh cell, it is the ﬂuid portion in the cell being advected in the study by Ashgriz et al. [31].

In view of the above survey, methods that are more accurate or more robust become more complicated. In this study, a simple method, within the framework of the VOF, is developed. It is robust because it is applicable to unstructured grids with arbitrary geometry and equally applicable to 3-D problems without causing any further complication. The interface is based on the contour surface of volume fraction 0.5, which can be easily constructed by interpolation. The advection step is fulﬁlled in a prediction-and-correction way. In the following, details of this solution method are described ﬁrst. This is followed by a brief description of the method for the velocity ﬁeld. Then, ﬁdelity of the algorithm and veriﬁcation of the accuracy are presented. Finally, concluding remarks are given.

SOLUTION METHOD FOR VOF

The present solution method is inspired by the work of Yeh [32, 33], in which the interface is represented by the contour line f¼ 0.5. In Yeh’s study, the VOF value is originally stored at the centroid of each mesh cell. To calculate the ﬂuid ﬂux, the VOF is sought on grid nodes where velocities are located. The VOF on the grid node is ﬁrst obtained using interpolation from the surrounding cells which share this node as a common vertex. The nodal value is then reset to 1 if the volume fraction is great-er than 0.5 and to 0 if it is less than 0.5. The advection equation is solved using the ﬁnite-element method. An obvious shortcoming of this method is that the VOF is not conserved. A conservative algorithm is proposed in the following.

Interface Reconstruction

In the solution, the contour surface, or contour line in 2-D problems, of f¼ 0.5 is used as an indicator of the interface. It can easily be reconstructed by an interp-olation practice. The original values of f are stored at the centroid of each grid cell. A linear weighting is employed to ﬁnd fvat each cell vertex,

fv¼ P jwjfj P jwj ð4Þ

where the subscript j denotes the cells surrounding the grid node and wjis the

weight-ing factor correspondweight-ing to the cell j. The calculation of the weightweight-ing factor can be

based on either the inverse of the volumes of the surrounding cells or the inverse of the distances between this node and those cell centroids.

After the nodal VOFs are obtained for the two end vertices of each edge of a cell, it is checked whether one of them is less than and the other greater than the value 0.5. If this is true, the interface intersects with this edge. The intersecting point, at which the VOF is 0.5, can then be determined by linear interpolation from the two nodes. After examining all the edges, the interface in the cell can be determined. This procedure proceeds in each cell of the mesh. In the end, a continuous piecewise-linear interface is constructed. This process is equally applicable to 2-D and 3-D cells.

Time Advancing

For a cell in a Eulerian mesh, the VOF Eq. (3) in the conservation form can be discretized using the ﬁnite-volume method as

DV Dt ðf n P f o PÞ þ X j Ffj ¼ 0 ð5Þ

Here the superscripts n and o denote the new and old time steps, respectively, the subscript P designates the considered cell, the subscript f indicates a face of the cell, Ff is the ﬂux through this face, and the subscript j denotes the face number. The

summation is over all the surrounding faces of the cell. The scheme being explicit or implicit depends on whether the ﬂuxes are evaluated using new time-level or old time-level values. In the following, the explicit scheme is employed. It is noted that the advection Eq. (2) rewritten in the following form was used for discretization in some studies.

qf

qtþ r ð~VV fÞ ¼ ðr ~VVÞf ð6Þ The reason for using this nonconservation expression is that the velocities on the cell faces do not satisfy the continuity constraint in these methods. The term on the right-hand side simply eliminates this error. For the method employed in this study, mass conservation is strictly obeyed in each cell.

In interface-capturing methods, the VOF ﬂux is approximated as

Ff ¼ Cfff ¼ ð~VVf ~SSfÞff ð7Þ

where ~VVf and ~SSf are the velocity and surface vectors at the face, respectively, Cfis the

volumetric ﬂuid ﬂux, and ffis the ﬂuxing value of the VOF. By these algorithms, the

issue of controlling numerical diffusion and dispersion becomes most important, which depends on the approximate estimation of the ﬂuxing value ffusing the VOFs

of neighboring cells [11–14].

In the present interface-reconstruction method, the face ﬂux is the volumetric ﬂux through the portion of the face wetted by the ﬂuid:

Ff ¼ ~VVf ~SSfw ð8Þ

where ~SSw

f represents the wetted area on the face. As an illustration, a 2-D cell is

shown in Figure 1. The intersecting of the interface with the edges of the cell leads to three wetted areas. These wetted areas function as inlets or outlets to the cell for the ﬂuid. The conservation law applied to this control volume results in the following equation: DV Dt ðf n P fPoÞ þ X j ~ V Vfj ~SSwfj ¼ 0 ð9Þ

It can be seen that there is no need to estimate the face value ff, as required in Eq. (7),

and the calculation is undertaken only in the cells where the interface is located. Thus, numerical diffusion does not prevail in this method.

Cell Adjustment

It is ideal for a simulation based on VOF methods that the interface region covers a cell width and the VOF is either zero or one in the rest of the regions. The solution of the above time-marching step does not meet this ideal situation. In general, four kinds of problems may be encountered which occur when the inter-face leaves a cell and moves into another cell.

Figure 1. A control volume partially ﬁlled with ﬂuid (color ﬁgure available online).

1. Overfilling (fp> 1). After solving Eq. (9), it is possible for the value of fPn

to be greater than 1. This situation can be identiﬁed from Figure 2a, in which the inter-face is advancing to the right in a 1-D ﬂow. In the beginning, the VOF value in the cell increases with the advancement of the interface, because ﬂuid ﬂows into the cell from the left face and no ﬂuid leaves from the right face. By the time the interface crosses the right face, the VOF becomes greater than 1 because the wetted area on the right face remains zero. To force the interface to move across the cell face, the excessive ﬂuidDf(¼ fP 1)must be reallocated to the downstream cell N1. For a

multi-dimensional ﬂow, as shown in Figure 2b, the excess must be portioned and assigned to the two downstream cells N1and N2. The weighting factor wNifor the portioning

is based on the ratio of the efﬂux through the corresponding face to the total efﬂux:

wN1¼ Cf 1 Cf 1þ Cf 2 wN2¼ Cf 2 Cf 1þ Cf 2 ð10Þ

where Cfi ð¼~VVfi ~SSfiÞ is the volumetric ﬂuid ﬂux through the face fi. The formulation

for the weighting factor can be written in a general form, which can be applied to all neighboring cells: wNi¼ maxðCfi;0Þ P jmaxðCfj;0Þ ð11Þ

Here the sum is over all the faces. It is obvious that the weighting factor is zero for upstream cells. The VOFs for the downstream cells are corrected by

fNi¼ fNiþ wNiðfP 1Þ

DvP

DvNi

ð12Þ Figure 2. Advancing of interface through the face of a control volume (color ﬁgure available online).

After the neighboring cells are adjusted, the VOF for this overﬁlling cell is reset to 1. 2. Overdepleting (fp< 0). When the interface retreats from a cell, as seen in

Figure 3a in a 1-D ﬂow, the VOF becomes less than zero at the time of crossing the right face. The overdepleted volume of ﬂuid must be retrieved from the downstream cell N1. Similar to the overﬁlling situation, the overdepleted volume needs to be portioned and reallocated for multidimensional problems such as the one shown in Figure 3b. The weighting factor is given by Eq. (11). The VOFs for the neigh-boring cells are adjusted in the following way:

fNi¼ fNiþ wNifP

DvP

DvNi

ð13Þ

Afterwards, zero value is assigned to fP.

3. Underfilling (fp< 1). A shear ﬂow has large velocity gradients in the

trans-verse direction. In such ﬂows, different from overﬁlling, the value of VOF may become less than 1 when the interface advances to leave the cell (Figure 2b). Since the interface does not lie in this cell any more, the VOF value remains less than 1 in later time. Hence, whenever it is found that fP<1 while the interface moves out

of this cell, i.e., fvi>0.5 (fvirepresents the VOFs of all the cell vertices), ﬂuid must

be retrieved from the downstream cells to ﬁll this cell so that fP is equal to 1. The

weighting factor and the adjustment of the VOF are the same as those shown in Eqs. (11) and (12).

4. Underdepleting (fp> 0). Another situation encountered in high-shear

ﬂows is that the VOF may remain greater than zero when the interface retreats to leave the cell (Figure 3b). Similar to the underﬁlling cell, the volume of ﬂuid remains unchanged and the ﬂuid in the cell becomes stagnant. Therefore, whenever it is found that fP>0 and the interface is not situated in this cell any more (fvi<0.5 for all cell

vertices), ﬂuid in this cell must be allocated to the downstream cells so that Figure 3. Retreating of interface through the face of a control volume (color ﬁgure available online).

fP becomes zero. The weighting factor and the adjustment are given by Eqs. (11)

and (13).

The solution procedure described above is summarized as follows. 1. The interface is reconstructed ﬁrst using the VOF from the last time step. 2. New values of VOF are sought via conserving ﬂuid mass in interface cells. 3. The VOFs are adjusted to ensure that no cells are overﬁlled, underﬁlled,

overdepleted, or underdepleted. It is noted that after the correction process for underﬁlling and underdepleting cells, it is possible for some of the neighboring cells to be overﬁlled or overdepleted. More similar correction steps can then be undertaken. Usually, two sweeps of such adjustment are enough to wipe out this problem.

This completes the calculation of the VOF in one time step. It can be seen that the conservation law is obeyed in each of the predictor and corrector steps. The resulting interface region is sharp and lies in one cell only.

SOLUTION METHOD FOR VELOCITY FIELD

In each time step of calculation, the VOF is obtained ﬁrst, followed by solving for the velocity ﬁeld. The ﬂow in both ﬂuids is assumed to be laminar and incom-pressible. The conservation of mass and momentum for the two ﬂuids can be written in one set of equations:

qVj qxj ¼ 0 ð14Þ qqVi qt þ q qxj ðqVjViÞ ¼ qP qxi þqsij qxj þ qgiþ Ti ð15Þ

where Vj is the velocity, q the density, P the pressure, sij the viscous stress, gi the

gravitational acceleration, and Ti the force due to surface tension. Surface tension

is a surface force at the free interface in two-phase ﬂows. With the CSF (continuum surface force) formulation [34], it can be modeled as a body force:

~ T T¼ rr rf rf j j rf ð16Þ

where r is the surface tension coefﬁcient.

Smoothing of Fluid Properties

In the one-equation model for two-ﬂuid ﬂows, the ﬂuid properties, such as density and viscosity, are expressed as functions of the local volume fraction:

q¼ f q1þ ð1 f Þq2 ð17aÞ

m¼ f m1þ ð1 f Þm2 ð17bÞ

where the subscripts 1 and 2 denote the two ﬂuids. When the difference between the ﬂuid properties, especially the density, of the two ﬂuids is large, the sharp gradient across the interface region may cause numerical instability. To alleviate this problem, the gradient of the properties needs to be smeared artiﬁcially. Various smoothing techniques can be found in [5, 35, 36]. In this study, an averaging smoother is used. First, the property values at all cell vertices are calculated using the interpolation method given by Eq. (4). A new value at the cell centroid is then obtained by aver-aging over these vertices. In this manner, the interface region is expanded to the two adjacent cells on either side of the interface cell. This smoothing process can be repeated to enhance the smearing effects. In our calculations, no more than two such smoothing steps are taken.

Numerical Method

Discretization of the equations is performed using the ﬁnite-volume method. The convective momentum ﬂux is expressed as

Fc¼ ðqfVV~f ~SSfÞ/f ð18Þ

where / designates the Cartesian velocity components. The face ﬂux value /f is

approximated using neighboring nodal values,
/_{f} ¼ /_{U}þcðrÞ

2 ð/D /UÞ ð19Þ where, as seen in Figure 4, the subscript U denotes the cell node upstream of the considered face and the subscript D is the one downstream. The symbol c represents the ﬂux limiter, which depends on the gradient ratio r deﬁned by

r¼/U /UU /D /U

ð20Þ

Figure 4. Illustration for determination of the far upstream node UU.

Here /UUis the value at a node far upstream (see Figure 4), estimated by

/UU ¼ /D 2 r/ð ÞU~ddUD ð21Þ

where (r/)Udenotes the gradient at the upstream node U and ~ddUD is the distance

vector directing from U to D. A number of TVD and NVD schemes are available by assigning different expressions to the function c(r) [37, 38]. In the following cal-culations, the Van Leer scheme is adopted, which is given by

cðrÞ ¼rþ rj j

rþ 1 ð22Þ

The diffusive momentum ﬂux, suitable for unstructured grids, is approximated by

Fd ¼ mS
2
f
~_{d}_{d}
PN ~SSf
ð/N /PÞ þ m r/ð Þf ~SSf
S2
f
~_{d}_{d}
PN ~SSf
~_{d}_{d}
PN
!
ð23Þ

where the subscripts P and N denote the considered cell and neighboring cell, respectively, ~ddPN is the vector connecting nodes P and N, and (r/)fis the gradient

at the face.

All variables, including all velocity components and pressure, are stored at the centroid of each grid cell. The coupling between momentum equations and the con-tinuity equation is tackled using the noniterative technique of the PISO algorithm [39]. This is a kind of predictor-corrector method. First, the momentum equations are solved to obtain velocities at cell centroids in the predictor step using the prevail-ing velocity and pressure ﬁelds. A momentum interpolation method is then employed to ﬁnd the velocities on the cell face and the mass ﬂux through the face. A pressure-correction equation can be obtained by forcing the corrected mass ﬂux to obey the conservation law in each cell. After the pressure correction is solved for, the pressure and velocity are upgraded accordingly. This adjustment procedure takes place in two corrector steps, after which both the momentum and continuity equa-tions are better satisﬁed. The way to derive the face mass ﬂux and the pressure-correction equation can be found in [40].

RESULTS AND DISCUSSION

The solution method for the VOF is ﬁrst tested in two model problems with given velocities. One is the advection of a square or a circular cylinder in a uniform velocity ﬁeld and the other is the advection of a circular cylinder in a shear ﬂow. Real ﬂows are then considered. The realistic ﬂows selected for testing include collapse of a water column, rising of a gas bubble, and the Rayleigh-Taylor instability problem.

Advection of Cylinders in a Uniform Velocity Field

The computational domain is a unit square with a constant velocity ~

V

V ¼ u~iiþ v~jj being imposed. Either a square cylinder or a circular cylinder is initially

placed at a position x¼ 0.25 and y ¼ 0.25. The cylinder has a side of length (or
diam-eter for the circle) 0.3. It is moved to the position corresponding to x¼ 0.75. The
domain is partitioned into various rectangular and triangular grids. Theoretically,
the geometry of the cylinders remains unchanged during translation by the uniform
velocity ﬁeld. To quantify the error caused by the numerical method, the normalized
L1norm is used.
E¼
P
f_{j}nDv fe
jDv
P
f_{j}iDv ð24Þ

where the superscript n designates the numerical solution, the superscript e the exact solution, and the superscript i the initial condition. The sum is over all the cells in the domain.

The L1errors for the case u¼ v ¼ 1 calculated by using rectangular grids with 5050, 100100, and 200200 cells and various sizes of time step, corresponding to Courant numbers 0.1, 0.25, 0.5, and 0.75, are shown in Figure 5 for the square and circle translations. The Courant number is deﬁned as

Cn¼

P_{maxðC}

fj;0ÞDt

Dv ð25Þ

where Cfjis the volumetric ﬂux through the face fjand the sum is taken over all the

faces of the cell. It is obvious that the error decreases as the grid size or time interval is reduced. In view of the ﬁgure, the accuracy of the calculation for the transport of the circle is generally higher than that for the square. Comparing the resulting shapes of the two cylinders given in Figure 6a (for Cn¼ 0.5) and Figure 6b (for Cn¼ 0.1)

demonstrates this point. Large errors occur at the corners of the square where the slope is discontinuous. Wavy forms appear at the two corners on the upper left and lower right. This situation can be soothed by reducing the Courant number and grid spacing. For the case of circle the interface is much smoother. However,

Figure 5. Numerical error against grid number for advection of cylinders in uniform ﬂow.

as seen from Figure 6a (Cn¼ 0.5), nonsmoothness still can be identiﬁed when the

Courant number is large. It was mentioned that the present method is applicable to grids of arbitrary geometry. Figure 6c presents the results using a grid with 22,474 triangular cells.

Advection of a Circle in a Shear Flow

In real ﬂows, the ﬂow velocity is nonuniform. The interface is subject to strain-ing and, thus, deforms continuously. To mimic this situation, the followstrain-ing velocities are assumed:

~ V

V¼ sin x cos y~ii cos x sin y~jj ð26Þ A circle of diameter 0.4p is placed in a square with side length p. The initial center of the circle is located at (p=2, (1þ p)=5). It is advected by the shear ﬂow in 16 units of time ﬁrst, followed by reversing the velocity ﬁeld for another 16 units of time. The circle will be stretched and deformed by the ﬂow straining in the for-ward step and recovers its original shape at the end of the backfor-ward step.

In the calculation, the time step is determined in the way that the maximum cell
Courant number is less than 0.1. The resulting geometries after the forward and
backward steps are given in Figure 7. The circle is gradually stretched to become
serpentine-like during the forward stage. With low grid densities the tail is shortened.
It is broken and separated form the main body for the 50_{50 grid. It can also be seen}

that the recovered circle becomes irregular for the lowest level of grid by the end of the backward stage. The prediction accuracy deteriorates by using triangular grids.

Figure 6. Resulting cylinders in uniform ﬂow (color ﬁgure available online).

The tail of the serpentine becomes shorter and the circle is not so round even for the grid with 89,776 cells. The same ﬂow problem has been tested by the present authors in [14] using the ﬂux blending schemes FBICS and CICSAM along with a 100100 rectangular grid. Comparing the results indicates that the recovered shape obtained by the present method is closer to a round circle.

Collapse of a Water Column

This case becomes a classical problem to test numerical methods for tracing a water front because experimental data are available in [41]. A water column 0.146 m in width and 0.292 m in height is initially placed at the left corner if a tank of size 0.584 m 0.340 m. The large density difference between the two phases often causes computational instability. Without a smoothing activity, calculations simply fail unless the highly diffusive, ﬁrst-order upwind difference scheme is used. Figure 8 presents the ﬂow ﬁeld obtained by using the smoothing procedure described in Section 3.1. It is clear that with only one smoothing sweep the ﬂow on the gas side Figure 7. Advection in shear ﬂow. The results in the ﬁrst row are obtained after the forward stage and those in the second row after the backward stage.

Figure 8. Flow ﬁeld for the collapse of water column: (a) with one smoothing step; (b) with two smoothing steps (color ﬁgure available online).

of the interface is unstable. The velocities vary in an irregular manner on the part of the interface where the water front is retreating and the gas ﬂows toward it. A secondary vortex appears near the left boundary. These oscillating phenomena can be suppressed when one more smoothing step is taken, as seen in Figure 8b.

Comparison of the predicted interface position with experiments [41] is made in Figure 9. The one shown in Figure 9a is the edge position of the advancing leading front along the ﬂoor, which is slightly overpredicted. This may be attributed to the uncertainties in determining the exact location of the front edge in experiments. Such an overprediction can also be found in other studies [11, 14]. The agreement in terms of the height of the collapsing column along the side wall shown in Figure 9b is good. However, as is not surprising, the curve obtained by one smoothing step exhibits oscillating behavior.

Rising of a Gas Bubble

A circular gas bubble of diameter 0.01 m is released in a chamber of size 0.05 m0.15 m. The densities of the gas and liquid phases are assumed to be 1,000 and 1 kg=m3 and the viscosities to be 5.56e-1 and 5.56e-4 Ns=m2. Three different values of surface tension coefﬁcient are under consideration: 9.79e-01, 9.79e-2, and 9.79e-3 N=m. This leads to Eotvos numbers 1, 10, and 100 and Morton numbers 0.001, 1, and 1,000. The Eotvos number is deﬁned as Eo¼ gd2

eðql qgÞ=r and the

Morton number as Mo¼ gm4

l=qlr3, where the subscripts l and g denote the liquid

and gas phases, respectively, and de is the initial bubble diameter. The evolution of

the bubble pattern at t¼ 0.05, 0.15, and 0.25 s for the three cases is shown in Figure 10. It is recognized that the interfacial tension force tend to minimize the surface area of the bubble. Therefore, the bubble pattern remains nearly a circle when the Eotvos number is sufﬁciently small. It is transformed into the shape of a circular cap as Eo is enlarged. Further increase of this number brings about a crescent shape. Figure 9. The collapse of water column: (a) position of the leading front edge along the ﬂoor; (b) height of the retreating front along the wall.

The time variation of the position of the bubble center is shown in Figure 11. The relation becomes linear after a short initial stage. The uprising velocity of the bubble, represented by the slope of the curve, is lower initially, followed by a gradual increase.

Figure 10. Evolution of the rising bubble for the three cases.

Figure 11. Center positions of rising bubbles.

Finally, it reaches a constant value in the linear stage, which is termed terminal velocity. The terminal velocities are 0.166, 0.12, and 0.115 m=s for the three cases. The gradual decrease of this speed with the Eotvos number is expected, because the deformation of the bubble results in an increase of drag.

Rayleigh-Taylor Instability

The computational domain is of dimensions 1 3. Initially, the top one-third is ﬁlled with a heavier ﬂuid of density 1.2 and the lower two-thirds with another ﬂuid of lower density 1.0. The governing equations are solved with the Froude number [Fr¼ U=(gL)1=2

] set at 0.5 and the Reynolds number (Re¼ qUL=m) at 500. Symmetric condition is speciﬁed on the left boundary, while no-slip condition on the other boundaries. A small perturbation of the form 0.02cos(1þ x)p is imposed on the initial interface between the two ﬂuids. The evolution of the unstable ﬂow is shown in Figure 12 for times t¼ 4, 6, and 8 s. In the early stage, the perturbation is enlarged in the linear sense. Then the nonlinear interaction of the descending ﬂow of the heavier layer and the ascending ﬂow of the lighter layer results in a mushroom-like structure when time reaches 4 s. The ﬂow patterns are further complicated by wrapping around the mushroom cap at t¼ 8 s. This ﬂow ﬁeld is very similar to the calculations of Rudman [9].

CONCLUSIONS

A conservative interpolation scheme for interface tracking (CISIT), within the framework of the VOF, has been introduced to overcome the difﬁculties encountered

Figure 12. Evolution of Rayleigh-Taylor ﬂow obtained using 96288 grid.

at the interface embedded in two-ﬂuid ﬂows. The formulation of the method is very simple. It is easy to code and applicable to grids of arbitrary geometry without caus-ing any further complication. The resultcaus-ing distribution of ﬂuid fraction is a good approximation to the Heaviside function. It has been seen from real-ﬂow calcula-tions that stable solucalcula-tions are obtained with two sweeps of the smoothing process.

REFERENCES

1. S. O. Unverdi and G. Tryggvason, A Front-Tracking Method for Viscous, Incompress-ible, Multi-ﬂuid Flows, J. Comput. Phys., vol. 100, pp. 25–37, 1992.

2. G. Tryggvasson, B. Bunner, A. Esmaeeli, D. Juric, N. Al-Rawahi, W. Tauber, J. Han, S. Nas, and Y.-J. Jan, A Front-Tracking Method for the Computations of Multiphase Flow, J. Comput. Phys., vol. 169, pp. 708–759, 2001.

3. S. Osher and J. A. Sethian, Fronts Propagating with Curvature-Dependent Speed: Algorithms Based on Hamilton-Jacobi Formulations, J. Comput. Phys., vol. 79, pp. 12–49, 1988.

4. M. Sussman, P. Smereka, and S. Osher, A Level Set Approach for Computing Solutions to Incompressible Two-Phase Flow, J. Comput. Phys., vol. 114, pp. 146–159, 1994. 5. M. Sussman, E. Fatami, P. Smereka, and S. Osher, An Improved Level Set Method for

Incompressible Two-Phase Flows, Comput. Fluids, vol. 27, pp. 663–680, 1998.

6. M. Sussman, and E. Fatami, An Efﬁcient, Interface-Preserving Level Set Redistancing Algorithm and Its Application to Interfacial Incompressible Fluid Flow, SIAM J. Sci. Comput., vol. 20, pp. 1165–1191, 1999.

7. R. R. Nourgaliev, S. Wiri, N. T. Dinh, and T. G. Theofanous, On Improving Mass Conservation of Level Set by Reducing Spatial Discretization Errors, Int. J. Multiphase Flow, vol. 31, pp. 1329–1336, 2005.

8. C. W. Hirt and H. D. Nichols, Volume of Fluid (VOF) Method for the Dynamics of Free Boundaries, J. Comput. Phys., vol. 39, pp. 201–225, 1981.

9. M. Rudman, Volume-Tracking Methods for Interfacial Flow Calculations, Int. J. Numer. Meth. Fluids, vol. 24, pp. 671–691, 1997.

10. S. T. Zalesak, Fully Multidimensional Flux-Corrected Transport Algorithms for Fluids, J. Comput. Phys., vol. 31, pp. 335–262, 1979.

11. S. Muzaferija, M. Peric, P. Sames, and T. Schellin, A Two-Fluid Navier-Stokes Solver to Simulate Water Entry, Proc. 22th Symp. on Naval Hydrodynamics, Washington, DC, pp. 638–649, 1998.

12. O. Ubbink and R. I. Issa, A Method for Capturing Sharp Fluid Interfaces on Arbitrary Meshes, J. Comput. Phys., vol. 153, pp. 26–50, 1999.

13. M. Darwish and F. Moukalled, Convective Schemes for Capturing Interfaces of Free-Surface Flows on Unstructured Grids, Numer. Heat Transfer B, vol. 49, pp. 19–42, 2006. 14. Y.-Y. Tsui, S.-W. Lin, T.-T. Cheng, and T.-C. Wu, Flux-Blending Schemes for Interface Capture in Two-Fluid Flows, Int. J. Heat Mass Transfer, vol. 52, pp. 5547–5556, 2009.

15. P. H. Gaskell and A. K. C. Lau, Curvature-Compensated Convective Transport: SMART, A New Boundedness-Preserving Transport Algorithm, Int. J. Numer. Meth. Fluids, vol. 8, pp. 617–641, 1988.

16. W. F. Noh and P. Woodward, SLIC (Simple Line Interface Calculation), Lecture Notes Phys., vol. 59, pp. 330–340, 1976.

17. D. L. Youngs, Time-Dependent Multi-material Flow with Large Fluid Distortion, in K. W. Morton and M. J. Baines (eds.), Numerical Methods for Fluid Dynamics, vol. 24, pp. 273–285, Academic Press, New York, 1982.

18. W. J. Rider, and D. B. Kothe, Reconstructing Volume Tracking, J. Comput. Phys., vol. 141, pp. 112–152, 1998.

19. J. E. Pilliod, Jr., and E. G. Puckett, Second-Order Accurate Volume-of-Fluid Algorithms for Tracking Material Interfaces, J. Comput. Phys., vol. 199, pp. 465–502, 2004. 20. N. Ashgriz and J. Y. Poo, FLAIR: Flux Line-Segment Model for Advection and Interface

Reconstruction, J. Comput. Phys., vol. 93, pp. 449–468, 1991.

21. F. Mashayek and N. Ashgriz, Advection of Axisymmetric Interfaces by the Volume-of-Fluid Method, Int. J. Numer. Meth. Volume-of-Fluids, vol. 20, pp. 1337–1361, 1995.

22. F. Mashayek and N. Ashgriz, A Hybrid Finite-Element-Volume-of-Fluid Method for Simulating Free Surface Flows and Interfaces, Int. J. Numer. Meth. Fluids, vol. 20, pp. 1363–1380, 1995.

23. M. Sussman and E. G. Puckett, A Coupled Level Set and Volume-of-Fluid Method for Computing 3D and Axisymmetric Incompressible Two-Phase Flows, J. Comput. Phys., vol. 162, pp. 301–337, 2000.

24. G. Son and N. Hur, A Coupled Level Set and Volume-of-Fluid Method for the Buoyancy-Driven Motion of Fluid Particles, Numer. Heat Transfer B, vol. 42, pp. 523–542, 2002. 25. S. P. van der Pijl, A. Segal, C. Vuik, and P. Wesseling, A Mass-Conserving

Level-Set Method for Modelling of Multi-Phase Flows, Int. J. Numer. Meth. Fluids, vol. 47, pp. 339–361, 2005.

26. X. Yang, A. J. James, J. Lowengrub, X. Zheng, and V. Cristini, An Adaptive Coupled Level-Set=Volume-of-Fluid Interface Capturing Method for Unstructured Triangular Grids, J. Comput. Phys., vol. 217, pp. 364–394, 2006.

27. D. L. Sun and W. Q. Tao, A Coupled Volume-of-Fluid and Level Set (VOSET) Method for Computing Incompressible Two-Phase Flows, Int. J. Heat Mass Transfer, vol. 53, pp. 645–655, 2010.

28. R. Tavakoli, R. Babaei, N. Varahram, and P. Davami, Numerical Simulation of Liquid=Gas Phase Flow during Mold Filling, Comput. Meth. Appl. Mech. Eng., vol. 196, pp. 697–713, 2006. 29. S. J. Mosso, B. K. Swartz, D. B. Kothe, and S. P. Clancy, Recent Enhancement of Volume Tracking Algorithm for Irregular Grids, Los Alamos Natl. Lab. Rep. LA-CP-96-227, Los Alamos, NM, 1996.

30. K. Shahbazi, M. Paraschivoiu, and J. Mostaghimi, Second Order Accurate Volume Tracking Based on Remapping for Triangular Meshes, J. Comput. Phys., vol. 188, pp. 100–122, 2003.

31. N. Ashgriz, T. Barbat, and G. Wang, A Computational Lagrangian-Eulerian Advection Remap for Free Surface Flows, Int. J. Numer. Meth. Fluids, vol. 44, pp. 1–32, 2004. 32. J.-T. Yeh, Simulation and Industrial Application of Inkjet, 7th Natl. CFD Conf.,

Pingtung, Taiwan, 2000.

33. J.-T. Yeh, A VOF-FEM, and Coupled Inkjet Simulation, Proc. ASME FEDSM, New Orleans, LA, pp. 1–5, 2001.

34. J. U. Brackbill, D. B. Kothe, and C. Zemach, A Continuum Method for Modeling Surface Tension, J. Comput. Phys., vol. 100, pp. 335–354, 1992.

35. C. S. Peskin, Numerical Analysis of Blood Flow in the Heart, J. Comput. Phys., vol. 25, pp. 220–252, 1977.

36. J. J. Monaghan, Simulating Free Surface Flow with SPH, J. Comput. Phys., vol. 110, pp. 399–406, 1994.

37. Y.-Y. Tsui and T.-C. Wu, A Pressure-Based Unstructured-Grid Algorithm Using High-Resolution Schemes for All-Speed Flows, Numer. Heat Transfer B, vol. 53, pp. 75–96, 2008.

38. Y.-Y. Tsui and T.-C. Wu, Use of Characteristic-Based Flux Limiters in a Pressure-Based Unstructured-Grid Algorithm Incorporating High-Resolution Schemes, Numer. Heat Transfer B, vol. 55, pp. 14–34, 2009.

39. R. I. Issa, Solution of the Implicitly Discretised Fluid Flow Equations by Operator-Splitting, J. Comput. Phys., vol. 62, pp. 40–65, 1986.

40. Y.-Y. Tsui and Y.-F. Pan, A Pressure-Correction Method for Incompressible Flows Using Unstructured Meshes, Numer. Heat Transfer B, vol. 49, pp. 43–65, 2006.

41. J. C. Martin and W. J. Moyce, An Experimental Study of the Collapse of Liquid Columns on a Rigid Horizontal Plane, Phil. Trans. Roy. Soc. Lond., Ser. A, vol. 244, pp. 312–324, 1952.