Flux-blending schemes for interface capture in two-fluid flows
Yeng-Yung Tsui
*, Shi-Wen Lin, Tong-Ting Cheng, Tian-Cherng Wu
Department of Mechanical Engineering, National Chiao Tung University, Hsinchu 300, Taiwan, ROC
a r t i c l e
i n f o
Article history: Received 4 March 2009
Received in revised form 17 June 2009 Accepted 17 June 2009
Available online 31 July 2009 Keywords:
Volume of fluid (VOF) method Two-fluid flows
Free surface flows Interface capturing Unstructured meshes
a b s t r a c t
This paper is aimed at developing a robust algorithm (FBICS) to capture the interface between two immis-cible fluids without the need of interface reconstruction. The advection equation of the volume fraction is solved using the fully conservative finite volume method. Determination of the convective flux through each cell face is based on blending of high resolution schemes and compressive schemes to preserve the sharpness and boundedness of the interface. The flux-blending practice is fulfilled with the use of flux limiters. Test on simple advection flow problems indicates that the well-known CICSAM and HRIC schemes lose accuracy as the Courant number increases. In contrast, the present method maintains high-accuracy performance for Courant numbers up to one. The capability of the method to cope with the complicated dynamics of free surface flows is demonstrated via calculation of the collapsing flow of a water column with an obstacle.
Ó 2009 Elsevier Ltd. All rights reserved.
1. Introduction
The flow involving two immiscible fluids with a free interface has been of interest to many engineering applications, including biochemical engineering, marine engineering, and casting, injection or extrusion processes. In numerical simulation of this kind of flow, an important issue is to track the motion of the interface. Meanwhile, the sharpness of the interface must be maintained. The methods developed in the past for free surface flows can be divided into two categories: Lagrangian and Eulerian. In Lagrangian methods, the flow field of the considered fluid is covered by a mesh moving with the fluid. The fluid boundaries always coincide with the grid bound-aries and the fluid inside each cell of the grid always remains in that cell[1,2]. This method is not suited for flows undergoing large distor-tions because the mesh will be greatly deformed, which degrades the accuracy of the solution and causes instability of the solution procedure. To soothe these problems, Muzaferija and Peric [3] allowed the grid lines inside the flow field to move in an Eulerian– Lagrangian manner. However, it still can not cope with breaking or overturning of the interface.
In the category of Eulerian methods, the grids used for fluid flow calculations are fixed without motion. The interface is treated as a sharp front moving through the grid. Two basic approaches were adopted to track the interface. In the front-tracking method, the interface is represented by a connected set of points, which forms a moving boundary, whereas a stationary grid is constructed for the solution of the Navier–Stokes equations[4,5]. This method is rather difficult to implement due to the interaction between the
interface grid and the Eulerian grid. The complexity is further enhanced by the necessity of remeshing the interface grid during time marching. Another problem results from the interaction of a front with another front when both appear in a grid cell simultaneously.
A different front-tracking approach is the level set method[6,7]. A continuous function (the level set function) is defined as the shortest distance between the considered point and the interface. The interface is, thus, located at a level set value of zero. The func-tion is transported by the fluid via solving an advecfunc-tion equafunc-tion, similar to the VOF method addressed below. One of the disadvan-tages is that it suffers numerical errors of interface smearing. The main weakness of this method is loss of accuracy in highly dis-torted flows because mass conservation is not guaranteed[8].
An alternative to front-tracking is tracking. In volume-tracking methods, different fluids must be marked by different indicators. In the marked-and-cell (MAC) method of Harlow and Welch [9], the considered fluid is covered by massless particles. These marked particles are advected with the local fluid velocity. Their distribution determines the fluid configuration. Because the number density of particles is usually low, quantitative informa-tion on interface orientainforma-tion is poor. This problem becomes even severe in regions subject to high shear flow. Thus, a large number of marker particles are needed, resulting in significant computa-tional overheads, especially for three-dimensional flows. It needs to be addressed here that in MAC, the fluid flow velocities required for the transportation of particles are determined by solving the equations in a prescribed mesh. Particle methods without need of mesh were developed by Monaghan[10]and Koshizuka et al. [11]. A major disadvantage of these particle methods is that inflow and outflow through boundaries can not be handled. To tackle this 0017-9310/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved.
doi:10.1016/j.ijheatmasstransfer.2009.06.026
* Corresponding author. Tel.: +03 5131 556; fax: +03 5720 634. E-mail address:[email protected](Y.-Y. Tsui).
Contents lists available atScienceDirect
International Journal of Heat and Mass Transfer
j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / i j h m tproblem, a hybrid scheme, consisting of Lagrangian and Eulerian phases, was proposed by Yoon et al. [12]. After the Lagrangian calculation, a one-dimension local grid is generated for convection interpolation.
Another way to mark the fluid is the use of an indicator func-tion, known as the volume fraction. This function represents the fraction of a local cell volume occupied by one of the fluids. The transport of the indicator function through the domain is governed by a hyperbolic equation. This approach is generally referred to as the VOF method. In the past, a variety of techniques have been developed to maintain a well defined interface within the frame of the VOF method. In one class of such schemes, the exact location of the front is discarded and the interface is reconstructed in an approximate manner using information of local fluid volumes. The SLIC (simple line interface calculation) of Noh and Woodard [13]approximate the interface in each cell as a vertical or horizon-tal line, depending on the sweep direction, in an operator-splitting scheme. In the PLIC (piecewise linear interface calculation) method of Youngs[14], an oblique line is used in each cell, which gives a better approximation to the interface and leads to higher accuracy. After the interface is approximately reconstructed, new fluid frac-tion values are computed by moving the fluid volume according to the local velocity field[15,16].
In application of the above VOF methods, the computational cells adopted are usually rectangular mainly due to the difficulty of interface reconstruction for irregular meshes. Another approach, requiring no explicit interface reconstruction, is to capture the sharp interface by directly solving the hyperbolic equation of the indicator function. As in the shock-capturing calculations for high-speed compressible flows, it needs to overcome the problem of numerical diffusion, which smoothes out the sharp gradient of the interface, and that of numerical dispersion, which may cause the volume fraction unbounded. A scheme based on the flux-cor-rected transport (FCT) concept has been proposed by Rudman [17,18]for interfacial flow calculations. The applied FCT schemes were developed mainly based on one-dimensional theory, which renders them not suitable for unstructured grids.
Other choices are the high-resolution schemes of TVD (total variation diminishing) and NVD (normalized variable diagram). However, direct use of these schemes cannot eliminate the
numer-ical diffusion and unboundedness of the interface. An effective way to alleviate these problems is to blend high-resolution schemes and compressive schemes together. These blending schemes switch in a continuous manner from the high-resolution scheme to the compressive scheme, depending on the angle between the interface and the grid lines. Several such schemes have been intro-duced, such as the CICSAM (compressive interface capturing scheme for arbitrary meshes) of Ubbink and Issa[19], the HRIC (high resolution interface capturing) of Muzaferija and Peric[20], and the STACS (switching technique for advection and capturing of surfaces) of Darwish and Moukalled[21].
In recent years, a pressure-based algorithm within the frame-work of unstructured mesh was developed for incompressible flows by the group of the present authors [22]. It was further extended to deal with all-speed flows, including incompressible, subsonic, transonic and supersonic flows[23,24]. To capture the shock waves embedded in high-speed compressible flows without causing smear and unboundedness, the amount of the convective flux through cell faces is controlled with the use of flux limiters. A variety of TVD and NVD schemes can easily be implemented in the form of the limiting function. Based on the unstructured-grid algorithm and the limiter technique, a methodology is introduced in this study to capture the sharp front of the free-surface flows. 2. Mathematical model
The different fluid flows separated by an interface are assumed to be incompressible and obey the same form of equations describ-ing the conservation of mass and momentum:
@Vj @xj ¼ 0 ð1Þ @
q
Vi @t þ @ @xj ðq
VjViÞ ¼ @P @xi þ@s
ij @xj þq
gi ð2Þwhere Vjis the velocity,
q
the density, P the pressure,s
ijthe viscousstress tensor, and githe gravitational acceleration. The mixture of
fluids is considered as a single continuum. The fluid properties can be expressed as a function of the volume fraction f. The density and viscosity at a location are evaluated by
Nomenclature
CN Courant number
Fc convective flux of momentum
Fd diffusive flux of momentum
g gravitational acceleration _
m mass flux
P pressure
P0 boundary pressure at the node next to boundaries
Pb pressure on the boundary face
r gradient ratio
~S surface vector
t time
Dt time step size ~
V; Vj velocity vector
DV volume of the cell xj Cartesian coordinates
Greek symbols
~d distance vector
/ Cartesian velocity components
c
(r) flux limiter depending on the gradient ratiol
viscosityh angle between grid line and fluid interface
q
densitys
ij viscous stressesx
(h) weighting factor depends on the angle h Subscripts b boundary C neighboring node D downstream node f control surface P principal node U upstream nodeUU far upstream node
Superscripts
n new time level
o old time level
HR high-resolution scheme
BD bounded downwind scheme
q
¼ fq
1þ ð1 f Þq
2 ð3aÞl
¼ fl
1þ ð1 f Þl
2 ð3bÞwhere the subscripts 1 and 2 denote the two fluids.
The volume fraction f is advected as a Lagrangian invariant and has zero material derivative:
Df Dt¼ @f @tþ Vj @f @xj¼ 0 ð4Þ
With the help of the continuity equation, the advection equation can be cast into the divergence form:
@f @tþ
@ @xj
ðVjf Þ ¼ 0 ð5Þ
3. Flux-blending interface-capturing scheme (FBICS)
For discretization the differential equations are integrated over a control volume. Those terms in divergence form are transformed into surface integrals according to the divergence theorem. The convective flux through each face of the control volume can be approximated using the mean value theorem. Thus, Eq. (5) for the volume fraction can be approximately expressed as
D
VD
tðf n P fPoÞ þ X f Ffff ¼ 0 ð6ÞHere the superscripts n and o denote the new and old time values, the subscripts P and f designate the centroids of the control volume and the surrounding faces, and Ffis the volumetric flux through a
face defined by
Ff ¼ ~Vf ~Sf ð7Þ
where ~Sfis the surface vector of the face (seeFig. 1). The sum is over
all the surrounding faces of the control volume.
To complete discretization, the volume fraction at each face re-quires estimate from neighboring nodes:
ff¼ fUþ
c
ðrÞ2 ðfD fUÞ ð8Þ
where, as seen inFig. 2a, the subscript U denotes the node upstream of the considered face and the subscript D is the one downstream. The function
c
(r) is the flux limiter depending on the gradient ratio: r ¼fU fUUfD fU
ð9Þ
where fUUis the value at a node far upstream. This far upstream node
is easily identified in structured grids (seeFig. 2a). The determination of the value at this node in unstructured grids is referred later to Eq. (20). It is obvious that the expression represents the upwind differ-ence scheme for
c
= 0, the central difference scheme forc
= 1 and the downwind difference scheme forc
= 2. Thus the second term on the right of Eq.(8)represents anti-diffusion to the upwind scheme.In solving the advection equation, the upwind part of the con-vection flow is treated by the Crank–Nicolson scheme whereas the anti-diffusion part in the explicit manner. Thus, it is obtained
D
VD
tðf n p f o pÞ þ X f Ff 1 2ðf n Uþ f o UÞ þc
ðrÞ 2 ðf o D f o UÞ ¼ 0 ð10ÞThe limiters for linear difference schemes can be expressed as sim-ple linear functions of the gradient ratio. For examsim-ples,
c
(r) = r is the linear upwind scheme,c
(r) = r/4 + 3/4 the QUICK scheme, andc
(r) = r/2 + 1/2 the Fromm’s scheme. Nonlinear schemes like TVD and NVD schemes can also be easily implemented by expressing the limiting functions in terms of r[23]. Tests of a variety of such schemes on model problems had shown that the smear of interface is serious even with the use of the SUPERBEE[25]which was re-garded as the most non-diffusive among the TVD schemes.It was recognized that in contrast to the upwind differencing, the downwind scheme could cause significant oscillations and, thus, instability of the solution due to the fact the sign of its artifi-cial viscosity is negative. This anti-diffusion feature makes it capa-ble to compress a smooth function into a step profile[26]. To take advantage of this compression nature, a strategy is to blend slightly diffusive high-resolution schemes (HR) with bounded downwind schemes (BD):
c
¼ ½1x
ðhÞc
HRþ
x
ðhÞc
BDð11Þ This approach is referred to as Flux-Blending Interface-Capturing Scheme (FBICS). The flux limiters for the high-resolution scheme and the bounded downwind scheme used in the present article are given by
Fig. 1. A typical control volume and its neighboring cells.
Fig. 2. Illustration of the downwind, upwind and far upwind nodes for (a) structured grid and (b) unstructured grid.
c
HR¼ max 0; min 4r;1 2r þ 1 2;2 ð12aÞc
BD ¼ max½0; minð4r; 2Þ ð12bÞThe Fromm’s scheme
c
(r) = r/2 + 1/2 is employed as the basic scheme in the HR.The weighting factor
x
depends on the angle h between the grid lines and the interface. A rule of thumb is that the bounded down-wind scheme is used when the interface is parallel with the consid-ered face and it is gradually switched to the high-resolution scheme until the interface becomes orthogonal to the face. In the pre-vious studies[19–21], schemes in the NVD form were blended in the manner similar to Eq.(11)with various expressions forx
. The differ-ent functions have been tested using the presdiffer-ent blending scheme. It was found that the results are not sensitive to the choices ofx
. In the following, the formx
= cos4h[21]is employed, where h is given byh¼ cos1
r
ff ~dUDj
r
ffjj~dUDjð13Þ Here ~dUDis the distance vector directing from the upwind node U to
the downwind node D (seeFig. 2b).
The flux limiters given in Eqs.(12a) and (12b)represent combi-nations of different linear schemes, which can be more easily understood using the normalized variables formulation (NVF). ~
f ¼ f fUU
fD fUU ð14Þ
Eq.(8)can then be rewritten as ~ ff ¼ ~fUþ
c
ðrÞ 2 ð1 ~fUÞ ð15Þ where r ¼ ~ fU 1 ~fU ð16Þ It can be detected that in the NVF form, the face value ~ffdepends onthe upwind value ~fUonly.
NVD schemes must satisfy the convection boundedness criteria (CBC) [27]which can be represented by the combination of the upper triangle bounded by the lines ~fU¼ 0; ~ff¼ 1, and the diagonal
line ~ff¼ ~fUinFig. 3.
The constraints on TVD schemes are more stringent, which is bounded by an oblique line ~ff¼ 2~fU instead of the vertical line
~
fU¼ 0 in the NVD schemes as shown inFig. 3.
The aabove two schemes can then be reinterpreted using the NVF as
HR scheme ~ff ¼ 3~fU 0 < ~fU18 ~ fUþ14 18< ~fU34 ðFDSÞ 1 3 4< ~fU 1 ðDDSÞ ~ fU ~fU 0 or ~fU 1 ðUDSÞ 8 > > > > > < > > > > > : ð17aÞ BD scheme ~ff ¼ 3~fU 0 < ~fU13 1 1 3< ~fU 1 ðDDSÞ ~ fU ~fU 0 or ~fU 1 ðUDSÞ 8 > > < > > : ð17bÞ
In the above, FDS, DDS and UDS denote the Fromm’s difference scheme, downwind difference scheme and upwind difference scheme, respectively. These schemes are drawn on the normalized variable diagram inFig. 4a. Note that for schemes in the NVD range, but out of the TVD regime, which is bounded by the line ~ff¼ 2~fUas
shown in Fig. 3, are less diffusive. It is readily seen that both schemes represented by Eqs.(17a) and (17b)fall out of the TVD range for small values of ~fU because of the higher slope ~ff ¼ 3~fU,
equivalent to
c
= 4r in the flux limiter.It is emphasized that in the sketch of NVD, the diagonal repre-sents the upwind scheme and the upper bound reprerepre-sents the downwind scheme. Therefore, scheme curves closer to the
diagonal line have the character of higher diffusion and those clo-ser to the upper line and the vertical line ~fU¼ 0 are more
compres-sive. The line ~ff¼ ~fU=CNis adopted as the left boundary in the NVD
diagram in the CICSAM scheme of Ubbink and Issa[19], where CNis
the local Courant number in the considered cell. When the Courant number approaches zero, it becomes the vertical line ~fU¼ 0 and
the CICSAM scheme is compressive. However, for CNclose to one
it degenerates into the diagonal line and results in significant diffusion.
To allow the present scheme to be more compressive for small Courant numbers, the above schemes shown in Eqs. (12a) and (12b)are modified as
c
HR¼ max 0; min max 2 1CN 1 r; 4r ;1 2r þ 1 2;2 ð18aÞ
c
BD¼ max 0; min max 2 1CN 1 r; 4r ;2 ð18bÞ In the NVD form, they become
HR schemef ~ff¼ max 1 CN;3 ~ fU 0 < ~fU min 4ð1CCN NÞ; 1 8 ~ fUþ14 min CN 4ð1CNÞ; 1 8 < ~fU34 ðFDSÞ 1 3 4< ~fU 1 ðDDSÞ ~ fU ~fU 0 or ~fU 1 ðUDSÞ 8 > > > > > > < > > > > > > : ð19aÞ BD scheme ~ff¼ max 1 CN;3 ~ fU 0 < ~fU min CN;13 1 min CN;13 < ~fU 1 ðDDSÞ ~ fU ~fU 0 or ~fU 1 ðUDSÞ 8 > > > < > > > : ð19bÞ
In the above, the constraint ~ff ¼ 3~fU is utilized when the Courant
number is greater than 1/3 and replaced by ~ff ¼ ~fU=CNfor CN< 1/3
(seeFig. 4b).
For non-rectangular meshes used in unstructured-grid methods the far upstream value fUUis not available. An estimate can be
ob-tained with a pseudo node UU located at a distance of 2~dUDaway
from the node D (seeFig. 2b). The solution around the cell U is as-sumed to be distributed linearly. Thus, it can be estimated as
fUU¼ fD 2ð
r
/ÞU ~dUD ð20ÞThe volume faction obtained from the above linear extrapolation may exceed the bounds 0 and 1. To ensure boundedness, the follow-ing step is carried out[19]
fUU¼ max½minðfUU;1Þ; 0 ð21Þ
4. Solution method
Similar to the volume fraction equation, the convective flux of momentum through a face can be expressed as
Fc
¼
q
fFf/f ð22Þwhere / designates the Cartesian velocity components and, similar to the above, the face value is approximated by
/f ¼ /Uþ
c
ðrÞ2 ð/D /UÞ ð23Þ
The Van Leer scheme is used in the momentum calculation, which is given by
c
ðrÞ ¼r þ jrjr þ 1 ð24Þ
In NVD form, it can be written as ~ /f ¼ 2~/U ~/2U 0 ~/U 1 ~ /U /~U<0 or ~/U>1 ( ð25Þ
It is seen that unlike most other high-resolution schemes, this scheme features a quadratic form and is, thus, continuously differentiable.
An approximation to the diffusive flux of momentum through a face of a control volume, applicable to the unstructured grid arrangement, is given by Fd ¼
l
S 2 f ~d PC ~Sf ð/C /PÞ þl
ðr
/Þf ~Sf S2 f ~d PC ~Sf ~d PC ! ð26Þ where the subscripts P and C denote the considered control volume and the neighboring cell, ~dPCis the vector connecting nodes P and C(seeFig. 1), and ðr/Þf represents the gradients at the face obtained
by linear approximation from the gradients at nodes P and C. The solution procedure is to solve the volume fraction equation first to advance the interface. After the volume fraction is updated, velocities are predicted by solving the momentum equation using prevailing pressure. However, this velocity field does not satisfy the continuity constraint and the pressure has not been upgraded. A correction step is then conducted to adjust these variables. In this step, a pressure-correction equation is derived by forcing the corrected velocities to satisfy the mass conservation. Details about this procedure are referred to reference[22]. To make the pressure field get rid of the mass residual left by the predictor step and to obtain better approximation to the momentum conservation, a fur-ther correction to the velocities and pressure is performed, as in the PISO algorithm[28]. This completes the calculation in one time step and the solution procedure is advanced to the next time step. It is common in interfacial flows to have open boundaries on which pressures are specified. Mass flux through this boundary must be derived using the prescribed pressure. In general, there are two approaches to estimate the mass flux. One is to make an approxima-tion to the momentum equaapproxima-tion in a manner similar to the momen-tum interpolation method for internal faces[29]. However, with this method the mass is not conserved unless iteration on the solution is undertaken. Therefore, it is not appropriate in the non-iterative pro-cedure of the present method. An alterative is to make use of the principle of mass conservation.Fig. 5illustrates a control volume next to an open boundary. The boundary pressure P0is prescribed
at the centroid P of this boundary cell. The pressure on the boundary face Pbis calculated using the relation
Pb¼ P0þ
r
PP ~dPb ð27Þwhere the pressure gradient at node P is approximated by
r
PP¼ 1D
V X f Pf~Sf ð28ÞCombining the above two equations yields Pb¼
P0þ ðPf –bPf~Sf ~dPbÞ=
D
V1 ð~Sb ~dPbÞ=
D
Vð29Þ
where the summation is taken over all the faces of the cell except the boundary face. With this boundary face pressure Pb, the velocity
at the node P can be calculated in the momentum predictor step in the same way as the other internal nodes. After the mass fluxes through all the internal faces are calculated using the momentum interpolation method, the mass flux _mbthrough the boundary face
is obtained according to the mass conservation. _ mbþ X f –b _ mf ¼ 0 ð30Þ
where the summation is taken over all the internal faces apart from the boundary face.
The velocities and the volume fraction need to be determined on the open boundary. The convective boundary condition[30]is imposed, in which a wave equation is solved.
@/ @tþ uc
@/
@n¼ 0 ð31Þ
where n designates the normal direction and ucis the local flow
velocity at the boundary. After the velocities on the boundary are obtained, they are adjusted in a way that the mass flux through the boundary face is consistent with that calculated from Eq.(30).
5. Results and discussion
The algorithm described above is tested via comparison against analytical solutions and experimental measurements if available. The cases selected for testing include (1) advection of hollow cylin-ders in a uniform velocity field, (2) advection of a circle in a shear flow, (3) collapse of a water column with or without an obstacle. In the first two cases, the velocities are prescribed in advance and no coupling exists between the volume fraction field and the velocity field. They provide benchmark testing for evaluation of different schemes.
5.1. Advection of hollow cylinders in a uniform velocity field The velocities are assumed to be uniformly distributed with ~V ¼ 2~iþ~j. The computational domain is a 4 4 square. Two differ-ent cylinder shapes are under consideration. One is a hollow square with outer side length 0.8 and inner side length 0.4 and the other is a hollow circle with outer diameter 0.8 and inner diam-eter 0.4. The centers of the cylinders are initially placed at a loca-tion (0.8, 0.8). They are advected to the posiloca-tion (2.8, 1.8) after 1 unit of time. The domain is partitioned into either 100 100 rect-angular cells or 22478 trirect-angular cells. The trirect-angular mesh is illus-trated inFig. 6.
Theoretically, the geometry of the cylinders remains unchanged during advection due to the uniform velocity field. To quantify the errors produced by numerical methods, the solution error is de-fined by E ¼ P ijfin
D
V f e iD
Vj P ifiiD
V ð32Þ where fnis the numerical solution, fethe exact solution and fithe initial condition. The sum is over all the cells in the computational domain. The contours of the volume fraction over the range 0.05–0.95 in intervals of 0.1 are displayed inFig. 7for cell Courant numbers 0.75 and 0.1. The schemes adopted for comparison are the CICSAM of Ubbink and Issa[19], the HRIC of Muzaferija et al.[20], and the two schemes represented byEqs. (12) and (18), where the former is termed scheme A and the latter scheme B. The CICSAM tested here is the version without the correction step. The cell Courant number is defined byFig. 6. A mesh with 22478 triangular cells.
Rectangular mesh
Triangular mesh
Schemes
0.75
NC
=
0.1
C
N=
C
N=
0
.
7
5
CICSAM
HRIC
FBICS-A
FBICS-B
CN¼
PmaxðF
f;0Þ
D
tD
V ð33Þwhere Ffis the volumetric flux through the face and the sum is
ta-ken over all the faces of the cell. It is revealed from the figure that with the rectangular mesh, the results obtained by the CICSAM and HRIC are not acceptable for CN= 0.75. It was mentioned in the
last section that the CICSAM scheme is limited by a bound ~f
f ¼ ~fU=CN. It is reduced to the upwind scheme as CNapproaches
1. Thus, it is highly diffusive at Courant numbers close to 1. Even higher diffusive phenomenon can be found in the HRIC as a result of the use of upwind scheme for CNgreater than 0.7 in this scheme.
It needs to be pointed out that in implementation of the HRIC in the present calculations, it is the cell Courant number, instead of the face Courant number (=FfDt/DV, which is usually smaller than the
cell Courant number) as used in the original study, adopted in determining the flux blending. The contour plots for the schemes A and B of the present method are exactly the same because these two schemes are identical as the Courant number is greater than 1/ 3. The accuracy of the results for CN= 0.1 is much more improved.
The HRIC scheme blends the upwind difference and a bounded downwind difference for CN< 0.3. As a consequence of the effects
of the embedded upwind difference, smear of the interface is still visible. Both the predictions by the CICSAM and the scheme B are very satisfactory. Comparing with these two schemes, the contours obtained by the scheme A is only slightly more diffusive. Similar observations can be found for the case with circular cylinder. How-ever, the circle tends to transform into an octagon for CN= 0.1 due to
the large compression effects of the downwind difference. This is especially true for the CICSAM and HRIC schemes.
With use of the triangular mesh, the diffusion character of the CICSAM and HRIC at the high Courant number of 0.75 is not signif-icant. It should be noted that unlike the rectangular mesh, the Cou-rant numbers for the cells of the triangular mesh are not constant. Here the Courant number 0.75 simply represents the maximum value of its distribution. An analysis of the cell Courant number distribution indicates that nearly all the values fall in the range 0.4–0.5. It is also needs to be understood that the cell volumes of the triangular mesh are about half of those of the rectangular mesh. Keeping the same Courant numbers results in smaller time steps for the runs on the triangular mesh.
The effects of Courant number on the prediction accuracy for the considered schemes are shown inFig. 8a and b for the square cylinder and the circular cylinder, respectively. Due to the use of upwind difference in the flux blending in the HRIC, the numerical errors are much greater than the other schemes even at low Cou-rant numbers. The sharp increase of the errors at high CouCou-rant numbers for both the HRIC and CICSAM is owing to the approach to the upwind difference scheme. In general, the CICSAM is accu-rate for CN< 0.4. It can be identified that the present schemes A
and B perform satisfactorily, regardless of the Courant numbers. As expected, the scheme A is identical to the scheme B for CN> 0.33. For small Courant numbers the scheme B is slightly
infe-rior to the scheme A and the CICSAM. 5.2. Advection of a circle in a shear flow
In the above case, the cylinders are transported by the uni-form velocity field without changing its shape. In realistic prob-lems, the interface is subject to flow straining and deforms continuously. To mimic this situation, the following velocities are assumed.
~
V ¼ sin x cos y~i cos x sin y~j ð34Þ
with (x, y)
e
(0,p
). Although triangular meshes were also adopted in calculations, the results presented in the following are obtainedwith a 100 100 rectangular mesh. The initial condition is a circle of diameter 0.4
p
with its center located at (p
/2, (1 +p
)/5). The time marching of the computation is carried out over N units of time. It is followed by another N units of time in which the velocity field is re-versed. The circle will be stretched by the velocity straining in the forward step and recovers its original shape by the end of the back-ward step.Selected results for the CICSAM and the scheme A at CN= 0.25
and 0.75 for N = 8 and 16 are shown inFig. 9. For CN= 0.25 the
scheme A is slightly better than the CICSAM by comparing the final circles after the backward step. As the Courant number is increased to 0.75, the wide spread of the contours lines for the CICSAM indi-cates the smear caused by the numerical diffusion.
Fig. 10presents the solution error against the Courant number for the case with N = 16. As expected, the errors for the HRIC increases sharply for CN> 0.3 whereas the corresponding point for the CICSAM
occurs at CN= 0.5. The schemes A and B are not sensitive to the
Cou-rant number with the scheme A as the winner. The two schemes have less error throughout the entire range considered.
Courant number C Er ro r 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 1 1.2 1.4 CICSAM HRIC FBICS-A FBICS-B N
(a)
Courant number C Er ro r 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 1 1.2 1.4 CICSAM HRIC FBICS-A FBICS-B N(b)
Fig. 8. Numerical error against Courant number for (a) square cylinder and (b) circular cylinder in uniform flow obtained using rectangular mesh.
5.3. Collapse of a water column with or without an obstacle Experiments have been made to investigate the collapse of a water column[31,32], which has been widely used as a classic test case in the modeling of free surface flows. We consider a water col-umn with 0.146 m in width and 0.292 m in height, standing at the left corner in a tank with 0.584 0.340 m in size. The top side is an open boundary with the atmospheric pressure imposed. Since the velocity is variable with time and space, the Courant number for real flow calculations can not be fixed at a specific value. The time step is determined by setting the Courant numbers for all cells not greater than 0.25. Calculations were taken using the scheme A incorporating a 80 50 or 120 70 rectangular mesh and a 6218 or 9698 triangular mesh.
The variation of the position of the leading front and the reduc-tion of the height of the column are shown inFig. 11a and b. Both the distances and the time are nondimensionalized as shown in the figures. The predictions follow the measurements [31] closely. However, the calculated leading front moves faster than the exper-imental one. This may be owing to the difficulty to determine the exact position of the front in the experiments. Such a situation can also be found in the calculations of other studies[32,33].
The collapsing flow becomes much more complicated when an obstacle is placed in the way of the flow front. A rectangular block in the size 0.024 0.048 m is placed at the center of the floor. The height of the tank is increased to 0.584 m. The evolution of the re-sulted interface and fluid velocity, predicted with the Scheme A on
a 200 200 rectangular mesh, is shown inFig. 12for the times t = 0.1, 0.3, 0.4, and 0.6 s. In the early stage, the front of the collapsing water column proceeds along the floor. The movement of the front is then obstructed by the block. A tongue is formed due to the water 0.25 N C = CN=0.75 Schemes N=8 forward 8 N= backward 16 N= forward 16 N= backward 16 N= forward 16 N= backward CICSAM FBICS-A
Fig. 9. Results for the advection in shear flow obtained using rectangular mesh (contour range: 0.05–0.95 in intervals of 0.1).
Courant number C Err or 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 1 1.2 CICSAM HRIC FBICS-A FBICS-B N
Fig. 10. Numerical error against Courant number for the advection in shear flow obtained using rectangular mesh.
(a)
(b)
t*sqrt (2g/a) z/a 0 0.5 1 1.5 2 2.5 3 3.5 1 1.5 2 2.5 3 3.5 4 Exp. 1.125" Exp. 2.25" 80*50 rectangular mesh 120*70 rectangular mesh 6218 triangular mesh 9698 triangular mesh t*sqrt (g/a) h/(2a) 0 1 2 3 4 0.2 0.4 0.6 0.8 1 Experiment 80*50 rectangular mesh 120*70 rectangular mesh 6218 triangular mesh 9698 triangular meshFig. 11. Variation of (a) the leading edge and (b) the height of the collapsing water column without obstacles.
splash as seen at t = 0.3 s. It is followed by a strike of the tongue on the side wall. It can be detected that at t = 0.4 s, the water sheet formed on the wall has started to fall under the action of gravity. Air is trapped underneath the water tongue. This results in Ray-leigh–Taylor instability and, thus, a wavy form is seen on the lower surface, where a number of small ‘fingers’ can be identified. This kind of flow is very unstable. By the time t = 0.6 s, the structure of the water tongue is completely destroyed. The interface has broken and become extremely irregular. Photographs taken by Koshizuka at the corresponding times have been made available in reference [33]. Good qualitative agreement between the two studies can be seen, which justifies the current methodology to cope with compli-cated interfacial flows. However, small water droplets produced in the splash procedure are not found in the present predictions. To predict this behavior correctly, the mesh size needs to be refined to a level smaller than the droplets, which is prohibited in the cur-rent stage. The plots of velocity vectors indicate that the splash has a great effect on the air flow. A number of vortex flows are generated during evolution of the free surface flow.
6. Conclusions
A solution algorithm has been developed to deal with interfacial flows in two-fluid flow systems. It is based on the finite volume method and is applicable to unstructured meshes of arbitrary topology. It features the use of flux limiters to blend high resolu-tion schemes and compressive schemes to determine the convec-tive fluxes through cell faces. The blending factor is a function of the angle between the interface and the mesh lines. Two such schemes have been introduced. The limiters of the first scheme are independent of Courant number and the second scheme re-duces to the first one as the Courant number is greater than 1/3. Rigorous tests on these schemes reveal that they outperform the well-known CICSAM and HRIC schemes, especially for large Cou-rant numbers. To avoid the diffusion effects encountered by the last two schemes, the cell Courant numbers are suggested to be less than 0.4 for them. On the other hand, the present schemes maintain their accuracy throughout the range 0 < CN< 1.
Acknowledgment
The support of National Science Council of the Republic of China is acknowledged.
References
[1] B. Ramaswamy, M. Kawahara, Lagrangian finite element analysis applied to viscous free surface fluid flow, Int. J. Numer. Methods Fluids 7 (1987) 953– 984.
[2] J. Fukai, Z. Zhao, D. Poulikakos, C.M. Megaridis, O. Miyatake, Modeling of the deformation of a liquid droplet impinging upon a flat surface, Phys Fluids A: Fluid Dyn. 5 (1993) 2588–2599.
[3] S. Muzaferija, M. Peric, Computation of free-surface flows using the finite-volume method and moving grids, Numer. Heat Trans. B: Fund. 32 (1997) 369– 384.
[4] S.O. Unverdi, G. Tryggvason, A front-tracking method for viscous, incompressible, multi-fluid flows, J. Comput. Phys. 100 (1992) 25–37. [5] G. Tryggvasson, B. Bunner, A. Esmaeeli, D. Juric, N. Al-Rawahi, W. Tauber, J.
Han, S. Nas, Y.-J. Jan, A front-tracking method for the computations of multiphase flow, J. Comput. Phys. 169 (2001) 708–759.
[6] S. Osher, J.A. Sethian, Fronts propagating with curvature-dependent speed: algorithms based on Hamilton–Jacobi formulations, J. Comput. Phys. 79 (1988) 12–49.
[7] M. Sussman, P. Smereka, S. Osher, A level set approach for computing solutions to incompressible two-phase flow, J. Comput. Phys. 114 (1994) 146–159. [8] D. Lakehal, M. Meier, M. Fulgosi, Interface tracking towards the direct
simulation of heat and mass transfer in multiphase flows, Int. J. Heat Fluid Flow 23 (2002) 242–257.
[9] F.H. Harlow, J.E. Welch, Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface, Phys. Fluids 8 (1965) 2182– 2189.
[10] J.J. Monaghan, Simulating free surface flows with SPH, J. Comput. Phys. 110 (1994) 399–406.
[11] S. Koshizuka, A. Nobe, Y. Oka, Numerical analysis of breaking waves using the moving particle semi-implicit method, Int. J. Numer. Methods Fluids 26 (1998) 751–769.
[12] H.Y. Yoon, S. Koshizuka, Y. Oka, A particle-gridless hybrid method for incompressible flows, Int. J. Numer. Methods Fluids 30 (1999) 407–424. [13] W.F. Noh, P. Woodward, SLIC (simple line interface calculation), Lecture Notes
Phys. 59 (1976) 330–340.
[14] D.L. Youngs, Time-dependent multi-material flow with large fluid distortion, in: K.W. Morton, M.J. Baines (Eds.), Numerical Methods for Fluid Dynamics, Academic Press, New York, 1982, pp. 273–285.
[15] W.J. Rider, D.B. Kothe, Reconstructing volume tracking, J. Comput. Phys. 141 (1998) 112–152.
[16] M.R. Davidson, M. Rudman, Volume-of-fluid calculation of heat or mass transfer across deforming interfaces in two-fluid flow, Numer. Heat Transf. B: Fund. 41 (2002) 291–308.
[17] M. Rudman, Volume-tracking methods for interfacial flow calculations, Int. J. Numer. Methods Fluids 24 (1997) 671–691.
[18] M. Rudman, A volume-tracking method for incompressible multifluid flows with large density variations, Int. J. Numer. Methods Fluids 28 (1998) 357– 378.
[19] O. Ubbink, R.I. Issa, A method for capturing sharp fluid interfaces on arbitrary meshes, J. Comput. Phys. 153 (1999) 26–50.
[20] S. Muzaferija, M. Peric, P. Sames, T. Schellin, A two-fluid Navier–Stokes solver to simulate water entry, in: Proceeding of Twenty-Second Symposium On Naval Hydrodynamics, Washington, DC, 1998, pp.638–649.
[21] M. Darwish, F. Moukalled, Convective schemes for capturing interfaces of free-surface flows on unstructured grids, Numer. Heat Transf. B: Fund. 49 (2006) 19–42.
[22] Y.-Y. Tsui, Y.F. Pan, A pressure-correction method for incompressible flows using unstructured meshes, Numer. Heat Transf. B: Fund. 49 (2006) 43– 65.
[23] Y.-Y. Tsui, T.-C. Wu, A pressure-based unstructured-grid algorithm using high-resolution schemes for all-speed flows, Numer. Heat Transf. B: Fund. 53 (2008) 75–96.
[24] Y.-Y. Tsui, T.-C. Wu, Use of characteristic-based flux limiters in a pressure-based unstructured-grid algorithm incorporating high-resolution schemes, Numer. Heat Transf. B: Fund. 55 (2009) 14–34.
[25] P.K. Sweby, High resolution schemes using flux limiters for hyperbolic conservation laws, SIAM J. Numer. Anal. 21 (1984) 995–1011.
[26] B.P. Leonard, The ULTIMATE conservative difference scheme applied to unsteady one-dimensional advection, Comp. Methods Appl. Mech. Eng. 88 (1991) 17–74.
[27] P.H. Gaskell, A.K.C. Lau, Curvature-compensated convective transport: SMART, a new boundedness-preserving transport algorithm, Int. J. Numer. Methods Fluids 8 (1988) 617–641.
[28] R.I. Issa, Solution of the implicitly discretised fluid flow equations by operator-splitting, J. Comput. Phys. 62 (1986) 40–65.
[29] Y.Y. Tsui, S.P. Jung, Analysis of the flow in grooved pumps with specified pressure boundary conditions, Vacuum 81 (2006) 401–410.
Fig. 12. Evolution of the collapsing flow of the water column with an obstacle (contour range: 0.4–0.6).
[30] I. Orlanski, A simple boundary condition for unbounded hyperbolic flows, J. Comput. Phys. 21 (1976) 251–269.
[31] J.C. Martin, W.J. Moyce, An experimental study of the collapse of liquid columns on a rigid horizontal plane, Philos. Trans. Roy. Soc. Lond. A: Math. Phys. 244 (1952) 312–324.
[32] S. Koshizuka, H. Tamako, Y. Oka, A particle method for incompressible viscous flow with fluid fragmentation, Comput. Fluid Dyn. J. 4 (1995) 29–46. [33] O. Ubbink, Numerical prediction of two fluid systems with sharp interfaces,
![diagram in the CICSAM scheme of Ubbink and Issa [19] , where C N is](https://thumb-ap.123doks.com/thumbv2/9libinfo/7630351.134875/5.892.61.439.290.551/diagram-cicsam-scheme-ubbink-issa-c-n.webp)
