VOF method

Advances in Water Resources 34 (2011) 1320–1334

Contents lists available at ScienceDirect

Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres

Numerical simulation of overflow at vertical weirs using a hybrid level set/VOF method Xin Lv a, Qingping Zou b,⇑, Dominic Reeve a a b

School of Marine Science and Engineering, University of Plymouth, PL4 8AA, UK Department of Civil and Environmental Engineering, University of Maine, Orono, Maine 04469, USA

a r t i c l e

i n f o

Article history: Received 17 October 2010 Received in revised form 30 April 2011 Accepted 18 June 2011 Available online 26 June 2011 Keywords: Overflow Broad-crested weir Sharp-crested weir Free surface Unstructured grid VOF

a b s t r a c t This paper presents the applications of a newly developed free surface flow model to the practical, while challenging overflow problems for weirs. Since the model takes advantage of the strengths of both the level set and volume of fluid methods and solves the Navier–Stokes equations on an unstructured mesh, it is capable of resolving the time evolution of very complex vortical motions, air entrainment and pressure variations due to violent deformations following overflow of the weir crest. In the present study, two different types of vertical weir, namely broad-crested and sharp-crested, are considered for validation purposes. The calculated overflow parameters such as pressure head distributions, velocity distributions, and water surface profiles are compared against experimental data as well as numerical results available in literature. A very good quantitative agreement has been obtained. The numerical model, thus, offers a good alternative to traditional experimental methods in the study of weir problems. Ó 2011 Elsevier Ltd. All rights reserved.

1. Introduction Weirs are a small overflow-type dams commonly used to raise the level of a river or stream and cause a large change of water level behind them. The overflow at a weir is used in many engineering scenarios for safety reasons, such as to discharge excess water from a river in times of flooding. The flow over the weir, in the mean time forms a jet which in turn hits the lower level water downstream of weir structures. The strong impact of the water jet after the weir usually causes the entrainment of air bubbles, which is also often used in water treatment processes to alter the concentration of dissolved gases, to strip volatile organics, or to reduce tastes and odours. An experimental study assessing the effects of different shapes of weirs on the air entrainment rate can be found in Emiroglu and Baylar [1]. It is obvious that the accurate prediction of the flow over weirs is of practical interest to hydrologists and engineers. Generally, this can be done either by deriving empirical formulae based on extensive experimental analysis, or by relatively cheaper numerical simulations. The experimental study of the hydraulic characteristics of different types of weirs can be dated back to the 19th and 20th centuries. To date, numerous studies have already been carried out in laboratories in order to reproduce overflows, measuring discharges and the free surface profile, such as for example Rouse and

⇑ Corresponding author. E-mail address: [email protected] (Q. Zou). 0309-1708/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.advwatres.2011.06.009

Reid [2], Bagheri and Heidarpour [3] and Rajaratnam and Muralidhar [21]. Compared to the experimental methods, a properly validated numerical model has many attractive merits, such as permitting one to obtain the flow characteristics of a weir for a wide range of hydraulic parameters without recourse to expensive and much more time consuming experiments in the laboratory. Furthermore, the numerical model also permits one to incorporate small changes in the geometric parameters, including changes in inlet and outlet conditions and study their impact on the weir flow characteristics. However, only very few numerical simulations of this topic have been published (e.g. [4,5]). This is due to the high complexity of this particular type of flow, requiring rather complex physical models and their accurate and robust numerical discretization. The classical shallow water equations (SWE) based models that are often used in hydraulic applications are not appropriate to simulate the flow over a weir. SWE models assume hydrostatic pressure distribution in the flow and neglect vertical velocity to justify depth-averaging approach. These assumptions are not valid in this problem since the flow has a strong variation in the vertical and since the flow over the weir exhibits a jet-type structure. A more sophisticated numerical model, such as a Navier–Stokes based solver and a sophisticated surface capturing scheme are thus required to reproduce accurately the physics of free surface flow downstream of the weir. Simulating free surface flow, however, remains a formidable challenge because non-linear boundary conditions are required to be satisfied on an arbitrarily moving surface whose position is

X. Lv et al. / Advances in Water Resources 34 (2011) 1320–1334

not prescribed a priori. Two main approaches have been adopted for tracking the free surface: the Lagrangian and Eulerian. Both approaches have advantages and weaknesses, and the focus of this paper will be on the latter option. Tracking the free surface using the Eulerian description employs an equation for the density:

@q þ u  rq ¼ 0; @t

ð1Þ

which is referred to as the transport equation. The density qðx; tÞ can be either continuous or discontinuous. By tracking the change of q it is possible to identify the location of the free surface. Direct numerical solution of above equation leads to excessive numerical diffusion. As a result, the free surface will become thicker and thicker and eventually become diffused into the whole domain. Hence other approaches have been used for such a purpose. Nichols et al. [23] developed the volume of fluid (VOF) method to track the free surface. By using the transport equation and assuming that density is everywhere constant ðq ¼ q0 Þ in the flow domain and zero ðq ¼ 0Þ in the air domain, it is possible to normalize the transport equation by q0 . Defining F ¼ q=q0 as the fractional volume of fluid function (VOF function), the transport equation becomes:

@F þ u  rF ¼ 0: @t

ð2Þ

F is a step function, where F ¼ 1 in cells containing fluid and F ¼ 0 in cells containing air. In free surface cells, the value of F ranges between (0–1). Therefore, by tracking the VOF function F, one can identify the free surface elements at any time step. VOF methods have proved to be robust and relatively easy to code, but still have drawbacks [24–26]. The overall accuracy of this method relies heavily on the performance of its interface reconstruction scheme. Another technique for interface tracking is the level set method [27]. By substituting q with the smooth level set function u Eq. (1) becomes

@u þ u  ru ¼ 0: @t

ð3Þ

The function u gives the normal distance from the interface, in which its original value is set to zero at the interface. By tracking the zero value of u based on Eq. (3), one may identify the interface during the computation. Because the set level function u is a smooth function, it is much easier to compute and it induces much less numerical diffusion. But it has been shown [28] that this method is subjected to significant mass loss under complicated situations since the scheme does not explicitly impose mass conservation. The ghost fluid method (GFM) was developed based on level set method in its original form [29]. Essentially, GFM exploits the concept of ghost and real fluid cells and manages them with an overlapping Schwarz like numerical procedure. In the material interface region, it sets the values of the pressure and normal velocity in the ghost fluid cells to those in the real fluid cells. And to eliminate an otherwise spurious ‘‘over-heating’’ phenomenon, this method computes the density of the ghost fluid using an isobaric fix technique. As explained in [29], this isobaric fix requires the solution of yet another auxiliary partial differential equation and therefore increases further the computational complexity of the method. Also as pointed out by some researchers (for example [30]), the original GFM fails to solve some air/water problems of interest due to the large density ratio. In [31], Kang and the co-authors extended GFM to multiphase incompressible flow including the effects of viscosity, surface tension and gravity. Lv et al. [6] developed a new hybrid Level set/Volume of Fluid, (LS/VOF), method based on Navier–Stokes equations for free surface flows using three dimensional unstructured tetrahedral grids. At each time step, the level set function is evolved by solving the

1321

level set advection equation using a high resolution characteristic based finite volume method and the volume fraction advection is performed using a bounded compressive Normalised Variable Diagram (NVD) based scheme (Leonard, [8]). The interface is reconstructed based on both the level set and the volume fraction information. In particular, the interface normal vector is calculated from the level set function while the intercepts are determined by enforcing mass conservation based on the volume fraction. The novelty of the method lies in the fact that an analytic method is used to define the intercepts on tetrahedral grids, which makes interface reconstruction efficient and conserves volume of fluid exactly. The level set function is then reinitialized to the signed distance to the reconstructed interface. Furthermore, the adaptive combination of high resolution discretization schemes ensures the preservation of the sharpness and shape of the interface while retaining boundedness of the volume fraction field. Since the level set function is continuous, the interface normal vector calculation is straightforward and accurate compared to a classic volume-offluid method, while tracking the volume fraction is essential for enforcing mass conservation. The method is also coupled to a well validated finite volume based Navier–Stokes incompressible flow solver [9]. In their latest work [32], Lv et al. improved the algorithm by introducing a fully implicit method so as to further reduce the total CPU time, as well as take the air phase effects into account. In this paper the overflow over broad- and sharp-crested weirs is investigated using the numerical technique described in Lv et al. [6]. The vortical motions, air entrainment and pressure variations due to violent collision after the overtopping will be examined as well as the evolution of the vertical distributions of pressure, velocity and turbulence downstream the weir crest. The outline of the paper is as follows: the governing partial differential equations (PDE) are described in Section 2. The numerical methods used for PDE discretization and free surface capture are introduced in Section 3. The results of our numerical simulations and the detailed comparison with experimental as well as numerical data available in literature are presented and discussed in Section 4. Concluding remarks are given in Section 5.

2. Governing equations Here we consider incompressible flows with two immiscible fluids. The density of one fluid is qa and the density of the second fluid is qg . The non-dimensional governing 3D equations (including incompressible Navier–Stokes equations and level set equation), modified by the Artificial Compressibility Method (ACM) [11], are given as

1 @p ~ ¼ 0; þrU b @s

ð4Þ

au u @p @u @u ~ ¼  1 @p þ 1  r2  u þ 1 F Sx þ F gx ; þ r  ðuUÞ þ þ qb @ s @ s @t q @x qRe q ð5Þ au v @p @ v @ v ~ ¼  1 @p þ 1  r2  v þ 1 F Sy þ F gy ; þ þ þ r  ðv UÞ qb @ s @ s @t q @y qRe q ð6Þ   au w @p @w @w ~ ¼  1 @p þ 1  r2  w þ 1 F Sz þ F gz ; þ r  wU þ þ qb @ s @ s @t q @z qRe q ð7Þ @u @u ~ ¼ 0; þ þ r  ðuUÞ @s @t

ð8Þ

1322

X. Lv et al. / Advances in Water Resources 34 (2011) 1320–1334

where s is the pseudo time introduced to perform dual-time stepping [9,10] and t the physical time. The velocity vector is ~ ¼ u ~ j þ w ~ k. Eqs. (4)–(8) can be expressed in non-dimenU i þ v ~ sional vector form as follows:

C

@W @W þ r ~ þK Fc ¼ r  ~ Fv þ ~ S; @s @t

ð9Þ

where

2

p

3

6 7 6u7 6 7 7 W¼6 6 v 7; 6 7 4w5

u

2

0 6 60 6 K¼6 60 6 40 0 2

3 ~ U 7 6 ~ 6 uU þ p=qdij 7 7 6 7 6 ~ ~ Fc ¼ 6 v U þ p=qdij 7; 7 6 7 6 wU 4 ~ þ p=qdij 5 2

~ uU

0

0

0

0

7 07 7 0 1 0 07 7; 7 0 0 1 05 0 0 0 1 3 0 7 61 6 q F Sx þ F gx 7 7 6 7 6 ~ S ¼ 6 q1 F Sy þ F gy 7: 7 6 61F þF 7 gz 5 4 q Sz 0 1 0

2

3

0

C¼

1 b 6 au u 6 qb 6 6 au v 6 qb 6 6 au w 4 qb

0

3 0 7 6 1 6 qRe  r  u 7 7 6 7 6 ~ Fv ¼ 6 q1Re  r  v 7; 7 6 6 1  r  w7 5 4 qRe 0 2

0 0

0

1 0

0

0 0

0

0

q ¼ He ðuÞqa þ ð1  He ðuÞÞqg ; l ¼ He ðuÞla þ ð1  He ðuÞÞlg :

The total density conservation equation is replaced by the level set equation, and they are equivalent if qa and qg are constant. What is presented above constitutes the standard level set method, where the advection of u, including a re-initialization step to retain u as a signed distance function, is not done in a conservative way, not even for divergence free velocity fields. This implies that the total mass bounded by the zero level set is not conserved. This drawback is mitigated by introducing the volume fraction (VOF) equation Eq. (2), which is repeated as follows after some simple mathematical derivations:

@F ~ ¼ 0: þ rðF  UÞ @t

ð13Þ

VOF methods have proved to be robust and accurate in mass conservation, but the overall accuracy of this method relies heavily on the performance of its interface reconstruction scheme. The problem of how to advect the interface with minimum numerical diffusion and dispersion was discussed in detail by Lv et al. [6], and the problem of coupling the level set function u and the volume fraction F will be discussed later in this paper. The overall target is to maintain the boundedness of the F function to avoid any nonphysical deformation of the interface shape.

3

7 07 7 7 0 1 0 0 7; 7 0 0 1 07 5 1

3. Numerical schemes 3.1. Navier–Stokes and level set equations discretization

In all the equations above, W is the vector of dependent variables, p and q are pressure and density, respectively, b the constant parameter introduced by ACM (typically set to 10). ~ Fc and ~ Fv are the ~ advective and viscous flux vectors. S contains the surface tension and gravity terms. The first term on the left-hand side of Eqs. (5)– (7) is the partial derivative of pressure with respect to pseudo-time s (the artificial compression term), which is introduced to couple velocity and pressure fields for the calculation of pressure based on the divergence-free condition. C is a preconditioning matrix with an arbitrary constant au (=1 as proposed by Turkel [12]). K is the unit matrix with its first element being zero. Since the surface tension and gravity play significant roles during the development of a free surface flow process, the combined effects for surface tension and gravity are included. F g is the gravity force per unit mass and ~ n is given as ~ F g ¼ Frg2 ; where Fr is the Froude number and ~ ng the unit vector along the prescribed direction of gravity. ~ F S is the surface tension force per unit volume, given by Brackbill et al. [13]. The level set function u is defined to be a signed distance function

juð~ xÞj ¼ dð~ xÞ ¼ minðj~ x ~ xI jÞ;

Following the work by Zhao and Tai [6,9,10], Eqs. (4)–(8) are firstly recast in integral form as

C

@ @s þ

ð11Þ

where e corresponds to half the thickness of the interface. It should not be too small in order to maintain robustness of the numerical scheme. In our proposed method, 1  2 times of the smallest grid size would be enough to keep the code stable and accurate. Finally, the constitutive relations for the density and dynamic viscosity are defined to ensure closure of the system;

Z Z Z

WdV þ K

V

@ @t

Z Z Z

WdV þ

V

I Sc v

ð~ Fc  ~ Fv Þ  ~ ndS

~ SdV ¼ 0;

ð14Þ

where ~ n is the surface unit normal vector pointing out of the control volume. Eq. (14) is then discretized on an unstructured tetrahedral grid and a cell-vertex scheme is adopted here, i.e., all computed variables in vector W are stored at the vertices of tetrahedral cells.

P

b 3

a c

where I is the interface, u > 0 on one side of the interface and u < 0 on the other. In practical computations, to achieve numerical robustness, a smeared out version of the Heaviside function is used,

e 6 u 6 e;

Z Z Z

V

ð10Þ

xI 2I

8 > < 0; u < e;   He ðuÞ ¼ 12 þ 2ue þ 21p sin peu ; > : 1; u > e;

ð12Þ

2

1

O

A

C Fig. 1. Construction of control volume within a tetrahedron for a node P.

B

1323

X. Lv et al. / Advances in Water Resources 34 (2011) 1320–1334

For every vertex, as shown in Fig. 1, a control volume is constructed using the median dual of the tetrahedral grid: nodes A, P, B and C form the vertex of the tetrahedral cell and O is the centre of the element APBC; a, b and c are the median duals of the edges AP, BP and CP, respectively; and 1, 2 and 3 are the centre points of triangles APC, CBP and ABP, respectively. In the cell-vertex scheme, the computed variables are stored at vertices A, P, B and C. Triangles O1a and O3a form the control volume surface for the convective term of the governing equation on edge AP. Likewise, the rest of the convective term for different edges is computed in an analogous manner. Triangles APC, CBP, ABC and ABP form the corresponding control volume surfaces for the calculation of viscous terms. Finite volume numerical discretizations are based on the solution of the governing equations in integral form and spatial discretization is performed by using Eq. (14). And the partial implicit matrix-free method proposed in Lv et al. [6] is employed to March the solution in physical time. 3.2. VOF equation discretization and free surface reconstruction The evolution of the volume fraction field F is governed by Eq. (13). The upwind-biased characteristic method developed by Zhao et al. [10] is too diffusive for the evolution of step function. In this work, the CICSAM (Compressive Interface Capturing Scheme for Arbitrary Meshes) discussed in Lv et al. [6] is employed for its good performance to maintain the sharpness of the interface while keeping reasonable accuracy. It makes use of the NVD concept [8] and switches between different high resolution differencing schemes to yield a bounded scalar field, but one which preserves both the smoothness of the interface and its sharp definition (over one or two computational cells). The derivation of the scheme is based on the recognition that no diffusion of the interface (whether physical or numerical) can occur; thus it is only appropriate for sharp fluid interfaces. The scheme is based on the finite volume technique and is fully conservative. The special implicit implementation embedded inside makes it applicable to unstructured meshes. Interface reconstruction is the essential part of the coupled VOF method, i.e., given the interface normal vector in a cell, the unique linear interface segment (line segment in 2D and planar surface in 3D) which also truncates the cell by the given volume fraction must be found. Analytical relations connecting linear interfaces and volume fractions in triangular and tetrahedral grids have been presented in Yang et al. [14]. The merit of this method is that the analytic formulations eliminate the need for iteration, and thereby improve accuracy and reduce computation time. (The volume of fluid is exactly conserved during this form of interface reconstruction). After the planar interface in each free surface tetrahedron (defined by F ¼ 0:5) has been reconstructed analytically, the level set field needs to be re-distanced, which is also an essential step in normal level set method. This can be easily done for each grid point in the computational domain by finding the shortest distance to the reconstructed VOF interface, and an efficient lookup algorithm within the group of VOF planar interfaces is needed to save CPU time. For details on the implementation of the finite volume discretization of VOF equation and free surface reconstruction, the reader is referred to Lv et al. [6]. 3.3. Coupling between VOF and level set field The level set function u and the volume fractions F are coupled as follows: 1. The normal vectors (or slopes) used in the volume-of-fluid reconstruction step are determined from the level set function;

2. The volume fractions are used, together with the slopes from the level set function, to construct the ‘volume preserving’ zero level set interface as well as a distance function along with providing ‘closest point’ information to the zero level set; 3. The smooth level set function is used to express the interfacial curvature; 4. The ‘‘re-distanced’’ level set function is then used to clip the VOF function within the narrow band near the zero level set. 3.4. Air phase deactivation and bubble treatment It is well known that solving a two phase flow can cause severe problems when the density ratio between the liquid and gaseous phase is large. This difficulty mainly comes from the fact that in such a case, the pressure gradient distribution is discontinuous across the interface. This implies that any small pressure gradient that is erroneously transferred from the liquid to the air region due to numerical error will accelerate the air considerably. This in turn will lead to loss of divergence, causing more spurious pressures. The whole cycle may, in fact, lead to a complete divergence of the solution. It should also be pointed out that, some semi-implicit schemes have the merit handling very large density ratios without losing accuracy and stability. Among others, the projection method has received much more attention recently, and it has been widely used by many researchers in practical applications of free surface flow computations [15,16]. The advantage of this type of multiphase model lies in the fact that the interaction between gaseous phase and liquid phase can be accurately predicted without resorting to approximation or empirical formulae. But this comes with the cost of more computational resources and CPU time in that all of the gaseous phase grid points must be involved in the computation even they are far away the interface and of little interest. Also, for large scale three dimensional computations, the parallelization of projection method is not easy to implement. Following the work done in Lv et al. [6], a switch subroutine has been designed and its function is to monitor the density ratio between the two phases. If it finds the ratio is too large, say above 200, it automatically deactivates the computation for the air phase. All of the necessary information needed in the computation of liquid phase is extrapolated from inside the liquid phase itself. The proposed method is single phased in its nature. We only solve the water phase and the effects of surrounding air are usually ignored. On the free surface velocity extrapolation is performed in such a way that the tangential stress boundary condition is fulfilled. For entrapped air bubbles we assume that at each instant the pressure in the bubble is (spatially) constant. As long as the bubble is not in contact with the atmospheric air (see Fig. 2 for instance), the pressure can be obtained from the isentropic relation:

Air Phase

Liquid Phase

Free Surface Entrapped Bubble

Fig. 2. Sketch of bubble in the water.

1324

Vb V b0

c ð15Þ

;

where pb ; V b denote the pressure and volume in the bubble and pb0 ; V b0 the reference values (e.g. those at the beginning of the simulation). And c is the ideal gas constant (typically set to 1.4). The gas in the bubble is marked by solving a scalar advection equation of the form given as follows [17]:

@b þ v a  rb ¼ 0; @t

ð16Þ

where b is a scalar defined to be 1 in the bubble and 0 elsewhere. At the beginning of every time step the total volumes occupied by gas is added. From this volume, the pressure is computed from Eq. (15). At the end of every time step, a check is performed to see if the bubble has reached contact with the air. Should this be the case, the pressure in the open bubble is set to the atmospheric pressure. According to our experiences, this treatment can preserve the total mass in the bubble accurately. 3.5. Turbulence modelling The Large-Eddy simulation (LES) has been chosen to calculate turbulent phenomena. In this work, a dynamic form of Smagorinsky SGS model has been implemented to calculate the SGS stress tensor [7]. This model relies upon the Germano identity, which has been generalized in order to be applied to other subgrid-terms arising in the filtered energy equation. Furthermore, an improved formulation of the dynamic mixed model has been proposed for better representing the backscatter of turbulence energy which has been proven to be important for predicting flow characteristics. The unstructured grid filtering was adapted from a new filtering approach based on the least-squares technique. This approach can filter a function to any given level of commutation error on unstructured grids. 4. Model verification In what follows, the proposed model will be applied to simulate the flow over a weir. Two different weir shapes are considered in this work, namely broad and sharp-crested weirs. The computations are validated against empirical results, as well as experimental data. 4.1. Overflow at a rectangular broad crested weir In this section the numerical model described above is applied to reproduce a classic weir problem. In Fig. 3 a schematic view in the x–y plane of the computational domain is plotted. As can be seen from the figure, the moving material interface between air and water freely evolves in space and in time. The weir contours

and the bottom consist of solid constraints on the fluid motion. The flow is forced by gravity and viscous effects can be neglected. From the theory of potential flow for an inviscid fluid, one can easily conclude that the pressure distribution in the nappe (the part of the water overflows the weir) is not hydrostatic. Furthermore, the nappe has a contraction and moves like a jet. Simulating correctly such overflow requires a numerical model able to handle nonhydrostatic pressure fields and rapidly varying gravity-driven free surface flows. Thus, it becomes obvious that the shallow water equations (SWE) are ill-suited to reproduce the physics of such weir flow. Instead, the Navier–Stokes equations, described in Section 2, seem to be more appropriate for this particular application since they are directly derived from conservation principles, without any of the simplifying SWE assumptions. The domain consists of a 40 m long reservoir, divided at x = 25 m into two parts by a vertical wall of 5.5 m height. The wall consists of a fixed broad crested weir of height Rd ¼ 4 m (see Fig. 1) and a removable gate on the top. On the left of the reservoir, water (with an initial depth that ranges from 0.1  0.8 m in different testing runs) is at rest. On the right, the tank is empty. At the initial conditions, the gate is instantaneously removed and the water head over the crest of the weir, i.e. Rc in Fig. 3, ranges from 0:1  0:8 m. The domain is discretized by 304,497 nodes and 1,443,577 tetrahedral elements with a local refinement in the region surrounding the vertical wall. This corresponds to the grid size Dh ¼ 0:01 m within the refined region and Dh ¼ 0:04 m elsewhere. Due to the 2D nature of this testing case, only one cell length was utilized in the transverse direction. A fixed time step size Dt ¼ 0:0015 s was used in the tests and the total simulation time for each run was 35 s. The simulation was run on a workstation with 8 Intel Xeon CPUs running at 2.66 GHz, and the total CPUs time was 7,651 s (2.125 h). Figs. 4 and 5 plot the predicted results of the unsteady overflow evolution when Rc ¼ 0:8 m. The colours (Fig. 4) denote the pressure field. The inflow boundary condition with prescribed velocity was imposed on the left edge to keep the water level in the left tank constant. Instantaneously after removing the top wall, water flows over the weir driven by gravity and inertial forces and moves freely in air. After around 1.5 s the flow stream impacts on the right bottom of the tank with large deformations. The nappe downstream of the weir has a contraction and looks like a jet. The local pressure field increases considerably due to the violent collision, as shown by the colours in Fig. 4 (T = 0.35 s). Although it keeps the hydrostatic distribution in the tank on the left, the unsteady overflow produces a perturbation in the pressure field that progressively moves upstream through the nappe according to our simulation, as well as the findings by other researchers [4]. As the simulation reveals, the impact produces a bifurcation in the flow and part of the overflowing discharge moves backward against the weir, and in a recirculation flow will finally be formed. Water adheres to the bottom and on the weir yielding an

Atmosphere

Removed top wall 6

4

Inflow

RC Rd=4

U=0.01 m/s

Y

pb ¼ pb0

X. Lv et al. / Advances in Water Resources 34 (2011) 1320–1334



2

BL=1

0 10

20

30

40

Z

X

Non-Slip wall Fig. 3. Sketch of the computational domain for overflow at broad crested weir. The basic boundary condition setup is also shown in the figure.

X. Lv et al. / Advances in Water Resources 34 (2011) 1320–1334

1325

Fig. 4. Snapshots of the free surface profile at the initial stages of the overflow. The pressure contour field is plotted. The numerical results refer to output time t = 0.15 s, 0.35 s, 0.55 s, 1.05 and 1.25 s respectively. The pressure field was non-dimensionalized using the reference pressure pref ¼ qwater U 2ref .

overturning breaking wave. This is well represented by Fig. 5, where strong vortices due to recirculation may be easily identified. It is interesting to notice the formation of a big bubble due to the

interactions between the forward moving overflowing discharge and the wave bounced back by the right wall of the tank. During its movement towards the right, the bubble keeps changing its

1326

X. Lv et al. / Advances in Water Resources 34 (2011) 1320–1334

Fig. 5. Snapshots of the free surface profile at the later stages of the overflow. The streamlines and velocity contour are shown. The numerical results refer to output time t = 13.125, 13.5, 13.875 and 14.25 s. The velocity magnitude was normalised by the maximum value in the flow field.

shape under the combined effects of gravity, friction and surface tension. It should be pointed out that due to the nature of the single-phased treatment within the model, the velocity field in the air neighbouring to the free surface was constructed by extrapolating from within the water phase. This was reflected by the continuous streamlines in Fig. 5.

In their studies, Chadwick and Morfett [18] expressed the discharge formula over a broad crested weir as follows:

qweir ¼ 1:705  C d jRc j3=2 ;

ð17Þ

where Rc is the overflow depth (freeboard) and C d the discharge coefficient.

1327

X. Lv et al. / Advances in Water Resources 34 (2011) 1320–1334

C d ¼ 0:848  C F ;

  Rc Rc  0:35 ; C F ’ 0:91 þ 0:21  þ 0:24 BL Rc þ Rd

ð18Þ ð19Þ

where, BL and Rd the weir width and height respectively. The comparison between discharge rate of weir equation (Eq. (17)) and the current prediction is presented in Table 1 and Fig. 6. It can be seen from Fig. 6 that the results produced from the current numerical model compare well with the weir equation for all of the negative freeboard. The model under predicts the overflow discharge as much as 10% compared to the weir equation. This may be related to the uncertainty in discharge coefficient ðC d Þ value. Another reason for the difference between the weir equation and current numerical discharge rates is that the weir equation is based on the Bernoulli equation as a starting point. The actual flow over a weir is complex, usually being unsteady and involving hydrodynamic pressure distribution. These effects are not covered in the Bernoulli equation. To determine the accuracy of the proposed methods, we carried out a grid convergence study for the overflow problem. Four successively finer mesh sizes corresponding to grid point spacing Dh ¼ 0:03 m; Dh ¼ 0:02 m; Dh ¼ 0:01 m, and Dh ¼ 0:005 m within the refined region near the vertical wall, respectively, were used for error analysis, and the weir equation (Eq. (17)) solution was taken as the baseline solution. On all the grids the same physical time step ðDt ¼ 0:0015 sÞ was employed in order to concentrate on the spatial resolution of the method. The results of the grid convergence study are summarized in Table 2. From the data provided in Table 2, it may be seen that the grid with Dh ¼ 0:01 m produces a reasonable result compared to the rest coarse grids if taking the result on the finest grid ðDh ¼ 0:005 mÞ as the baseline solution. Considering the increase of total CPU time (which is 107.7% according to our test) on the finest grid against the next finest grid ðDh ¼ 0:01 mÞ, the grid with Dh ¼ 0:01 m is considered to produce a good balance in terms of accuracy and efficiency. It should be noted that the error measurement is based on the overflow discharge, which is a general data measuring a volume effect rather than providing the answer whether the solutions are point-wise convergent. However, in the case of lacking analytical solutions, as well as experimental flow field measurements, the discharge, qweir , is considered a good option in evaluating the overall performance of the numerical model. The readers are referred to [6] for the comprehensive convergence study of our model. 4.2. Overflow at a sharp crested weir In this section, the model is applied to reproduce the overflow at a sharp-crested weir, which serves as a convenient and accurate

Table 1 Comparison between numerical prediction and weir equation Eq. (17). qweir was calculated with Eq. (17). Run No.

Rc ðmÞ

qweir ðm =m=sÞ

qpredicted ðm =m=sÞ

ðqpredicted  qweir Þ=qweir ð%Þ

1 2 3 4 5 6 7 8 9

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

0.000 0.039 0.114 0.215 0.341 0.490 0.662 0.857 1.075

0.000 0.035 0.111 0.220 0.339 0.484 0.637 0.812 1.001

0.0 10.2 2.6 2.3 0.6 1.2 3.8 5.3 6.9

3

3

1.5 Overflow discharges (m^3/m/s)

A number of empirical discharge formulae have been developed which incorporate C d . In [18] the equations for C d are in the form:

1.5 Weir equation solution Numerical prediction by current model

Le=-6.9%

1

1 Le=-3.8%

Le=-5.2%

Le=-1.2%

0.5

0.5

Le=-0.5% Le=2.3% Le=-2.6% Le=-10.2%

0

0

0.1

0.2

0.3 0.4 0.5 Overflow depth (m)

0.6

0.7

0.8

0

Fig. 6. Comparison between the numerical prediction and the weir equation for the case of overflow without waves. The error was defined as Le ¼ ðqpredicted  qweir Þ=qweir  100ð%Þ, where qweir was calculated with Eq. (17).

Table 2 Grid convergence study for the case of overflow without waves. The error was defined as Le ¼ ðqpredicted  qweir Þ=qweir  100ð%Þ. And the relative error was defined as Lr ¼ ðLe c  Le f Þ=Le f , where Le c and Le f are the errors on coarse grid and finest grid respectively. Run No.

Rc ðmÞ

1 2 3 4 5 6 7 8 9

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

Le ð%Þ=Lr Þð%Þ

Dh ¼ 0:03 m

Dh ¼ 0:02 m

Dh ¼ 0:01 m

Dh ¼ 0:005 m

0.0/NA 19.2/102.6 8.7/255.1 6.7/185.6 4.0/594.8 6.1/384.9 8.7/125.6 9.8/86.5 10.2/46.0

0.0/NA 12.2/29.7 4.7/96.2 3.1/34.4 1.7/240.0 3.3/166.6 5.8/50.6 7.5/41.5 8.7/24.4

0.0/NA 10.2/7.4 2.6/4.1 2.3/1.0 0.6/4.7 1.2/0.1 3.8/0.0 5.3/1.9 6.9/0.1

0.0/NA 9.4/NA 2.4/NA 2.3/NA 0.5/NA 1.2/NA 3.8/NA 5.2/NA 6.9/NA

device to control and regulate open channel flows. A large number of theoretical and experimental studies have been carried out to determine weir characteristics. Among many others, Rouse and Reid [2] made an analytical investigation of the design of spillway crests based on the investigation of sharp-crested weir flow characteristics. Kandaswamy and Rouse [19] experimentally investigated the weir discharge coefficient C d as a function of H=P, where H is the driving head and P is the height of sharp-crested weir. This relationship was further studied by Ramamurthy et al. [20] based on experimental results and simplified theoretical considerations. Kindsvater and Carter [22] presented a comprehensive solution for the weir discharge characteristics based on experimental results and dimensional analysis. Rajaratnam and Muralidhar [21] determined the detailed distributions of velocity and pressure in the region of the weir crest using experimental methods. In a more recent work by Qu et al. [3], this classic problem was studied numerically using a two dimensional Reynolds Averaged Navier– Stokes (RANS) solver. The Froude number (F r ¼ pVffiffiffiffi, where V is the characteristic gL

velocity, g the acceleration due to gravity and L the characteristic length), the Reynolds number (Re ¼ qlVL, where q is the water den2

sity, and l the water viscosity), and the Weber number (W e ¼ qVr L, where r is the surface tension coefficient) are the key flow parameters for this problem. In this study, however, the Froude number is in a very narrow range, and the Reynolds number is very high (>106) to reduce viscous effects to negligible values, plus that the Weber number is high enough to ignore surface tension effects as H > 6:5 cm. Therefore, H=P is the only significant parameter characterizing the flow [4].

1328

X. Lv et al. / Advances in Water Resources 34 (2011) 1320–1334

The computational domain is shown in Fig. 7. The channel upstream of the weir (ab) is 6.0 m long in the x-direction. The channel downstream of the weir (df) is 1.0 m long. The weir is 0.297 m high (bc) and 0.005 m wide (bd). The flow domain is meshed with 350,256 cells (average Dh ¼ 0:002 mÞ and decomposed into 16 sub-domains to perform parallel computations. Locally refined prism layers are applied well within the boundary layers on all of the solid walls. The model execution time was around 4 h using 16 CPUs. A grid convergence study was also performed using a coarser grid with double the grid spacing as well as a finer grid with half the grid spacing. The results of the coarser grid size were in less agreement than the results using the other two grids. Further, the results for the finer grid were almost the same as those of the final grid chosen. The deviations of parameters (velocities and pressures), were generally much less than 1.5% between the results obtained from the final grid chosen and the finer grid. In order to better resolve the sharp geometry at the top of the weir, the grid was locally refined within this region, see Fig. 8 for example. In Figs. 9 and 10 the simulation results of the instantaneous flow evolution are plotted. The colours in Fig. 9 denote the pressure field. As can be seen from the snapshots, the water flows over the weir due to gravity, and after about 0.3 s the jet front impacts on the right bottom of the tank with large deformations. As in the broad-crested weir simulation, the fluid is fully driven by gravity and inertial forces, and the pressure distribution exhibits a similar pattern. Fig. 10 shows the streamlines distribution at later stages of the simulation, where the details of the recirculation flow formed by the backward-moving overflow discharge can be easily identified. The current numerical predictions were validated quantitively against the full set of experimental data related to H/P = 0.625 [21,20], as well as the numerical results by Qu et al. [3]. The validation included comparisons for the water surface profile, pressure profiles, and velocity distributions at two different sections (t-c and s0 -s). Due to the same reasons explained in Qu et al. [3], the velocity readings based on Pitot tube measurements were used for the model validation. Fig. 11 shows the distributions of pressure and velocity as well as the velocity angle u ¼ tan1 v =u. In Fig. 11(a), Y e is nappe thickness at crest ‘c0 (see Fig. 7) and y is distance above the crest. The non-dimensional pressure distribution obtained by the present model is compared with the experimental data of Rajaratnam and Muralidhar [21], Ramamurthy et al. [20] and the numerical result from Qu [3]. In Fig. 11(b), the non-dimensional x-component velocity u=U 0 is plotted against the non-dimensional

Fig. 8. Schematic view of local mesh refinement near the top of the sharp-crested weir.

flow depth pffiffiffiffiffiffiffiffi above the sharp crest y=Y e . Here, the reference velocity U 0 ¼ 2gh. The predicted variation of the velocity angle u with y=Y e is given in Fig. 11(c). The results of the simulation are in generally good agreement with the measured data and show improvements against the numerical results by Qu et al. [3]. These improvements are ascribed to the additional accuracy of the numerical schemes adopted in this study. Fig. 12 shows the comparisons of the distributions of pressure and velocity as well as the velocity angle u at section s-s0 (see Fig. 7). Again, the predicted values are close to the measured data. In Fig. 13, the predicted surface profiles near the top of the weir are compared with the experimental profiles [21]. The two profiles agree extremely well which again confirms the accuracy of the current model. It is also interesting to note the differences in the free surface profiles for various freeboard scenarios. Fig. 14 shows the predicted free surface profiles for three freeboard values, H/ P = 0.625, 0.4 and 0.2 from top to bottom, respectively. The nappe impact location moves closer to the weir as freeboard decreases. Different vortical motion patterns are observed for the three values of freeboard.

Locations: 1

2 3 4 5

Ya

Fig. 7. Schematic view of computation domain for flow past sharp-crested weirs ðslopeab ¼ 0Þ. Domain size: ab = 6.0 m, df = 1.0 m. Weir height bc = P = 0.297 m, width bd = 0.005 m. The basic boundary condition setup is also shown in the figure.

X. Lv et al. / Advances in Water Resources 34 (2011) 1320–1334

1329

Fig. 9. Snapshots of the free surface profile at the initial stages of the simulation. The pressure contour field is plotted. The numerical results refer to output time t = 0.15 s, 0.3 s, 0.495 s and 0.9 s respectively. The pressure field was nondimensionalized using the reference pressure pref ¼ qwater U 2ref .

1330

X. Lv et al. / Advances in Water Resources 34 (2011) 1320–1334

Fig. 10. Streamlines of the flow field at the later stages of the simulation. The flow field is colour contoured by velocity magnitude. The numerical results refer to output time t = 2.4 s, 2.415 s, 2.43 s and 2.49 s, 2.505 s, 2.52 s respectively. The velocity magnitude was normalised by the maximum value in the flow field.

While Fig. 15 shows the quantitative comparisons of turbulence kinetic energy and dissipation vertical distributions at location s-s0 for these different freeboard senarios. The turbulence plays an import role when the nappe impacts on the channel bed downstream the weir and generates complex vortical motions and air entrapment processes. LES is particularly suitable for capturing the 3-D unsteady turbulent coherent structure and was therefore used in this study. The main focus of the model validation was to reproduce the pressure, velocity vertical distributions and overflow profiles near the weir crest where turbulence is negligible. The choice of turbulence characterisation has little effect on the computed flow fields in this region, as demonstrated in Fig. 16 which shows almost identical predictions of pressure distributions at the weir crest t-t0 using both LES and k < epsilon > turbulence models. Fig. 17 shows the vertical distributions of velocity and turbulence characteristics at 5 different locations for H/P = 0.625. Fig. 17a indicate that the turbulence is generated at the weir crest near location 1 and extend upwards and increase in intensity further downstream from location 1 to 5.

5. Conclusion An unstructured Navier–Stokes turbulence solver together with a newly developed hybrid Level Set/VOF method has been applied to reproduce the flow at broad and sharp crested weirs in a rectangular channel. The predictions of the numerical model agree well with the experimental measurements and show improvements against existing numerical results by other researchers in terms of water surface profiles, distributions of the pressure and velocity and overflow discharge. Furthermore, the model is capable of resolving the time evolution of very complex vortical motions, air entrainment and impact pressure variations due to violent collision downstream of the weir crest. These processes are difficult to measure in both the laboratory and field, but are important to our understanding of structure damage, water quality downstream the weir and scour at the channel bed near the weir. Our model results indicate that the characteristics of the vortical motions generated by the impact of the nappe and turbulence are dependent on the weir type and freeboard. These behaviours have not been compared to traditional

1331

X. Lv et al. / Advances in Water Resources 34 (2011) 1320–1334

(a)

1

(a) 0.9

Exp. (Ramamurthy et al., 1987) Exp. (Rajaratnam & Muralidhar, 1971) Num. (Qu et al., 2009)

0.8

0.9 0.8

Present Simulation

Present Simulation

0.6

0.6

0.5 0.4

0.3

0.2

0.2

0.1

0.1

0

0.2

0.4

p / ( γ Ye )

0.6

0.8

0

1

(b)

0.9

Exp. (Ramamurthy et al., 1987)

0.8

Exp. (Rajaratnam & Muralidhar, 1971) Num. (Qu et al., 2009)

0.4

p / ( γ Y0 )

0.6

0.8

1 Exp. (Ramamurthy et al., 1987) Exp. (Rajaratnam & Muralidhar, 1971) Num. (Qu et al., 2009) Present Simulation

0.7

0.6

0.6

y / Y0

y / Ye

0.2

0.8

Present Simulation

0.5 0.4

0.5 0.4

0.3

0.3

0.2

0.2

0.1

(c)

0

0.9

0.7

0

0.5 0.4

0.3

(b)

Hydrostatic Pressure

0.7

y / Y0

y / Ye

0.7

0

Exp. (Ramamurthy et al., 1987) Exp. (Rajaratnam & Muralidhar, 1971) Num. (Qu et al., 2009)

0.1 0.4

0.6

u / (U 0 )

0.8

1

1

0 0.4

0.6

0.8 u / (U 0 )

1

(c) 1

0.9

Exp. (Rajaratnam & Muralidhar, 1971) Num. (Qu et al., 2009)

0.8

Exp. (Rajaratnam & Muralidhar, 1971) Num. (Qu et al., 2009)

Present Simulation

0.8

Present Simulation

0.6

0.6

y / Y0

y / Ye

0.7

0.5 0.4

0.4

0.3 0.2

0.2 0.1 0

20

ϕ = tan v / u −1

40

60

Fig. 11. Quantitative comparisons at location c-t: (a) pressure distribution, (b) veloc ity distribution, and (c) angle of velocity.

0 -60

-40

-20

0

ϕ = tan −1 v / u Fig. 12. Quantitative comparisons at location s-s0 : (a) pressure distribution, (b) velo city distribution, and (c) angle of velocity.

1332

X. Lv et al. / Advances in Water Resources 34 (2011) 1320–1334

Fig. 13. Comparison between predicted free surface profiles against measured data near the weir crest.

(a)

1

y / Y0

0.8

H / P = 0.625 H / P = 0.4 H / P = 0.2

0.6

0.4

0.2

0

0.04

0.06

0.08

0.1

0.12

0.14

0.16

0.18

k

(b)

1

y / Y0

0.8

H / P = 0.625 H / P = 0.4 H / P = 0.2

0.6

0.4 Fig. 14. Snapshots of the free surface profile at the output time t = 0.9 s for three freeboards. In the above pictures, H/P equals to 0.625, 0.4 and 0.2 from top to bottom, respectively.

experimental methods in predicting the flow characteristics, numerical simulation of the weir flows based on extensively validated model provides an efficient while accurate alternative for various flow configurations encountered in practical engineering problems. In this work, the numerical model validation is limited to the weir crest region where the turbulence effect is negligible.

0.2

0

0

1

2

ε

3

4

5

Fig. 15. Quantitative comparisons at location s-s0 : (a) vertical turbulence kinetic energy distribution and (b) vertical turbulence dissipation distribution, for three different freeboards.

X. Lv et al. / Advances in Water Resources 34 (2011) 1320–1334

We expect that LES and Kappa-Epsilon turbulence models may produce different results near the nappe impact location, fo rinstance in their predictions of the air entrapment rate; and is the subject of ongoing research.

1 Exp. (Ramamurthy et al., 1987) Exp. (Rajaratnam &Muralidhar, 1971)

0.8

1333

Present Simulation (LES) Present Simulation (k-epsilon)

Acknowledgments

y / Ye

0.6

This research was supported by the Flood Risk from Extreme Events (FREE) Programme of the UK Natural Environment Research Council (NERC) (NE/E0002129/1), coordinated and monitored by Professors Chris Collier and Paul Hardaker. The numerical calculations have been carried out on the HECToR, the UK’s national supercomputing services funded by UK research agencies and the HPC facility at University of Plymouth.

0.4

0.2

References 0

0

0.2

0.4

0.6

p / ( γ Ye )

0.8

Fig. 16. Quantitative comparison of pressure distributions computed at location t-t0 between the current LES model and a standard two-equation k  e turbulence model.

(a)

1

location 1 location 2 location 3 location 4 location 5

y / Y0

0.8 0.6 0.4 0.2 0

(b)

0

0.5

1

1.5

k

1

location 1 location 2 location 3 location 4 location 5

0.8

y / Y0

2

0.6 0.4 0.2 0

(c)

0

5

10

15

ε

20

30

35

1

locations: 1 2 3 4 5

0.8

y / Y0

25

0.6 0.4 0.2 0 0.5

1

1.5

2

2.5

3

U magnitude Fig. 17. Vertical distributions of velocity and turbulence characteristics at 5 different locations for H/P = 0.625. (Locations 1 to 5 are indicated in Fig. 7): (a) turbulence kinetic energy; (b) dissipation distribution; (c) velocity magnitude.

[1] Emiroglu ME, Baylar A. Experimental study of the influence of different weir types on the rate of air entrainment. Water Qual Res J Can 2003;4:769–83. [2] Rouse H, Reid L. Model research on spillway crests. Civil Eng ASCE 1935;5:10–5. [3] Bagheri S, Heidarpour M. Flow over rectangular sharp-crested weirs. Irrigat Sci 2010;28:173–9. [4] Qu J, Ramamurthy AS, Tadayon R, Chen Z. Numerical simulation of sharpcrested weir flows. Can J Civil Eng 2009;36(9,1):1530–4. 5. [5] Ferrari A. SPH simulation of free surface flow over a sharp-crested weir. Adv Water Resour 2010;33(3):270–6. [6] Lv X, Zou QP, Reeve DE, Zhao Y. A novel coupled level set and volume of fluid method for sharp interface capturing on 3D tetrahedral grids. J Comput Phys 2009;229(7):2573–604. doi:10.1016/j.jcp.2009.12.005. [7] Lv X, Zhao Y, Huang XY, Su XH. An efficient parallel/unstructured-multigrid implicit method for simulating 3D fluid-structure interaction. Commun Comput Phys 2008;4(2):350–77. [8] Leonard BP. The ULTIMATE conservative difference scheme applied to unsteady one-dimensional advection. Comput Methods Appl Mech Eng 1991;88:17–74. [9] Zhao Y, Tai CH. Higher-order characteristics-based methods for incompressible flow computation on unstructured grids. AIAA J 2001;39(7):1280–7. [10] Zhao Y, Tan HH, Zhang BL. A high-resolution characteristics-based implicit dual time-stepping VOF method for free surface flow simulation on unstructured grids. J Comput Phys 2002;183:233–73. [11] Chorin AJ. A numerical method for solving incompressible viscous flow problems. J Comput Phys 1967;2:12–26. [12] Turkel E. Preconditioning methods for solving the incompressible and low speed compressible equations. J Comput Phys 1987;72(2):227–98. [13] Brackbill JU, Kothe DB, Zemach C. A continuum method for modelling surface tension. J Comput Phys 1992;100:335. [14] Yang X, James A, Lowengrub J, Zheng X, Cristini V. An adaptive coupled levelset/volume-of-fluid interface capturing method for unstructured triangular grids. J Comput Phys 2006;214:41–54. [15] Lin PZ, Liu, Philip L-F. A numerical study of breaking waves in the surf zone. J Fluid Mech 1998;359:239–64. [16] Sussman M, Smith KM, Hussaini MY, Ohta M, Zhi-Wei R. A sharp interface method for incompressible two-phase flows. J Comput Phys 2007;221:469–505. [17] Lohner R, Yang C, Onate E. On the simulation of flows with violent free surface motion. In: 44th Aerospace Sciences Meeting and Exhibit, Reno, NV, 2006. [18] Chadwick A, Morfett J. Hydraulics in civil and environmental engineering. London and New York: E & FN SPON; 1998. [19] Kandaswamy AM, Rouse H. Characteristics of flow over terminal weirs and sills. J Hydraul Div 1957;83(4):1–13. [20] Ramamurthy AS, Tim US, Rao MV. Flow over sharp-crested weirs. J Irrigat Drainage Eng 1987;113(2):163–72. doi:10.1061/(ASCE)0733-943, 113:2(163). [21] Rajaratnam N, Muralidhar D. Pressure and velocity distribution for sharpcrested weirs. J Hydraul Res IAHR 1971;9(2):241–8. [22] Kindsvater CE, Carter R. Discharge characteristics of rectangular thin-plate weirs. J Hydraul Div 1957;83(3):1–36. [23] Nichols BD, Hirt CW, Hotchkiss RS. SOLA-VOF: A Solution Algorithm for Transient Fluid Flow with Multiple Free Boundaries, Technical Report, LA8355, Los Alamos National Laboratory, 1980. [24] Cummins SJ, Francois MM, Kothe DB. Estimating curvature from volume fractions. Comput Struct 2005;83:425–34. [25] Francois MM, Cummins SJ, Dendy ED, et al. A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework. J Comput Phys 2006;213:141–73. [26] Afkhami S, Bussmann M. Height functions for applying contact angles to 2D VOF simulations. Int J Numer Methods Fluids 2008;57:453–72. [27] Sussman M, Smereka P, Osher S. A level set approach for computing solutions to incompressible 2-phase flow. J Comput Phys 1994;114(1): 146–59.

1334

X. Lv et al. / Advances in Water Resources 34 (2011) 1320–1334

[28] Rider J, Kothe DB. Stretching and tearing interface tracking methods. In: 12th AIAA Computational Fluid Dynamics Conference, June 20, 1995, San Diego, Paper Number AIAA-95-1717 or LA-UR-95-1145, 1995. [29] Fedkiw R, Aslam T, Merriman B, Osher S. A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method). J Comput Phys 1999;152(2):457–92.

[30] Liu TG, Khoo BC, Wang CW. The ghost fluid method for compressible gas– water simulation. J Comput Phys 2005;204:193–221. [31] Kang M, Fedkiw R, Liu X-D. A boundary condition capturing method for multiphase incompressible flow. J Sci Comput 2000;15:323–60. [32] Lv X, Zou QP, Zhao Y, Reeve DE. A Preconditioned implicit free-surface capture scheme for large density ratio on tetrahedral grids. Commun Comput Phys 2011. in Press.