• 沒有找到結果。

A VOF-Based Conservative Interpolation Scheme for Interface Tracking (CISIT) of Two-Fluid Flows

N/A
N/A
Protected

Academic year: 2021

Share "A VOF-Based Conservative Interpolation Scheme for Interface Tracking (CISIT) of Two-Fluid Flows"

Copied!
22
0
0

加載中.... (立即查看全文)

全文

(1)

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

(2)

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-fluid (VOF) formulation and applicable to unstructured grids with arbitrary topology, is developed to track the interface in free-surface flows. 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 overfilled (f > 1) or underfilled (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 fluid 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 fluid 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 flows, reveal that this method is accurate and robust.

INTRODUCTION

Multifluid or multiphase flow features an interface separating fluids with differ-ent fluid properties. Accurate prediction of this interfacial flow presdiffer-ents a challenging task in numerical simulation because there exists a jump in fluid density and viscosity across the interface. This jump phenomenon can be identified as a Heaviside (step) function in theory. Different from the high-speed compressible flow, in which the flow 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

(3)

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

The geometric shape of the interface continues to change during fluid flow. 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 (fixed grids) with a method to track the movement of the interface. Tryggvason and co-workers [1, 2] developed a front-tracking method to fulfill this job. A fixed grid is assigned to the velocity field and a grid of lower dimen-sion to the fluid 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 fluid 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 fluids 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 defined 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 fluid flux Cn Courant number Eo Eotvos number f volume fraction F flux 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 flow velocity vector

w weighting factor c(r) flux limiter function ~dd distance vector Dt time-step size Dv cell volume m fluid viscosity q fluid density

r surface tension coefficient

sij viscous stress

/ Cartesian velocity components U level set function Subscripts

D downstream cell f cell face fi 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

(4)

qU

qt þ ~VV rU ¼ 0 ð1Þ where ~VV is the velocity of the flow. 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-fluid (VOF) function, which is defined as the fraction of volume in the mesh cell occupied by one of the two fluids. The VOF function f is either unity or zero in the cells con-taining a single fluid 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 difficult 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 classified into two cate-gories. One is to solve the VOF in a way similar to the shock capturing for super-sonic flows. 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 first-order upwind scheme is most stable but causes large diffusion, while the first-order down-wind scheme has negative artificial 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 first to combine the upwind and downwind schemes in such a way as to minimize diffusion in interfacial flow calculations. Rudman [9] adopted a method based on the flux-corrected transport algorithm [10]. It consists of two stages. In the first 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 first 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],

(5)

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 fluid fluxes 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 diffi-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 configurations that need to be considered for determination of fluid fluxes through the cell faces. Different from the PLIC, in which the interface is fitted 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 fluid flux 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 flow 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 flow 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 difficult 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 fluid 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 fluxing across cell faces, it is necessary to consider 16 interface configurations for 2-D flows and 64 configura-tions for 3-D flows, though these numbers of configuraconfigura-tions can be reduced after reg-ulating the cell geometry [28]. The extension of the method for unstructured meshes

(6)

with irregular cell geometry is not easy, partly due to the increased difficulty in reconstructing the interface and also due to the increased complexity in estimating the flow fluxes. 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 fluid 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 fulfilled in a prediction-and-correction way. In the following, details of this solution method are described first. This is followed by a brief description of the method for the velocity field. Then, fidelity of the algorithm and verification 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 fluid flux, the VOF is sought on grid nodes where velocities are located. The VOF on the grid node is first 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 finite-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 find 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

(7)

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 finite-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 flux 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 fluxes 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 flux 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 fluid flux, and ffis the fluxing 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 fluxing value ffusing the VOFs

of neighboring cells [11–14].

In the present interface-reconstruction method, the face flux is the volumetric flux through the portion of the face wetted by the fluid:

(8)

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 fluid. 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 filled with fluid (color figure available online).

(9)

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 identified from Figure 2a, in which the inter-face is advancing to the right in a 1-D flow. In the beginning, the VOF value in the cell increases with the advancement of the interface, because fluid flows into the cell from the left face and no fluid 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 fluidDf(¼ fP 1)must be reallocated to the downstream cell N1. For a

multi-dimensional flow, 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 efflux through the corresponding face to the total efflux:

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

where Cfi ð¼~VVfi ~SSfiÞ is the volumetric fluid flux 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 figure available online).

(10)

After the neighboring cells are adjusted, the VOF for this overfilling 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 flow, the VOF becomes less than zero at the time of crossing the right face. The overdepleted volume of fluid must be retrieved from the downstream cell N1. Similar to the overfilling 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 flow has large velocity gradients in the

trans-verse direction. In such flows, different from overfilling, 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), fluid must

be retrieved from the downstream cells to fill 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

flows is that the VOF may remain greater than zero when the interface retreats to leave the cell (Figure 3b). Similar to the underfilling cell, the volume of fluid remains unchanged and the fluid 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), fluid 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 figure available online).

(11)

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 first using the VOF from the last time step. 2. New values of VOF are sought via conserving fluid mass in interface cells. 3. The VOFs are adjusted to ensure that no cells are overfilled, underfilled,

overdepleted, or underdepleted. It is noted that after the correction process for underfilling and underdepleting cells, it is possible for some of the neighboring cells to be overfilled 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 first, followed by solving for the velocity field. The flow in both fluids is assumed to be laminar and incom-pressible. The conservation of mass and momentum for the two fluids 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 flows. 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 coefficient.

Smoothing of Fluid Properties

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

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

(12)

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

where the subscripts 1 and 2 denote the two fluids. When the difference between the fluid properties, especially the density, of the two fluids 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 artificially. 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 finite-volume method. The convective momentum flux is expressed as

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

where / designates the Cartesian velocity components. The face flux 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 flux limiter, which depends on the gradient ratio r defined by

r¼/U  /UU /D /U

ð20Þ

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

(13)

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 flux, suitable for unstructured grids, is approximated by

Fd ¼ mS 2 f ~dd PN ~SSf ð/N /PÞ þ m r/ð Þf ~SSf  S2 f ~dd PN ~SSf ~dd 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 fields. A momentum interpolation method is then employed to find the velocities on the cell face and the mass flux through the face. A pressure-correction equation can be obtained by forcing the corrected mass flux 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 satisfied. The way to derive the face mass flux and the pressure-correction equation can be found in [40].

RESULTS AND DISCUSSION

The solution method for the VOF is first tested in two model problems with given velocities. One is the advection of a square or a circular cylinder in a uniform velocity field and the other is the advection of a circular cylinder in a shear flow. Real flows are then considered. The realistic flows 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

(14)

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 field. To quantify the error caused by the numerical method, the normalized L1norm is used. E¼ P fjnDv fe jDv     P fjiDv ð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 defined as

Cn¼

PmaxðC

fj;0ÞDt

Dv ð25Þ

where Cfjis the volumetric flux 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 figure, 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 flow.

(15)

as seen from Figure 6a (Cn¼ 0.5), nonsmoothness still can be identified 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 flows, the flow 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 flow in 16 units of time first, followed by reversing the velocity field for another 16 units of time. The circle will be stretched and deformed by the flow 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 5050 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 flow (color figure available online).

(16)

The tail of the serpentine becomes shorter and the circle is not so round even for the grid with 89,776 cells. The same flow problem has been tested by the present authors in [14] using the flux 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, first-order upwind difference scheme is used. Figure 8 presents the flow field obtained by using the smoothing procedure described in Section 3.1. It is clear that with only one smoothing sweep the flow on the gas side Figure 7. Advection in shear flow. The results in the first row are obtained after the forward stage and those in the second row after the backward stage.

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

(17)

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 flows 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 floor, 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 coefficient 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 defined 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 sufficiently 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 floor; (b) height of the retreating front along the wall.

(18)

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.

(19)

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 filled with a heavier fluid of density 1.2 and the lower two-thirds with another fluid 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 specified 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 fluids. The evolution of the unstable flow 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 flow of the heavier layer and the ascending flow of the lighter layer results in a mushroom-like structure when time reaches 4 s. The flow patterns are further complicated by wrapping around the mushroom cap at t¼ 8 s. This flow field 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 difficulties encountered

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

(20)

at the interface embedded in two-fluid flows. 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 fluid fraction is a good approximation to the Heaviside function. It has been seen from real-flow 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-fluid 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 Efficient, 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.

(21)

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.

(22)

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.

數據

Figure 1. A control volume partially filled with fluid (color figure available online).
Figure 2. Advancing of interface through the face of a control volume (color figure available online).
Figure 3. Retreating of interface through the face of a control volume (color figure available online).
Figure 4. Illustration for determination of the far upstream node UU.
+7

參考文獻

相關文件

When the relative phases of the state of a quantum system are known, the system can be represented as a coherent superposition (as in (1.2)), called a pure state; when the sys-

The Shannon entropy for a specific source X can be seen as the amount of our ignorance about the value of the next letter, or the amount of indeterminancy of the unknownm letter..

 Machine language ( = instruction set) can be viewed as a programmer- oriented abstraction of the hardware platform.  The hardware platform can be viewed as a physical means

 Machine language ( = instruction set) can be viewed as a programmer- oriented abstraction of the hardware platform.  The hardware platform can be viewed as a physical means

Have shown results in 1 , 2 &amp; 3 D to demonstrate feasibility of method for inviscid compressible flow problems. Department of Applied Mathematics, Ta-Tung University, April 23,

Based on [BL], by checking the strong pseudoconvexity and the transmission conditions in a neighborhood of a fixed point at the interface, we can derive a Car- leman estimate for

At least one can show that such operators  has real eigenvalues for W 0 .   Æ OK. we  did it... For the Virasoro

• When a number can not be represented exactly with the fixed finite number of digits in a computer, a near-by floating-point number is chosen for approximate