A high-resolution shallow water model using unstructured quadrilateral grids

Computers & Fluids 68 (2012) 16–28

Contents lists available at SciVerse ScienceDirect

Computers & Fluids 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 / c o m p fl u i d

A high-resolution shallow water model using unstructured quadrilateral grids Soumendra Nath Kuiry a,⇑, Dhrubajyoti Sen b, Yan Ding a a b

National Center for Computational Hydroscience and Engineering, The University of Mississippi, MS 38677, USA Dept. of Civil Engineering, Indian Institute of Technology, Kharagpur, WB 721 302, India

a r t i c l e

i n f o

Article history: Received 18 July 2011 Received in revised form 9 April 2012 Accepted 21 July 2012 Available online 10 August 2012 Keywords: Shallow water equations Numerical models Upwind schemes Source terms Data reconstruction Slope limiter

a b s t r a c t A two-dimensional cell-centred finite volume model for quadrilateral grids is presented. The solution methodology of the depth-averaged shallow water equations is based upon a Godunov-type upwind finite volume formulation, whereby the inviscid fluxes of the system of equations are obtained using the HLL Riemann solver. A simple yet precise analytical expression is presented to compute hydrostatic flux through an interface of a quadrilateral cell in order to achieve exact balance between flux gradient and bed slope source terms under still water condition. A multidimensional gradient reconstruction procedure and a continuously differentiable multidimensional slope limiter based on a wide computational stencil are proposed to maintain second-order spatial accuracy. The proposed second-order scheme is shown to be more accurate even when distorted grids are used and is therefore more suitable for practical applications. The presented model is verified and validated by solving a wide variety of test cases having analytical solutions and laboratory measurements. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction The two-dimensional (2D) shallow water equations are numerically solved to simulate many hydraulic and environmental engineering flow problems. A great amount of literature exists describing various computational techniques such as the finite difference method [1,2], the finite element method [3,4] and the finite volume method [5–7], which have been developed to obtain satisfactory solutions of the system of equations. The use of the finite volume technique has become more popular in recent times for simulating free surface flows because of its simplicity of implementation and good flexibility for space discretization. In addition, the finite volume method is based on the integral form of the conservation equations, and thus a scheme in conservation form can easily be constructed to capture shock waves (shock-capturing property). Godunov-type finite volume solvers of the shallow water equations have a shock-capturing property that is essential to preserve discontinuous or steeply varying gradients that occur in transcritical and sharp-fronted shallow flows. Typical examples of Godunov-type flow models can be found in literature [e.g. 5–11]. By upwinding the flux within the integral conservation form of the governing equations, Godunov-type methods represent physically correct propagation of information throughout the flow field by solving sets of Riemann problems over the entire flow domain. ⇑ Corresponding author. E-mail addresses: [email protected] (S.N. Kuiry), [email protected] (D. Sen), [email protected] (Y. Ding). 0045-7930/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compfluid.2012.07.018

Many Riemann solvers now exist, out of which the HLL solver is preferred here because it better describes the flux for dry bed situation and does not require any entropy fix [8,12]. Although, the HLL Riemann solver is robust, difficulties arise in solving the Riemann problem when source terms are included in the governing equations. Essentially, a numerical imbalance is created by artificially splitting the surface gradient into flux gradient and bed slope terms. Improper discretization of these terms does not result in a well-balanced scheme. Bermúdez and Vázquez-Cendón [13] proposed upwinding of the source terms to achieve equilibrium between flux gradients and source terms, and this method significantly improved the accuracy of the numerical solution, compared with earlier methods. Later, Vazquez-Cendon [14] modified the same idea to solve more general flow problems in the case of a one-dimensional channel with longitudinal width variations. LeVeque [15] proposed a well-balanced scheme by introducing a Riemann problem inside a cell to account for the propagation of source terms, but this has been reported to be less robust when predicting steady transcritical flows that contain shocks. Burguete and Garcia-Navarro [16] proposed conservative schemes with flux adjusted source term discretization technique using either a semi-implicit or upwinding method. Zhou et al. [10] developed the Surface Gradient Method (SGM) and pointed out that the main source of error is caused by inaccurate reconstruction of water depth. Valiani and Begnudelli [17] presented a novel method to compute the bed slope source terms from the pressure terms evaluated for all the faces of an n-sided cell. Another simple balancing approach was presented by Kuiry et al. [11,18] for triangular and Cartesian grids by introducing an improved treatment of pressure

17

S.N. Kuiry et al. / Computers & Fluids 68 (2012) 16–28

term. However, most of the literature on well-balanced schemes are limited to Cartesian and triangular grids. The development of well-balanced scheme for unstructured quadrilateral grids is relatively difficult and hence the simple analytical method of Kuiry et al. [11,18] has been modified and presented herein. The performance of Riemann based approaches greatly depends on temporal and spatial accuracy. The use of piecewise constant data for the first-order accurate schemes often does not fulfil the desired accuracy. Therefore, higher-order schemes are developed by reconstructing the conservative variables within a cell. The reconstruction procedure involves computations of gradients of the variables based on the cell-centred values in a defined stencil. While generating grids using a mesh generator distorted grids in some regions may be automatically generated. Nevertheless, a grid can be further distorted due to the presence of steep gradients in the bottom topography. The MUSCL based one-dimensional reconstruction approach [19] has been quite effective for structured grids but may produce poor results on distorted grids [20–22]. An alternative approach could be the use of multidimensional reconstruction method based on the cell-centred and cell-vertex formulations on unstructured grids. Such a reconstruction method possesses dependence on a wide computational stencil and does not strongly depend on vertex values to preserve stability even for distorted grids. However, the higher-order schemes often produce nonphysical oscillations. These oscillations can be suppressed by limiting the slopes of the reconstructed variables using a nonlinear function called a limiter. For unstructured grids, the limiter should be inherently multidimensional and a one-dimensional approach may not be suitable. Jawahar and Kamath [21] and Yoon and Kang [22] proposed multi-dimensional reconstruction procedure and multidimensional limiter for triangular grids but similar developments for quadrilateral grids are not found in literature. In the present study, a gradient reconstruction method and a multidimensional limiter suitable for regular and distorted quadrilateral grids are presented. The proposed second-order method is shown to be producing better results when distorted grids are used. The study presents a cell-centred finite volume model on quadrilateral grids. In order to achieve exact balance between flux gradients and source terms, a simple algebraic method is proposed. For higher-order accuracy, a multidimensional reconstruction technique and a continuously differentiable multidimensional limiter are introduced for unstructured quadrilateral grids. The model is applied to a number of analytical and experimental test problems and the computed results are investigated for comparative study. 2. Governing equations The two-dimensional, shallow water mathematical model is obtained by integrating the Navier–Stokes equations over the flow depth. The assumptions made in the process are: incompressible fluid, uniform velocity distribution in the vertical direction, hydrostatic pressure distribution and small bottom slope. Neglecting diffusion of momentum due to viscosity and turbulence, wind effects and the Coriolis term, the continuity and momentum equations in conservation form can be expressed as [23]:

@U @FðUÞ @GðUÞ þ þ ¼ Sðx; y; UÞ @t @x @y in which

U ¼ ðH; qx ; qy ÞT

F¼

2 q2 gh qx qy qx ; x þ ; h 2 h

!T

ð1Þ

G¼

qy ;

qx qy q2y gh2 ; þ h h 2

!T

 T S ¼ 0; ghðS0x  Sfx Þ; ghðS0y  Sfy Þ where U represents the vector of conserved variables, F and G are the fluxes associated with the conserved variables in the x- and y-directions. In addition, qx = uh and qy = vh are the unitary water discharges, u and v are depth-averaged velocities in the x- and ydirections, respectively, while Sox = @zb/@x and Soy = @zb/@y are the bed slopes along these directions. The variable H = h + zb is the water level, h represents depth of flow, zb defines bottom elevation, g is the acceleration due to gravity. Here, water level (H) is considered as one of the unknown variables instead of water depth (h) due to the fact that reconstruction of h for higher-order spatial accuracy directly from the cell average values will not guarantee a continuous reconstruction at the cell boundaries if the bed slope varies from cell to cell [10,24]. The friction slopes are estimated using the Manning relation [1,5,25] as given below.

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4=3 Sfx ¼ n2 u u2 þ v 2 h ;

Sfy ¼ n2 v

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4=3 u2 þ v 2 h

ð2Þ

where n is the Manning’s roughness coefficient. In general, the influence of bottom roughness prevails over the turbulent shear stress between cells. Therefore the effective stress terms are neglected in the computations. 3. Numerical solution The computational domain is divided into a finite number of unstructured quadrilateral cells which form the control volumes. Eq. (1) is then integrated over an elementary control volume and discretized by finite volume method. The dependent variables of the system are assumed to be stored at the centre of the cell and represented as piecewise constants. It is useful to rewrite Eq. (1) as

@U þ r  EðUÞ ¼ Sðx; y; UÞ @t

ð3Þ

where E = (F, G)T, the flux tensor, is introduced in order to express the integral form of the equations over a fixed volume V,

@ @t

Z V

UdV þ

Z

ðr  EÞdV ¼

V

Z

SdV

ð4Þ

V

For the computational domain defined by a set of quadrilateral cells, a discrete approximation to Eq. (4) is applied over each cell V so that the volume integrals represent integrals over the cell with the dependent variables represented as piecewise constants. The Gauss divergence theorem is applied to the second integral of Eq. (4) and the contour integral is approached via a mid-point rule, that is, a numerical flux is defined at the mid-point of each edge giving

I n

ðFnx þ Gny Þdn ¼

4 X Ee  ne dne

ð5Þ

e¼1

where e is the cell edge index and E⁄ is the numerical flux vector, n is the boundary of the plan area A of the control volume V, nx and ny are the components of unit normal in the x- and y-directions, respectively. The explicit expression of E⁄ depends upon the selected Riemann solver [5,7,12,26–31]. In the present work, the HLL Riemann solver of Harten, Lax, and van Leer [12,28] is used to compute the flux. It is preferred to Roe’s Riemann solver because it better describes the flux for dry bed situation and does not require any entropy fix [8,12]. The HLL scheme assumes one constant state between left and right waves and defines the flux at an interface as

18

S.N. Kuiry et al. / Computers & Fluids 68 (2012) 16–28

E  n ¼

8 > < ðEL Þ  n

if SL P 0

SR ðEL ÞnSL ðER ÞnþSL SR ðUR UL Þ SR SL >

:

ðER Þ  n

4.1. Second-order reconstruction

if SL < 0 < SR

ð6Þ

if SR P 0

in which EL = E(UL) and UL are the flux and conservative variable vectors evaluated at the left-hand side of each cell interface while ER = E(UR) and UR are the same quantities evaluated at the righthand side. The left hand side of an interface always refers to the inside of the cell. The symbols SL and SR represent the celerities of the waves, separating constant states of the local Riemann problem solution at cell interfaces. They can be estimated through the ‘‘two expansion’’ approach of Toro [32] as follows:

pffiffiffiffiffiffiffi pffiffiffiffiffiffiffi 8 n  ghL ; u  gh Þ if both sides are wet > < minðqL p ffiffiffiffiffiffiffi SL ¼ qL  n  ghL if the right side is dry > pffiffiffiffiffiffiffiffi : qR  n  2 ghR if the left side is dry ð7:aÞ pffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffi 8  ghR ; u þ gh Þ if both sides are wet > < maxðqR  n pffiffiffiffiffiffiffi SR ¼ qL  n þ 2 ghL if the right side is dry > pffiffiffiffiffiffiffiffi : if the left side is dry qR  n þ ghR ð7:bÞ

u ¼

1 ðq þ qR Þ þ 2 L

qffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffi ghL  ghR

qffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffi 1 1  ghL þ ghR þ ðqL  qR Þ gh ¼ 2 4

ð8Þ

ð9Þ

where q = [u, v]T. If the intermediate states UL and UR are defined as the cell-centred values, a first-order accurate scheme is obtained, which suffers from excessive numerical dissipation, and accuracy is undermined. A more accurate method can be obtained by reconstructing the variables for the left and right states, as described in the following sections.

Second-order accuracy is achieved by means of a multidimensional linear reconstruction of the cell variables. For a given cell with centroid i, for example, this requires the reconstruction of the cell variables in the form

Uðx; yÞ ¼ Ui þ rUli  ri

ð10Þ

where ri is the position vector relative to the cell centroid and rUli is the limited gradient of the variables. Fig. 1 shows a quadrilateral cell i surrounded by neighbours j, k, l and m. The gradients of U for the triangles at edge e = 1 formed by joining its end points 1 and 2 with the centroids of the corresponding cells i and j across the edge and denoted as (rU)1i2 and (rU)1j2, respectively, may be computed using the Green–Gauss theorem as under:

rU ¼

I

1 A

U  ndn

ð11Þ

n

where n is the integration path connecting the vertices of each triangle and A is the area of the triangle. The calculation of gradients (rU)1i2 and (rU)1j2 requires the values of conserved variables at the cell vertices. A second-order accurate method is employed to calculate the vertex values using the pseudo-Laplacian approach introduced by Holmes and Connell [36]. Due to one-sided triangulation at the cell boundaries, this interpolation method yields unrealistic values of weights that are either zero or negative. However, these unrealistic values can be modified simply by assigning unity for the weights if they all vanish and if any value becomes negative its absolute value is taken for simplicity and no weight is allowed to exceed unity. Such a modification has negligible effects on the solution and a detailed study on this topic can be found in [21]. 4.2. Unlimited gradient The unlimited gradient (rU)e=1 at the edge e = 1 common to cells i and j is computed, using the area-weighted average of the gradients of the two triangles across the edge formed by joining the edge vertices and the respective cell centroids having areas A1i2 and A1j2, in the following manner:

h

A1i2 ðrUÞ1i2 þ A1j2 ðrUÞ1j2   ¼ A1i2 þ A1j2

i

4. Linear reconstruction

ðrUÞe¼1

For a first-order accurate Riemann solver, the left and right states are nothing but the corresponding cell-centre values. This first-order accurate scheme, based on piecewise constant data within a cell is often inadequate to achieve desired accuracy. For this reason, implementation of a higher-order scheme for the variables inside each cell involving a gradient reconstruction is necessary. A number of higher-order schemes applied to shallow water equations on unstructured triangular and quadrilateral grids [7,9,22,33–35] are available. The comprehensive studies presented in [20–22] have shown that the higher-order schemes based on the extension of one-dimensional MUSCL approach to unstructured grids make the schemes strongly dependent on grid connectivity and therefore poor results are obtained on distorted grids. Therefore, MUSCL based second-order scheme may not be a good choice for unstructured quadrilateral grids when the presence of distorted cells in the flow domain cannot be avoided. In the present investigation, a gradient-reconstruction technique proposed by Jawahar and Kamath [21] for triangular cells, which includes a wide computational stencil and does not strongly depend on vertex values, is devised for the quadrilateral grids and discussed below.

The gradients at the other edges of the cell i, that is, (rU)e=2, (rU)e=3, (rU)e=4, are computed in a similar manner considering the neighbouring cells k, l and m of the other three edges,

u

A n

B

ð12Þ

H

l

t

G

m

1 C

p D

i

e=1 j q

e=2

o

e=4

e=3

s

k

2 r E

Fig. 1. Stencil for new reconstruction.

F

19

S.N. Kuiry et al. / Computers & Fluids 68 (2012) 16–28

respectively. Though viscous terms are not considered in this study, it can be noted that Eq. (12) can be used to evaluate such terms. Thus the reconstruction technique proposed here can reduce the computational cost by computing the gradients only once and using them for both inviscid and viscous fluxes. The unlimited cell gradient (rU)i for a quadrilateral cell i is then calculated using the area-weighted average of the corresponding face gradients given by   Ai1j2 ðrUÞe¼1 þ Ai2k3 ðrUÞe¼2 þ Ai3l4 ðrUÞe¼3 þ Ai4m1 ðrUÞe¼4   ðrUÞi ¼ Ai1j2 þ Ai2k3 þ Ai3l4 þ Ai4m1

ð13Þ In the above expression, the areas are the quadrangular regions formed by the two triangles on either side of each edge of the cell i. 5. Multidimensional limiter Unfortunately, higher-order schemes produce unphysical oscillations near discontinuities and some kind of gradient limiting procedure is required as a remedial measure. The construction of the limiter should be inherently multidimensional and a one-dimensional implementation on unstructured grids is not suitable [21,22]. Also, the limiter should be continuously differentiable, which helps in achieving a smooth transition at a discontinuity. The use of nondifferentiable functions such as max and min may adversely affect the convergence to steady state. The limited gradient ðrUlj Þ within a cell i can be obtained as weighted average of the four representative gradients of the surrounding cells j, k, l and m:

ð14Þ

where the weights xj, xk, xl and xm are given by the multidimensional function and rUj, (rU)k, (rU)l and (rU)m are the unlimited gradients of the four neighbouring cells as given by Eq. (13). Instead of the four unlimited face gradients chosen for limiting by Cueto et al. [33], an implementation more straightforward to the physics of finite volume formulation is adopted in the present study. Since each cell interacts with its neighbour through the common interface, the cell-centred gradients of four neighbouring cells should be more relevant choice and are used in the present work. The resulting stencil is shown in Fig. 1. The weights are prescribed as

xk ¼

xl ¼

g k g l g m þ e3 þ g 2k þ g 2l þ g 2m þ 4e3

ð15:aÞ

g j g l g m þ e3 2 2 g j þ g k þ g 2l þ g 2m þ 4e3

ð15:bÞ

g j g k g m þ e3 þ g 2k þ g 2l þ g 2m þ 4e3

ð15:cÞ

g 2j

g 2j

Z

2

Sfx dV ¼ n

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi!t u2 þ v 2 h

V

Z

Sfy dV ¼ n2

4=3

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi!t u2 þ v 2 h

V

4=3

ð16:aÞ

v tþ1

ð16:bÞ

where t is the time index. The bed slope treatment needs special attention as improper discretization leads to imbalance with hydrostatic flux under still water condition, that is, q = 0 as oH/ox = oH/ oy = 0. The present study discretizes the bed slope term in a simple way with a well-balanced scheme and is discussed below. In the absence of flows through the boundaries, the momentum equations reduce to

 Z I  1 2 gh nx dn ¼ g hSox dV n 2 V

ð17:aÞ

 Z I  1 2 gh ny dn ¼ g hSoy dV n 2 V

ð17:bÞ

This indicates that the sum of the hydrostatic pressure forces at the interfaces of a cell in the longitudinal direction is numerically equal to the weight component of the water contained within the cell in the same direction. As the four vertices of a quadrilateral cell may not lie on the same plane it is difficult to compute bed slope easily like a triangular cell. In this study, therefore, a quadrilateral cell is assumed to be subdivided into two triangular sub-cells, D123 and D134, as shown in Fig. 2, and an equivalent bed slope source term evaluated. The hydrostatic pressure force on a vertical face (e.g. 1220 10 ) of a triangular sub-cell is calculated considering the face to be trapezoidal

3

2 4

h3 h4

g j g k g l þ e3 xm ¼ 2 2 g j þ g k þ g 2l þ g 2m þ 4e3

utþ1

3

ð15:dÞ

where g j ¼ krUj k22 ; g k ¼ krUk k22 ; g l ¼ krUl k22 ; and g m ¼ krUm k22 ; and e is a small number introduced to avoid division by zero when all the four gradients vanish in regions of uniform flow. The limiter presented in Eq. (14) is an extension of the limiter proposed by Jawahar and Kamath [21]. However, their limiter was designed for triangular grids. After computing unlimited gradients for every cell using Eq. (13), the limited gradients are calculated with the help of Eq. (14). Finally, for the left and right states of each cell, the variables are reconstructed by Eq. (10).

1

h2

hc1

xj ¼

The source term consists of friction and bed slope. It is convenient to discretize the bed friction term semi-implicitly in order to make the diagonal dominant and subsequently to increase model stability:

h c2

rUli ¼ xj rUj þ xk rUk þ xl rUl þ xm rUm

6. Treatment of the source terms

h1

2 1

2

4

1 Fig. 2. Water prism used for bottom slope.

20

S.N. Kuiry et al. / Computers & Fluids 68 (2012) 16–28

p 2 2 with an equivalent water depth ðheq ¼ fðh1 þ þh1 h2 þ h2 Þ=3gÞ at each cell face [11]. The net hydrostatic force on the quadrilateral cell, say for the x-direction, can be expressed as the summation of forces on all the faces of the sub-cells D123 and D134 as

I

q

  I

I

1 2 1 2 1 2 gheq nx dn ¼ q gheq nx dn gheq Þnx dn þq ð 2 2 2 D123 D134 1 ¼ cfhc1 ½ðy2  y3 Þh1 þ ðy3  y1 Þh2 þ ðy1  y2 Þh3  2 þhc2 ½ðy3  y4 Þh1 þ ðy4  y1 Þh3 þ ðy1  y3 Þh4 g ð18Þ

where c = qg, hc1 = (h1 + h2 + h3)/3 and hc2 = (h1 + h3 + h4)/3 depths at the centroids of the sub-cells D123 and D134, respectively, {hk, k = 1, 4} and {yk, k = 1, 4} are, respectively, the water depths and y co-ordinates of the vertices. The water depths at the vertices of a quadrilateral cell are obtained by a second-order reconstruction of the water surface at the cell vertices {Hk, k = 1, 4} as described before and deducting the corresponding vertex elevations {zk, k = 1, 4} of the cell bottom. In the HLL solver, interface fluxes are calculated based on one of three possible formulae in Eq. (6), depending on the local wave speeds at the cell interface. However, because oH/ox = oH/oy = 0 and thus H = constant for a wet bed application, it is trivial to prove that the interface fluxes are nothing but the hydrostatic pressure forces at the interfaces no matter which calculation formula is used. Therefore, Eq. (18), when divided by q, represents the discretized form of the left hand side of Eq. (17.a). The component of weight of water contained within the prisms resting on the linear triangular planes conforming to the bed topography may be computed by multiplying the weight (qgAkhck) with the bed slope (Soxk) in that direction, for example, along the xdirection, as

Z

2 X hSox dV ¼ qg hck Soxk Ak

q g

V

ð19Þ

k¼1

The slope of a triangular sub-cell (Soxk) is calculated as the slope of a plane passing through the three vertices of the sub-cell. The bed elevations {zk, k = 1, 4} in the bed slope equations are substituted by water depths at the vertices {zk, = Hk  hk, k = 1, 4}. The substitution follows the weight component of the volume of water in a quadrilateral cell, as for example in the x-direction, as

Z

q g

V

1 hSox dV ¼ cfhc1 ½ðy2  y3 Þh1 þ ðy3  y1 Þh2 þ ðy1  y2 Þh3 2  ðDwx Þ123 þ hc2 ½ðy3  y4 Þh1 þ ðy4  y1 Þh3  þðy1  y3 Þh4  ðDwx Þ134 ð20Þ

The detailed derivation of the weight component for a triangular cell is described in [11]. Eq. (20), when divided by q, represents the discretized form of the right-hand side of Eq. (17.a). In Eq. (20), the terms (wx)123 and (wx)134 appear due to second-order spatial variation of the water surface within a cell [11]. A second-order approximation of the water surface, thus, captures this net difference of forces within a cell. The first-order approximation neglects this difference as it considers a horizontal water surface within a cell to evaluate the pressure force terms, though this difference may not always be negligible [37]. When these exact expressions for the hydrostatic pressure and bed slope terms are directly discretized, Eq. (20) reduces exactly to that given by Eq. (18) for x-momentum under still water condition as (wx)123 and (wx)134 vanish and a still water condition is guaranteed to be satisfied. A similar proof can be derived for the y-direction.

7. Time integration Eq. (1) can be solved explicitly at each time step with suitable initial and boundary conditions. Though a second- or higher-order accuracy can be followed for the time discretization in the present study the first-order forward Euler scheme is used and the discretized form of Eq. (1) for cell i can be written as

Unþ1 i

¼

Uni

" #n 4    Dt X  E  ne dne þ DtSni Ai e¼1 e

ð21Þ

i

where Un+1 and Un are the conserved variables at the new and old time steps, respectively, Dt is the time step, Ai is the projected area of the cell Vi on horizontal plane. The stability criterion adopted has followed the usual method in explicit finite volumes and the formulation presented in [38] is used in this study. 8. Boundary conditions When a face of an element coincides with the boundary of the flow domain or some physical boundary, it is necessary to solve a boundary Riemann problem [9,22,11]. The theory of characteristics provides sufficient information to establish a relation to find out the unknown variables at the boundary. In this study open and wall boundary conditions are implemented. The flow through open boundary can be subcritical or supercritical depending upon the local flow regime. One may refer to [11] for details of the implementation of boundary conditions. 9. Applications To verify the accuracy and illustrate its general applicability, the proposed model is applied to a wide variety of problems ranging from theoretical to experimental observations as discussed in the following sections. The unstructured quadrilateral grids are generated using freely available grid generation packages ‘‘AUTOMESH2D’’ developed by Zhao and Ma [39] and CCHE_Mesh3.0 mesh generator [40]. 9.1. Stationarity test The stationarity test over a 2D bump in a closed channel is conducted to test the accuracy of the proposed bed slope discretization technique. The spatial domain is represented by a 1  1m rectangular cross section channel, discretized into 3000 quadrilateral cells. The bump is defined by relatively more refined cells as shown in Fig. 3. The problem is quite similar to that solved in [41] but with a different peak value of the hump and flow variables. The channel

Fig. 3. Computational mesh for stationarity test.

S.N. Kuiry et al. / Computers & Fluids 68 (2012) 16–28

21

Fig. 5. Water surface comparison for subcritical flow over a hump. Fig. 4. Depth contour plot for stationarity test.

Table 1 Error measures for low- and high-resolution grids. Error

Dx = 0.05

Dx = 0.1

Dx = 0.125

Dx = 0.25

L1 L2 L1 RMSE

5.00475  1005 5.00475  1005 5.00475  1005 1.00090  1004

1.11765  1004 1.82502  1004 7.33282  1004 3.63086  1004

2.41179  1004 4.33201  1004 2.15015  1003 8.58499  1004

2.01449  1004 5.63337  1004 4.20026  1003 1.12151  1003

is frictionless and all four sides are closed by vertical walls. The bottom topography is described by the following function:

h

i zb ðx; yÞ ¼ max 0; 0:25  5  ðx  0:5Þ2 þ ðy  0:5Þ2 The model is run for sufficiently long time with the following initial conditions: H = 0.6 m, qx = 0.0 m2 s1 and qy = 0.0 m2 s1. The maximum deviations of the variables from the exact solutions are recorded as a measure of the artificial acceleration imposed by the scheme. The maximum deviations of the variables from the initial values are observed to be in the ranges of: H = 0.0, qx  1014 and qy  1014. Fig. 4 shows the depth contour plot obtained after running the model for sufficiently long time. As expected, the scheme preserves stationarity when the exact expressions for the bed slope and hydrostatic pressure force terms as proposed by the authors are utilized. 9.2. Grid convergence and order of accuracy test Evaluation of truncation error in formulating a numerical model is carried out by Taylor’s series analysis. According to Roach [42], the error of a numerical model can be formulated as a function of Dpc:

Error ¼ fn  fa ¼ C Dpc þ H:O:T: where fn and fa are numerical and analytical solutions for the mathematical model respectively, pc is the order of accuracy (convergence) and C is the evidence of a model’s convergence with a rate pc. The error can be computed by varying mesh size or time step size. In the following analysis, the error of the numerical solutions with respect to mesh spacing Dx will be evaluated by reporting error measurements and the coefficients C and pc. To conduct a grid convergence test, it is convenient to test it under steady state condition. A series of tests were performed in the

workshop on dam-break wave simulations [43] to test the accuracy of a numerical scheme. The spatial domain is represented by a 25  1 m rectangular cross section channel. The channel is frictionless and bottom topography is described by the following function:

8 for x 6 8m > <0 zðxÞ ¼ 0:2  0:05ðx  10Þ2 ; for 8 < x < 12m > : 0 for x P 12m According to the initial and boundary conditions, the flow may be subcritical, transcritical with or without shock, or supercritical. The analytical solution for the test cases can be obtained using Bernoulli’s theorem. The subcritical flow condition has been chosen herein for the grid convergence test. A discharge per unit width of qx = 4.42 m2 s1 is imposed as the upstream boundary condition and water level, h + zb = 2 m is specified at the downstream boundary. The model is run for sufficiently long time so that steady state is attained. The computed and analytical solutions at steady state are used for error analysis. For the convergence and order of accuracy test, four different grids with the mesh sizes Dx = 0.05, 0.1, 0.125 and 0.25 m are considered and different errors such as L1, L2, L1 and root mean square error (RMSE) as defined in [44] for the water surface are evaluated and reported in Table 1. The computed error limits show that the model is sufficiently accurate. The computed water surface on the coarse grid is compared with the analytical solution in Fig. 5. Jia [45] used RMSE and Ray and Nguyen [46] used L2 error for the convergence test. In this study both are evaluated and presented in Fig. 6. The values of the coefficients C, pc and regression coefficients are shown in Table 2. It can be concluded from Table 2 that the order of accuracy of the proposed model is slightly more than 1.54 though a second-order accurate scheme has been implemented. However, the order of accuracy may be improved by considering an ideal test case considering no bottom variation which maintains linear relationship in the numerical implementation so that error analysis by Taylor’s series is exactly valid. 9.3. Comparisons with experimental measurements by CERC The aim of this experimental test is to find out the performance of the proposed multidimensional limiter compare to the commonly used minmod and superbee limiters when the grids are distorted due to presence of steep bottom topography. The experiment conducted by Briggs et al. [47] consists of a two-dimensional

22

S.N. Kuiry et al. / Computers & Fluids 68 (2012) 16–28

Fig. 6. Order of accuracy by: (a) RMSE and (b) L2 error analysis.

sional limiter produce results with similar accuracy of minmod limiter but better than superbee limiter when the grids are distorted.

Table 2 Parameter coefficients for grid convergence and order of accuracy test. Error

C

pc

R2

L2 RMSE

0.00647 0.01279

1.5489 1.5458

0.8871 0.8884

problem in which a solitary wave passes around a conical island. The domain is 26 m long in the x-direction and 17.6 m long in the y-direction. The island is located at the centre of the domain and has a base diameter equal to 7.2 m and a top diameter of 2.2 m and is 0.625 m high. The still water depth over the flat portion of the domain is 0.32 m and bed friction is neglected. The domain is discretized by 300,000 quadrilateral cells with local refinement around the island. Due to the presence of the island with steep inclined surface many cells were distorted in this region. The mesh was generated using AUTOMESH-2D package. Such a huge number of cells were used in the computation in order to generate sufficiently accurate results at the observing gauge points. The details of the location of different gauges can be found in [48,49]. The solitary wave is generated at the inlet edge as suggested in Bradford and Sanders [48]. The water surface elevation was measured at several gauges around the island, but for the sake of brevity model results are compared with the measured data at stations G3, G6, G7 and G13 which are located closest to the shoreline. While simulating the wave propagation, three different limiters such as the proposed multidimensional, minmod and superbee were used. The computed resultant velocity at 12 s is analysed for the three simulations and are reported in Fig. 7. It can be seen that though the contour plots are almost similar, the proposed multi-dimensional limiter produces the smoothest results while the superbee limiter shows minor oscillations in the area marked by dotted circle in Fig. 7c. The results produces by minmod limiter are similar to that of multidimensional limiter. If the grid distortion is increased it is expected that the superbee limiter will produce more oscillatory results and a similar conclusion was made by Yoon and Kang [22] for their study on distorted triangular grids. The accuracy of the proposed model is investigated by comparing the time series of water surface elevation at different gauges in Fig. 8. Except at G6, at all other gauges the computed results are close to the measured data and the RMSE values at the four stations are 0.00537, 0.00706, 0.00377 and 0.00658 respectively. Similar results are reported by Kim et al. [48] and Bradford and Sanders [49]. The waves generated in the laboratory are dispersive and therefore the results can be further improved by considering dispersive terms which are not present in the governing shallow water equations. This test case proves that the proposed reconstruction method and the multidimen-

9.4. Supercritical flow in a channel with contraction Yoon and Kang [22] proved that the multidimensional gradient reconstruction and multidimensional limiter for triangular grids provide oscillation-free results when compared with the ‘‘superbee’’ limiter. A similar test case, in a symmetrical channel, is considered here but with quadrangular grids with distorted cells. Supercritical flow in a symmetrical channel is considered here. The initial width of 40 m is contracted from both sides with an angle of 15°. The 120 m long channel is straight after the contraction. The imposed inlet flow parameters are: Froude number, Fr = 3, and depth, h = 1 m. The steady state flow is simulated on three nonorthogonal meshes developed by using the CCHE_Mesh3.0 mesh generator developed at NCCHE [40] with 50, 60 and 65 grids in the cross-sectional direction and 150, 180 and 195 grids in the longitudinal direction, respectively. A small portion of the nonorthogonal mesh for 65  195 grids is shown in Fig. 9. It can be seen that the mesh at the contracted section is composed of distorted cells. Computations are carried out by using the present model, the ‘‘minmod’’ and the ‘‘superbee’’ limiters [50]. In all the simulations, the present model and the ‘‘minmod’’ limiter provided smooth solutions throughout the channel while the ‘‘superbee’’ limiter produced unphysical oscillations close to the first shock from the inlet. Fig. 10 shows the water surface and depth contour plots at steady state on 65  195 grids and remaining results are omitted due to space limitations. The ‘‘minmod’’ limiter produces almost same level of accuracy in this test case. However, it is expected that if the grid distortion is increased ‘‘superbee’’ limiter will produce more oscillations. This test case does not have analytical solutions but Cueto et al. [33] considered the same test and the results are very similar. This test case proves that the present model preserves accuracy even on distorted grids. 9.5. Partial dam-break in a laboratory flume An experimental partial dam-break test conducted by Fraccarolo and Toro [8] is considered in this section. The two-dimensional aspects of the flow, accuracy of the higher-order reconstruction and the ability of the model to correctly reproduce the flow behaviour can be validated by this test. The flume is rectangular in shape and connected to a tank at the upstream which can be treated as a reservoir. The flume which is considered as a floodplain, is open on all the three sides. The reservoir is 1 m long and 2 m wide while the floodplain is 3 m long and 2 m wide. The width of the gate which is

S.N. Kuiry et al. / Computers & Fluids 68 (2012) 16–28

23

Fig. 7. Velocity contour plot obtained using: (a) multidimensional limiter, (b) minmod limiter and (c) superbee limiter.

Fig. 8. Time histories of free surface displacements at different gauges: (a) G3, (b) G6, (c) G7 and (d) G13.

symmetrically centred, is 0.4 m. The gate is operated by a pneumatic cylinder and guarantees a very short opening time of less than 0.1 s. A number of gauges are located at different locations

to record the depth evolution using wave height metres and their locations are described in Table 3. Details of the experimental set-up and its operation can be found from the cited reference.

24

S.N. Kuiry et al. / Computers & Fluids 68 (2012) 16–28 Table 3 Locations of stage gauges. Stations

-5A

C

4

0

8A

x (m) y (m)

0.18 1.00

0.48 0.40

1.00 1.16

1.00 1.00

1.722 1.00

Fig. 9. Mesh for the flow through contracted channel (length is shown up to 54 m).

Fraccaorolo and Toro [8] conducted several tests out of which a specific case is selected in which the water level in the reservoir was kept at 0.6 m and the downstream floodplain was initially dry. Fig. 11 shows the computational mesh and the locations of the different gauges. The total number of quadrilaterals used is 6141 and finer cells are used in the vicinity of the gate where the flow pattern changes drastically. The time series of water depths at different gauges are compared in Fig. 12. In the graphs, the simulated results of Ying and Wang [51] are also presented to examine the accuracy of the present model. The time dependent curves are never monotonically decreasing but show oscillations. This means that though the water volume in the reservoir is strictly monotonically decreasing, the free surface is moving alternatively up and down. As soon as the gate is opened suddenly, a strong rarefaction wave is generated and moves backward into the reservoir while the front wave spreads laterally and flows onto the floodplain. The free surface variations at stations ‘‘-5A’’ and ‘‘C’’ compare well

Fig. 11. Computational mesh for the partial dam-break test and gauge locations.

with the experimental observations as well as the simulated results of [51] and the RMSE at these stations are 0.01757 and 0.01738. The gauges ‘‘0’’ and ‘‘4’’ shows some initial oscillations which are also confirmed by the experimental observations. At these stations, the water depth reduces drastically to a minimum soon after the gate is opened. After the minimum, the plots show a rising stage which is due to the contribution of the lateral propagation and represent, therefore, a two-dimensional effect. This aspect is well reproduced by the numerical model. The slight delay of the measured data can be attributed to the fact that actually the gate opening is not instantaneous as pointed out by Fraccarolo and Toro [8]. The RMSE values at these gauges are calculated as

Fig. 10. Water surface and depth contour plots: (a) minmod limiter, (b) superbee limiter and (c) present model.

S.N. Kuiry et al. / Computers & Fluids 68 (2012) 16–28

25

Fig. 12. Computed and measured depth evolutions at different gauges: (a) -5A, (b) C, (c) 4, (d) 0 and (e) 8A.

9.6. Dam-break flow in a sloped converging–diverging channel

Fig. 13. Converging–diverging section of the computational grid for the dam-break test case by Bellos et al. [50].

0.06386 and 0.05422 respectively. Nevertheless, the discrepancy at these stations can be attributed to the three-dimensionality of the flow which cannot be exactly captured by a depth-averaged twodimensional model. The depth prediction at station ‘‘8A’’ is close to the observed depth hydrograph and the wave arrival time to this station is quite accurate. The RMSE value at this station is 0.06571. The RMSE values are slightly large but the computed results are well compared with Ying and Wang (2006). In spite of the limitations of the shallow water representation, the overall predictive accuracy is quite satisfactory. Also, it can be noted that the proposed gradient reconstruction and limiter are accurate enough to capture discontinuous flow and the wet-dry situation.

Bellos et al. [52] conducted a series of dam-break test cases under various flow conditions. One of the experimental test cases is considered here. This test case can examine the accuracy of the model because it includes some difficulties such as an irregular flow domain and an initially dry bed. The width of the rectangular sloped channel varies like a converging–diverging channel. The detailed geometry of the test case can be found in [22,52]. The channel is 20.7 m long with uniform bottom slope 0.006. The Manning’s roughness co-efficient was reported as 0.012 m1/3 s. A gate is located at x = 8.5 m from the upstream and initial water level in the reservoir section, i.e. upstream of the gate, was set at 0.3 m. Depth evolutions were measured by eight probes located along the centre-line of the channel. All the three sides of the channel are defined as solid wall and the downstream is defined as freeoutfall. The simulation is performed on a grid of 5047 quadrilateral cells and the converging–diverging section is shown in Figs. 13 and 14 shows the measured and computed flow depths at the three locations at x = 0.0 m, x = 4.5 m and x = 18.5 m, respectively.

26

S.N. Kuiry et al. / Computers & Fluids 68 (2012) 16–28

Fig. 14. Computed and measured flow depth comparisons at: (a) x = 0.0 m, (b) x = 4.5 m and (c) x = 18.5 m.

Fig. 15. Definition sketch of the experimental set up for dam break over a triangular hump (distorted scale).

Fig. 14 shows that the present model reproduces the measured depths quite well. The computed RMSE at the three probes are within reasonable limits and are 0.00473, 0.00203 and 0.00766 respectively. In the downstream, the arrival time of the water front is accurately predicted. However, the present model slightly overestimates the flow depth at x = 18.5 m, shown in Fig. 14c. 9.7. Laboratory dam-break over a triangular hump The numerical model is applied to reproduce the complex laboratory dam-break flow over a triangular hump. The experiment was recommended by the EU CADAM project and conducted at the Recherches Hydrauliques Laboratory, Châtelet together with the University of Bruxelles (Belgium) under the supervision of J.M. Hiver. The experimental setup is illustrated in Fig. 15. A dam is located 15.5 m away from upstream end in a 38 m long rectangular channel. The dam has a reservoir in the upstream side with the still water surface elevation at 0.75 m and the downstream floodplain is initially dry. A triangular hump on the bed is located at 13 m away from the dam at the downstream. The symmetric triangular hump is 0.4 m high and the projected lengths of the normal and adverse slopes of the hump are both 3 m. The upstream end is solid wall and the downstream end is assumed to be open. A Manning coefficient 0.0125 m1/3 s1 is used throughout the domain. The simulation is run for 90 s on a grid with a resolution of 0.1 m size square cells. At t = 0, the dam suddenly breaks and the initial still water in the reservoir rushes onto the downstream floodplain. As soon as the wave front reaches to the hump, due to the interaction between the incoming flow and the bed topography a shock wave gets formed that propagates upstream. During this

time, a rarefaction wave is developed, moving downstream, which causes the water depth above the hump to decrease. The shock wave which was developed by the bed hits the solid wall at the upstream end at around t = 24 s and gets reflected. The reflection causes a second shock to form and propagate downstream. At the same time, the water on the hump starts receding and the drying-out process continues at the upstream side before the second shock hits the hump. Then the shock attacks the hump and the huge momentum carried by the shock pushes the water to climb up the hump so that the entire hump is again submerged. Another shock is immediately developed and starts to move upstream. These complex processes associated with wave interactions and wetting and drying continue until the momentum of flow is damped by friction effects and the loss of mass from the domain through the downstream boundary. Time histories of water depth recorded at five gauges located at 4 m (G4), 10 m (G10), 11 m (G11), 13 m (G13) and 20 m (G20) downstream of the dam (Fig. 15). The observed values are compared with the simulated results in Fig. 16. The RMSE at different gauges (G4: 0.02874, G10: 0.05266, G11: 0.04682, G13: 0.01702 and G20: 0.026488) are calculated for up to 65 s because after that results are very oscillatory. The maximum errors are contributed due to under prediction at the peaks. At all gauges, the wave arrival time is precisely predicted. However, computed results show that the time series at G20 is under predicted at all times and this is also predicted by other researchers [53] using different numerical schemes. Therefore, a possible explanation for the disagreement at G20 may be that the flow is highly complex and unstable beyond the hump and the shallow water equations based on hydrostatic assumption may not be quite suitable for this situation. However, the wave arrival time is accurately modelled, which is a more important factor for dam-break simulations. Overall, the comparison between the numerical predictions and measurements is satisfactory at most of the gauge points and the wet/dry transitions are simulated accurately by the current scheme. This confirms the effectiveness of the current method on flux and source term balancing and second-order implementation. 10. Conclusions A cell-centred finite volume model based on the HLL approximate Riemann solver is developed to solve the two-dimensional shallow water equations on unstructured quadrilateral grids. A

S.N. Kuiry et al. / Computers & Fluids 68 (2012) 16–28

27

Fig. 16. Time histories of water depth at different gauges: (a) G4, (b) G10, (c) G11, (d) G13 and (e) G20.

new discretization technique considered for the hydrostatic pressure term exactly satisfies the still water condition. The model is tested for quiescent flow condition with excellent results. The grid convergence and order of accuracy test shows that the multidimensional gradient reconstruction procedure and the continuously differentiable multidimensional limiter proposed for the model produces results close to second-order spatial accuracy. The viscous flux is computed as part of the gradient reconstruction and hence extra computations unlike other reconstruction methods are avoided. The proposed multidimensional limiter is proved to be better than the superbee limiter even when distorted grids are used. The model is further verified and validated against a number of analytical and experimental test problems reported by previous researchers including two-dimensional dam break tests and dam-break flow over a triangular hump. Based on the quantitatively satisfactory results of the test cases, it can be concluded that the proposed bed discretization, multidimensional gradient reconstruction procedure and limiter for unstructured quadrilateral grids make the shallow water model robust and capable of simulating various types of shallow water flows over wet/dry bed with complex topography and distorted grids. The proposed gradient-reconstruction technique for quadrilateral grids in combination with that for unstructured triangular grids would be a perfect tool to simulate complex flows like flows in rivers and floodplains, which may be spatially discretized arbitrarily according to the terrain description, possibly in conjunction with arbitrary triangular grids. It is also useful for viscous flow

computations on unstructured grids as the same computational stencil can be used for both the inviscid and viscous terms.

Acknowledgments The authors greatly acknowledge of Dr. S.K. Jahangir Hossain, IIT Kharagpur, India for his valuable suggestions. Dr. Ying is appreciated for supplying the required experimental and simulation data.

References [1] Fennema RJ, Chaudhry MH. Explicit methods for 2-D transient free-surface flows. J Hydraul Eng 1990;116:1013–34. [2] Molls T, Chaudhry MH. Depth-averaged open-channel flow model. J Hydraul Eng 1995;121:453–65. [3] Akanbi AA, Katapodes ND. Model for flood propagation on initially dry land. J Hydraul Eng 1988;114:689–706. [4] Tucciarelli T, Termini D. Finite-element modelling of floodplain flow. J Hydraul Eng 2000;126:416–24. [5] Zhao DH, Lai JS, Shen HW, Tabios III GQ, Tan WY. A finite-volume twodimensional unsteady flow model for river basins. J Hydraul Eng 1994;120:863–83. [6] Alcrudo F, Garcia-Navarro P. A high resolution Godunov-type scheme in finite volumes for the 2D shallow water equations. Int J Numer Methods Fluids 1993;16:489–505. [7] Anastasiou K, Chan CT. Solution of the 2D shallow water equations using the finite volume method on unstructured triangular meshes. Int J Numer Methods Fluids 1997;24:1225–45.

28

S.N. Kuiry et al. / Computers & Fluids 68 (2012) 16–28

[8] Fraccarollo L, Toro EF. Experimental and numerical assessment of the shallow water model for two-dimensional dam-break type problems. J Hydraul Res 1995;33:843–64. [9] Sleigh PA, Berzins M, Gaskell PH, Wright NG. An unstructured finite-volume algorithm for predicting flow in rivers and estuaries. Comput Fluid 1998;27:479–508. [10] Zhou JG, Causon DM, Mingham CG, Ingram DM. The surface gradient method for the treatment of source terms in the shallow-water equations. J Comput Phys 2001;168:1–25. [11] Kuiry SN, Pramanik K, Sen D. Finite volume model for the shallow water equations with improved treatment of the source terms. J Hydraul Eng 2008;134:231–42. [12] Toro EF. Riemann solvers and numerical methods for fluid dynamics. New York: Springer; 1997. [13] Bermúdez A, Vázquez-Cendón ME. Upwind methods for hyperbolic conservation laws with source terms. Comput Fluid 1994;23:1049–71. [14] Vazquez-Cendon ME. Improved treatment of source terms in upwind schemes for the shallow water equations. J Comput Phys 1999;148:497–526. [15] LeVeque RJ. Balancing source terms and flux gradients in high-resolution Godunov methods: the quasi-steady wave propagation algorithm. J Comput Phys 1998;146:346–65. [16] Burguete J, Garcia-Navarro P. Efficient construction of high-resolution TVD conservative schemes for equations with source terms: application to shallow water flows. Int J Numer Methods Fluids 2001;37:209–48. [17] Valiani A, Begnudelli L. Divergence form of bed slope source term in shallow water equations. J Hydraul Eng 2006;132:652–65. [18] Kuiry SN, Ding Y, Wang SSY. Two-dimensional modeling of coastal flows using shock-capturing finite volume method. Comput Fluids 2010;39: 2051–68. [19] van Leer B. Towards the ultimate conservative difference scheme, V. A second order sequel to Godunov’s method. J Comput Phys 1979;32:101–36. [20] Barth JT, Jespersen DC. The design and application of upwind schemes on unstructured meshes. AIAA paper 89-0366;1989. [21] Jawahar P, Kamath H. A high-resolution procedure for Euler and Navier-Stokes computations on unstructured grids. J Comput Phys 2000;164:165–203. [22] Yoon TH, Kang SK. Finite volume model for two-dimensional shallow water flows on unstructured grids. J Hydraul Eng 2004;130:678–88. [23] Cunge JA, Holly FM, Verwey A. Practical aspects of computational river hydraulics. London: Pitman Publishing Limited; 1980. [24] Nujic M. Efficient implementation of non-oscillatory schemes for the computation of free surface flow. J Hydraul Res 1995;33:101–11. [25] Chaudhry MH. Open channel flow. New Delhi: Prentice Hall of India Private Limited; 1994. [26] Roe PL. Approximate Riemann solvers, parameter vectors and difference schemes. J Comput Phys 1981;43:357–72. [27] Roe PL. A basis for upwind differencing of the two-dimensional unsteady Euler equations. Numer method fluid dynam II. Oxford University Press; 1986. [28] Harten A, Hyman JM. Self adjusting grid method for one-dimensional hyperbolic conservation laws. J Comput Phys 1983;50:235–96. [29] Alcrudo F, García-Navarro P, Saviron JM. Flux difference splitting for 1D open channel flow equations. Int J Numer Methods Fluids 1992;14:1009–18. [30] Valiani A, Caleffi V, Zanni A. Case study: Malpasset Dam-break simulation using two-dimensional finite volume method. J Hydraul Eng 2002;128:460–72. [31] Caleffi V, Valiani A, Zanni A. Finite volume method for simulating extreme flood events in natural channels. J Hydraul Res 2003;41:167–77. [32] Toro EF. Riemann problems and the WAF method for solving the 2dimensional shallow-water equations. Phil. Trans Roy Soc Lond Ser A: Phys Sci Eng 1992;338:43–68.

[33] Cueto-Felgueroso L, Colominas I, Navarrina JFF, Casteleiro M. High-order finite volume schemes on unstructured grids using moving least-squares reconstruction. Application to shallow water dynamics. Int J Numer Methods Eng 2006;65:295–331. [34] Castro Diaz MJ, Fernández-Nieto ED, Ferreiro AM, Parés C. Two-dimensional sediment transport models in shallow water equations. A second order finite volume approach on unstructured meshes. Comput Methods Appl Mech Eng 2009;198:2520–38. [35] Loukili Y, Soulaïmani A. Numerical tracking of shallow water waves by the unstructured finite volume WAF approximation. Int J Comput Methods Eng Mech 2007:1–14. [36] Holmes DG, Connell SD. Solution of the 2D Navier Stokes equations on unstructured adaptive grids. In: Proceedings of AIAA 9th CFD conference, AIAA paper 89-1932; 1989. [37] Kuiry SN, Sen DJ. Closure to finite volume model for shallow water equations with improved treatment of source terms by Soumendra Nath Kuiry, Kiran Pramanik, and Dhrubajyoti Sen. J Hydraul Eng 2009;135:1017. [38] Brufau P, García-Navarro P, Vazquez-Cendon ME. Zero mass error using unsteady wetting–drying conditions in shallow flows over dry irregular topography. Int J Numer Methods Fluids 2004;45:1047–82. [39] Zhao G, Ma X. A fully automatic adaptive quadrilateral mesh generator (AUTOMESH-2D). [accessed 01.07.10]. [40] Zhang Y, Jia Y. CCHE-MESH: 2D structured mesh generator user’s manual, version 3.0. NCCHE-TR-2009-01. [41] Brufau P, García-Navarro P. Unsteady free surface flow simulation over complex topography with a multidimensional upwind technique. J Comput Phys 2003;186:503–26. [42] Roache PJ. Verification and validation in computational science and Engineering. PO Box 9110, Albuquerque (New Mexico, USA): Hermosa Publishers; 1998. p. 87119–9110. [43] Goutal N, Maurel F. Technical report HE-43/97/016/A, Electricité de France, Département Laboratoire National d’hydraulique, Groupe Hydraulique Fluviale; 1997. [44] Smith PE. Analytical solutions for mathematical verification. USA: American Society of Civil Engineers; 2009. p. 45–57. [45] Jia Y. Method of manufactured solution (MMS). USA: American Society of Civil Engineers; 2009. p. 132–46. [46] Ray RK, Nguyen KD. An unstructured finite-volume technique for shallowwater flows with wetting and drying fronts. World Acad Sci Eng Tech 2010;71:682–7. [47] Briggs MJ, Synolakis CE, Harkins GS, Green DR. Laboratory experiments of Tsunami runup on a circular island. Pure Appl Geophys 1995;144:569–93. [48] Kim D-H, Cho Y-S, Yi Y-K. Propagation and run-up of nearshore tsunamis with HLLC approximate Riemann solver. Ocean Eng 2007;34:1164–73. [49] Bradford SF, Sanders BF. Finite-volume model for shallow-water flooding of arbitrary topography. J Hydraul Eng 2002;128:289–98. [50] Roberts R, Zoppu C. Robust and efficient solution of the 2d shallow water equation with domains containing dry beds. Austral Math Soc 2000:C1260–82. . [51] Ying X, Wang SSY. Modeling flood inundation due to dam and levee breach. In: Proceedings of the US–China workshop on advanced computational modelling in hydroscience & engineering, Oxford, Mississippi, USA; 19–21 September, 2006. [52] Bellos CV, Soulis JV, Sakkas JG. Experimental investigation of two-dimensional dam-break induced flows. J Hydraul Res 1992;30:47–63. [53] Brufau P, Vazquez-Cendon ME, García-Navarro P. A numerical model for the flooding and drying of irregular domains. Int J Numer Methods Fluids 2002;39:247–75.