A mass conservative scheme for simulating shallow flows over variable topographies using unstructured grid

Advances in Water Resources 28 (2005) 523–539 www.elsevier.com/locate/advwatres

A mass conservative scheme for simulating shallow flows over variable topographies using unstructured grid A. Mohamadian

a,b,*

, D.Y. Le Roux a, M. Tajrishi b, K. Mazaheri

c

a

De´partement de Mathe´matiques et de Statistique, Universite´ Laval, Que´bec, Canada, G1K 7P4 Department of Civil Engineering, Sharif University of Technology, P.O. Box 11365-8939, Tehran, Iran Department of Aerospace Engineering, Sharif University of Technology, P.O. Box 11635-8639, Tehran, Iran b

c

Received 7 August 2003; received in revised form 1 June 2004; accepted 22 October 2004

Abstract Most available numerical methods face problems, in the presence of variable topographies, due to the imbalance between the source and flux terms. Treatments for this problem generally work well for structured grids, but most of them are not directly applicable for unstructured grids. On the other hand, despite of their good performance for discontinuous flows, most available numerical schemes (such as HLL flux and ENO schemes) induce a high level of numerical diffusion in simulating recirculating flows. A numerical method for simulating shallow recirculating flows over a variable topography on unstructured grids is presented. This mass conservative approach can simulate different flow conditions including recirculating, transcritical and discontinuous flows over variable topographies without upwinding of source terms and with a low level of numerical diffusion. Different numerical tests cases are presented to show the performance of the scheme for some challenging problems.
2004 Published by Elsevier Ltd. Keywords: Shallow-water; Variable topography; Unstructured grid; Recirculating flow; Finite volume method; Mass conservative approach

1. Introduction Shallow recirculating flows are a wide category of fluvial flows. They occur in many hydraulic situations such as breakwaters, bridge piers, embayment and channel junctions. Modeling these flows is important for engineering applications and optimum design of hydraulic structures. On the other hand, for many fluvial flows, the flow regime changes from subcritical to supercritical and the employed numerical method should be able to solve sub, super and transcritical flows. In this paper, our primary interest is the simulation of the above-men-

* Corresponding author. Address: De´partement de Mathe´matiques et de Statistique, Universite´ Laval, Que´bec, Canada, G1K 7P4. Fax: +1 418 656 2817. E-mail address: [email protected] (A. Mohamadian).

0309-1708/$ - see front matter
2004 Published by Elsevier Ltd. doi:10.1016/j.advwatres.2004.10.006

tioned flows using unstructured grids. Unstructured grids are attractive because of their flexibility for representing irregular boundaries and for local mesh refinement. Extensive research has been performed in the area of shallow-water equations, and several upwind methods originally designed to solve the Euler equations have been extended to the shallow-water system. Those are for example, the RoeÕs method [10], the Beam-warming scheme [9], the monotonic upstream schemes for conservation laws (MUSCL) using curve-linear coordinate [4], the Osher and Salmon scheme [36], the essentially nonoscillatory (ENO) schemes [24] and the Harten, Lax and van Leer (HLL) solver [21]. Most of these methods have the capability of shock capturing with a high level of accuracy in few computational cells, and the flux vector is determined based on the wave propagation structure.

524

A. Mohamadian et al. / Advances in Water Resources 28 (2005) 523–539

Some of these methods perform well for particular types of flow like discontinuous or transcritical flows over flat topographies, but in the case of flows over variable topography and recirculating flows around structures, there is room for considerable improvement. Some modifications have been brought to the abovementioned methods in the case of variable topographies and structured grids. Bermudez and Vazquez [2] extended the van LeerÕs Q-scheme for variable topographies by using an upwind discretization of the source terms and they introduced the C property, which states that the scheme should preserve the stagnant conditions. Their work has been extended in the case of a general 1D channel with breath variation by Vazquez-Cendon [31]. Hubbard and Garcia-Navarro [14] then extended the Vasquez-Cendon approach to flux difference splitting schemes. However, upwinding of source terms is computationally expensive for practical applications because the source terms have to be projected on a basis of eigenvectors. LeVeque [20], introduced a Riemann problem inside a cell for balancing the source terms and the flux gradients, and the resulting method was found to preserve both stagnant and quasi steady state conditions. However, the LeVequeÕs scheme is not directly transportable to unstructured grids. Kurganov and Levy [18] extended their Central-Upwind scheme for the shallow-water equations and by using the water surface elevation instead of the depth, they proposed an adaptive algorithm for variable topographies. They also proved that their scheme preserves the steady state condition with a positive depth, without resorting to any artificial drying and wetting strategy. Unfortunately, their scheme poorly performs in the case of circulating flows, as shown in the following. Alcrudo and Benkhaldoun [3] defined the topography such that a sudden change in the topography occurs at the interface of two cells. They also developed a Riemann solver at the interface with a sudden change in the bed elevation. However, their approach leads to several cases of Riemann fan and it is numerically too expensive. Vukovic and Sopta [32] extended the ENO and weighted ENO (WENO) schemes to shallow-water equations with the source terms in 1D channels. Jin [16] developed the interface method, which preserves the steady state condition up to second order accuracy on structured grid, but his approach is not directly usable on unstructured grids. Xu [35] proposed a second order gas kinetic scheme for shallowwater flows over variable topographies. The gas kinetic schemes are basically different from characteristics based schemes and they should be tested for challenging test cases such as recirculating flows. Nujic [24] used the water level variable instead of the depth and he extracted the gravitational terms from flux functions. Then, using the Shu and Osher (SO) scheme [26], he computed the flux vector and obtained accurate results for dam break problems over variable topographies.

Unfortunately, the SO scheme [26] generates a high level of numerical diffusion in the case of recirculating flows, as shown in the following. Zhou et al. [38] introduced the surface gradient method and showed that interpolating the depth without considering the bed variations may lead to erroneous results. They used the interpolated water surface elevation to calculate the depth at the interface and they showed that this approach combined with the HLL flux function [12] satisfies the C property. Their scheme performs very well for variable topographies without any extra efforts for balancing the source terms and the flux gradients. However, the C property does not hold for unstructured grids and moreover, the HLL flux induces a high level of numerical viscosity in recirculating flows, as shown in the following. The objective of this paper is to present an efficient numerical method for recirculating flows over variable topographies on unstructured grids. In order to fulfill this goal, the NujicÕs method [24] is combined with the RoeÕs scheme, and the surface gradient method [38] is used for calculating the water depth at the interface. The resulting scheme is mass conservative, i.e. the convection terms are conserved but the gravity source term is discretized using a non-conservative method. Hence, it can simulate flows over variable topographies without upwinding of source terms. Several numerical tests are presented to show that the non-conservative discretization of the source terms induces a minor effect on the ability of the proposed scheme in simulating discontinuous flows. On the other hand, the mass conservation property ensures that the errors will not accumulate in time, which is crucial for the long-term simulations. This paper is organized as follows. The model equations and the numerical schemes are presented in Sections 2 and 3, respectively. In Section 4, some test case problems are numerically simulated, and they show the accuracy of the proposed method. Some concluding remarks complete the study.

2. Governing equations The shallow water equations (SW) in a conservative form are written as ~ d oU þ r~ F ¼~ S þ r~ F ; ot ~ and ~ ¼ ðh; uh; vhÞ, ~ with U F ¼ ð~ E; GÞ 0 0 1 uh vh B B C 2 2 ~¼@ ~ E ¼ @ u h þ 0:5gh A; G uvh uvh

2

v h þ 0:5gh

ð2:1Þ 1 2

C A;

ð2:2Þ

where h is the water depth, u and v are the velocity components (Fig. 1) and g is the gravitational acceleration.

A. Mohamadian et al. / Advances in Water Resources 28 (2005) 523–539

h

u

3. Numerical scheme

η

3.1. Finite volume methods on unstructured grids

z

In this paper, the variables are located at the geometric centers of the cells, and each triangle represents a control volume. Let A be the area of a triangle with boundary s. The SW equations are integrated over every control volume ! Z Z ~ oU d ~ ~ ~ þ rF  S  rF dA dt ¼ 0: ð3:1Þ ot t A

Fig. 1. Schematic diagram of an unsteady flow over an irregular bottom and the corresponding notation.

The superscript ÔdÕ refers to the diffusion and the diffusive flux has the following form: d d ~d Þ; ~ F ¼ ð~ E ;G

ð2:3Þ

where 0

0

1

d B C ~ E ¼ @ ðv þ vt Þh ou ox A; ðv þ vt Þh ov ox

0

0

1

ou C ~d ¼ B G @ ðv þ vt Þh oy A; ov ðv þ vt Þh oy

ð2:4Þ

and v and vt are water and eddy viscosity coefficients respectively. The source term ~ S is written as 0 1 0 pffiffiffiffiffiffiffiffiffiffiffiffiffiffi oz C B ~ ð2:5Þ S ¼ @ cf u u2 þ v2 þ gh ox A; pffiffiffiffiffiffiffiffiffiffiffiffiffiffi oz 2 2 cf v u þ v þ gh oy where zb is the distance between the bed surface and the reference level (Fig. 1) and cf is the friction coefficient. The SW equations may also be written in a non-conservative form ~ ~ ~ d oU oU oU þA þB ¼~ S þ r~ F ; ot ox oy

525

ð2:6Þ

where A and B are the Jacobian matrices 0 1 0 1 0 ~ oE B 2 C A¼ ¼ @ u þ c2 2u 0 A; ~ oU uv v u ð2:7Þ 0 1 0 0 1 ~ oG B C B¼ ¼ @ uv v u A; ~ oU v2 þ c2 0 2v pffiffiffiffiffi and c ¼ gh is the wave velocity. The eigenvalues ai and bi, (i = 1, 2, 3) corresponding to the matrices A and B respectively are a1 = u + c, a2 = u, a3 = u  c and b1 = v + c, b2 = v, b3 = v  c, respectively.

The time integration is performed by using the first order forward (explicit) Euler time-stepping scheme, which leads to ! Z ~ nþ1  U ~n d n U ~ ~ ~ ð3:2Þ þ ðrF  S  rF Þ dA ¼ 0: Dt A In (3.2) the accuracy of the time stepping scheme may be increased to second order by employing the Lax– Wendroff scheme as described in Section 3.5. By using the divergence theorem, the diffusive and convective flux integrals are transformed into boundary integrals Z I d d ~ ~ ðrF  rF ÞdA ¼ ð~ F ~ n~ F ~ nÞds; ð3:3Þ A

S

and the boundary integral is approximated by a summation over the triangle edges I 3 X d d ð~ F ~ n~ F ~ nÞds ¼ ð~ F k ~ nk  ~ F k ~ nk Þdsk : ð3:4Þ s

k¼1

The diffusive fluxes on the cell interfaces are usually approximated by a centered scheme d d d ~ F ¼ 0:5ð~ FR þ~ F L Þ;

ð3:5Þ ~ and the convective fluxes F in the Godunov type methods are calculated based on an exact or approximate Riemann solver. Most approximate Riemann solvers are written as

~ F ¼ 0:5ð~ FR þ~ F L  D~ F Þ;

ð3:6Þ

~ L Þ and ~ ~ R Þ are the left and F ðU FR ¼ ~ F ðU where ~ FL ¼ ~ right flux vectors. The subscripts ÔRÕ and ÔLÕ represent the evaluation of the right and left sides of the interface, respectively, and D~ F is the flux difference. When D~ F ¼ 0 in (3.6), the scheme is equivalent to a standard centered scheme. Hence, D~ F may be considered as an ‘‘artificial diffusive flux’’. The maximum possible (stable) artificial diffusive flux is employed in the Lax–Fredrichs scheme, which in the 1-D case is written as Dx ~ DU ; D~ F ¼ Dt ~ ¼ ðU ~R  U ~ L Þ. where DU

ð3:7Þ

526

A. Mohamadian et al. / Advances in Water Resources 28 (2005) 523–539

A less diffusive approximation for the artificial diffusive flux is that of the Rusanov scheme ~R  U ~ L Þ; D~ F ¼ jamax jðU

ð3:8Þ

where amax is the greatest eigenvalue of the Jacobian matrices corresponding to the left and right flux vectors in (3.6). The Rusanov scheme is the base of many central schemes, and it provides good results in the case of shock wave problems (see [19] for more details). How ever, D~ F in (3.8) excessively damps the solutions for recirculating flows. Another approach, which works well for both recirculating flows and shock wave problems, is the RoeÕs method, for which an approximate Riemann solver is obtained by solving a linearized system. Here, D~ F is calculated as e ~ ffi j Aj D~ F ¼ jAjDU

3 X

k

~ e~ ¼ ak~

k¼1

3 X

B u2 Þnx  ~ u~vny ¼ @ ð~c2  ~ 2 ~ u~vnx þ ð~c  ~v2 Þny

nx 2~ un þ ~vny ~vnx

1

ny ~ny u

C A;

~ unx þ 2~vny ð3:10Þ

where pffiffiffiffiffiffi pffiffiffiffiffi uR hR þ uL hL ~ u ¼ pffiffiffiffiffiffi pffiffiffiffiffi ; hR þ hL pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ~c ¼ gðhR þ hL Þ=2

pffiffiffiffiffiffi pffiffiffiffiffi v R hR þ v L hL ~v ¼ pffiffiffiffiffiffi pffiffiffiffiffi ; hR þ hL

ð3:11Þ

e are simply and the eigenvalues of A ~ unx þ ~vny þ ~c; a1 ¼ ~

~ a2 ¼ ~ unx þ ~vny ;

~ a3 ¼ ~ unx þ ~vny  ~c;

~v þ ~cny 1

~R  U ~ L Þ; D~ F ¼ jajðU

ð3:16Þ

where a = kjamaxj (with k 6 1), and a TVD Runge-Kutta time marching method.

In order to calculate the numerical fluxes in (3.5) and (3.6), the values of the variables at the left and right sides of the interface must be calculated. Zhou et al. [38] proposed the surface gradient method, in which the water surface elevation g is interpolated at the cell faces instead of the water depth. The water depth at the cell face is then calculated using the water surface elevation and the bed topography. For example, once gL is calculated, the water depth can be obtained as h L ¼ gL þ z e ;

ð3:17Þ

where ze is the distance between the bed surface and the reference level at triangle edge midpoints. The surface gradient method leads to Dh = 0 at the cell faces in the case of stagnant water conditions. The velocities at cell interfaces are also calculated from the depth integrated velocities uh and vh at cell faces as ðuhÞL : hL

ð3:18Þ

ð3:12Þ 3.3. Calculation of derivatives

~cnx 1

ð3:15Þ

Shu and Osher [26] showed that instead of using the ~ expressed in the RoeÕs approximate solver with DU eigenvector basis, the ENO schemes can be efficiently implemented using a simpler Rusanov-type approximate solver with artificial diffusive flux

uL ¼

with the corresponding eigenvectors 0 1 0 1 1 0 1 2 B C B C u þ ~cnx A; ~ ~ e~ ¼ @ ~ e~ ¼ @ ~cny A; 0

 1 ½DðhvÞ  ~vDhnx  ½DðhuÞ  ~uDhny : ~c

ð3:9Þ

k¼1

0

~a2 ¼

 Dh 1  DðhuÞnx þ DðhvÞny  ð~unx þ ~vny ÞDh ; 2 2~c ð3:14Þ

3.2. The surface gradient method

k

~ ak j~ e~ ; ak j~

~ is the Jacobian matrix of the flux where A ¼ oð~ F ~ nÞ=oU e is RoeÕs average projection in the normal direction, A e U ~ ~ with Jacobian matrix, which satisfies DF ¼ AD ~ nÞ e ¼ oðF  ~ A u 0 o~

~a1;3 ¼

ð3:13Þ

3 B C u  ~cnx A; ~ e~ ¼ @ ~

~v  ~cny respectively. The ~ak coefficients in (3.9) arise from the ~ in the basis ð~ decomposition of DU e1 ;~ e2 ;~ e3 Þ, and they depend on the jumps D = ( )R  ( )L with

The derivatives of a scalar variable c on unstructured grids are calculated using the divergence theorem, and we obtain   Z oc 1 oc c1 Dy 1 þ c2 Dy 2 þ c3 Dy 3 dA  ¼ ; ox i Ai A ox Ai ð3:19Þ   Z oc 1 oc c1 Dx1 þ c2 Dx2 þ c3 Dx3 dA   ¼ ; oy i Ai A oy Ai ð3:20Þ where (see Fig. 2)

A. Mohamadian et al. / Advances in Water Resources 28 (2005) 523–539

L Fig. 2. The triangular cell i with the quantities cR j and cj computed at the right (R) and left (L) sides respectively at a given face j (j = 1, 2, 3). The coordinates (xj, yj), j = 1, 2, 3, are located at the three vertices of cell i.

Dy 1 ¼ y 3  y 2 ;

Dx1 ¼ x3  x2 ;

c1 ¼ ðcL1 þ cR 1 Þ=2; ð3:21Þ

Dy 2 ¼ y 1  y 3 ;

Dx2 ¼ x1  x3 ;

c2 ¼ ðcL2 þ cR 2 Þ=2; ð3:22Þ

527

is guaranteed. On structured grids, as was shown by Nujic [24], this is easily possible, but on unstructured grids a compatible discretization of the BS term is not as straightforward as on structured grids. Nujic used the Shu and Osher [26] method (SO) with k ’ 0.4 and mentioned that the SO scheme is attractive for solving dam break problems. However, as shown in the following, the SO scheme produces significant numerical diffusion in the case of recirculating flows. In order to reduce it, we use here the RoeÕs approximate solver (3.9) instead of the SO scheme, employed in the NujicÕs method, as follows. The discretized SW equations using the RoeÕs method lead to  n ~ nþ1 ¼ U ~ n  0:5 Dt ð~ FL þ~ U F R  D~ F Þiþ1=2 Dx  n  ð~ FL þ~ F R  D~ F Þi1=2 þ Dt~ S; ð3:25Þ or, by rearranging the terms ~ nþ1 ¼ U ~n þ ~ ~ U L þ K;

ð3:26Þ

where

~¼ K



k1 k2

 ¼

Dy 3 ¼ y 2  y 1 ;

! 0  ; Dt Dt hi ðziþ1=2  zi1=2 Þ  0:5g Dx ðh2L þ h2R Þiþ1=2  ðh2L þ h2R Þi1=2 g Dx

Dx3 ¼ x2  x1 ;

c3 ¼ ðcL3 þ cR 3 Þ=2: ð3:23Þ

This approach may be used for calculating the gh(oz/ ox) and gh(oz/oy) terms and also the viscous terms, i.e.   ou u1 Dy 1 þ u2 Dy 2 þ u3 Dy 3 mh  mhi : ð3:24Þ ox i Ai The treatment of variable topographies now follows. 3.4. Variable topographies As mentioned in the previous sections, the bed slope (BS) term gh oz/ox is included in the source term, and it is usually discretized using a centered scheme. This leads to an imbalance between the depth gradient (DG) term 0.5g oh2/ox and the BS term and hence to an artificial source term in the numerical solution. In order to overcome this problem, the NujicÕs method [24] is adopted and modified here in the following manner. In the NujicÕs method, the DG term is extracted from the flux functions and it is discretized using a central scheme. However, the eigenvalues of the original system are used to calculate the fluxes using the SO scheme. The BS term is then dicretized in the NujicÕs method such that the compatibility between the DG and BS terms

ð3:27Þ

and ~ L includes all the remaining terms of the right hand side of (3.25). In the case of stagnant conditions g is constant, and Dh is thus zero at the cell faces using the surface gradient method. Because the stagnant conditions, we also have uh = 0 and D(uh) = 0. Therefore, the ~ ak coefficients calculated using (3.14) and (3.15) vanish and hence, D~ F ¼ 0. All other terms in ~ L also vanish in the case of stagnant conditions. In order to satisfy ~ must also vanish, which is not the case the C property K when using k2 as defined in (3.27). By rewriting k2 in (3.27) as Dt Dt hi ðziþ1=2  zi1=2 Þ  0:5g Dx Dx h    ðhL Þiþ1=2 þ ðhL Þi1=2 ðhL Þiþ1=2  ðhL Þi1=2   i þ ðhR Þiþ1=2 þ ðhR Þi1=2 ðhR Þiþ1=2  ðhR Þi1=2 ,

k2 ¼ g

ð3:28Þ using the following approximations   hi ¼ 0:5 ðhR Þiþ1=2 þ ðhR Þi1=2 þ OðDx2 Þ;

ð3:29Þ

  hi ¼ 0:5 ðhL Þiþ1=2 þ ðhL Þi1=2 þ OðDx2 Þ;

ð3:30Þ

528

A. Mohamadian et al. / Advances in Water Resources 28 (2005) 523–539

and neglecting the second order terms, Eq. (3.28) leads to ðgL þ gR Þiþ1=2  ðgL þ gR Þi1=2 : ð3:31Þ k 2 ¼ gDthi 2Dx

The components k2 and k3 in (3.36) and (3.37), are obviously zero in the case of stagnant conditions and hence, the C property is satisfied. Note (3.36) and (3.37) reduce to (3.31) in the 1-D case.

The k2 component in (3.31) obviously vanishes in the case of stagnant conditions and hence the C property is now satisfied. In the tests presented below, it is shown that neglecting the second order terms in (3.29) and (3.30) has a minor effect on the TVD property of the RoeÕs method, while considerably improves its performance in the case of variable topographies. This approach may be viewed as a mass-conservative approach. Although the discretization of the gradient term is not conservative, the convection terms are treated in a conservative manner, hence the mass is conserved, which is crucial for the long term simulations. A similar approach is used for the 2-D case. The bed slope term in a triangular control volume i may be discretized using the central scheme (3.19) as       oz oz z1 Dy 1 þ z2 Dy 2 þ z3 Dy 3 gh  ghi  ghi ; ox i ox i Ai

3.5. Entropy and high resolution corrections

ð3:32Þ ~ ¼ ðk 1 ; k 2 ; k 3 Þ in (3.26), with k1 = 0 and and K Dt h k2 ¼ g hi ðz1 Dy 1 þ z2 Dy 2 þ z3 Dy 3 Þ Ai      0:5 h2L þ h2R 1 Dy 1 þ h2L þ h2R 2 Dy 2 i   þ h2L þ h2R 3 Dy 3 , Dt h hi ðz1 Dx1 þ z2 Dx2 þ z3 Dx3 Þ Ai      0:5 h2L þ h2R 1 Dx1 þ h2L þ h2R 2 Dx2 i   þ h2L þ h2R 3 Dx3 :

The RoeÕs method violates the entropy condition in the case of sonic rarefaction (transcritical flow), and produces a shock inside a sonic rarefaction, which is obviously incorrect for the SW equations. Many methods have been proposed to overcome this problem. Here, we use the approach proposed by Harten and Hyman [12] with D~ F ¼

3 X

  k ~ak w ~ak ~ e~ ;

ð3:38Þ

k¼1

instead of (3.9). Here, the entropy correction function w is defined as  jwj if jwj P e; wðwÞ ¼ ð3:39Þ e if jwj < e: Harten and Hyman [12] proposed to calculate e as   ð3:40Þ e ¼ max 0; ~ak  akL ; akR  ~ak : k

ð3:33Þ

k 3 ¼ g

ð3:34Þ

Then, Dy3 is replaced in (3.33) by Dy3 = (Dy1 + Dy2) and using ðhR Þ1 þ ðhR Þ3 ðhL Þ1 þ ðhL Þ3  2 2 ðhR Þ2 þ ðhR Þ3 ðhL Þ2 þ ðhL Þ3  ; ð3:35Þ  2 2 we obtain    Dt k 2 ¼ 0:5g hi ðgL þ gR Þ1 Dy 1 þ ðgL þ gR Þ2 Dy 2 Ai  ð3:36Þ þðgL þ gR Þ3 Dy 3 : hi 

Similarly, the same procedure leads to    Dt k 3 ¼ 0:5g hi ðgL þ gR Þ1 Dx1 þ ðgL þ gR Þ2 Dx2 Ai  ð3:37Þ þðgL þ gR Þ3 Dx3 :

The first order RoeÕs scheme can be extended to a high resolution method by adding a correction term (based on the Lax–Wendroff scheme) to D~ F in the first order flux (3.5). The resulting D~ F is written in the following form:    3  X k Dt 1  1  k j~ak j uk ~ak wð~ak Þ~ e~ ; ð3:41Þ D~ F ¼ Dx k¼1 where uk is the flux limiter associated with the kth com ponent of the decomposition. Here, D~ F vanishes in the case of stagnant conditions (since ~ak vanishes) and hence, the C property still holds with the high resolution scheme. The flux limiter approach (3.41) has also been employed by Hubbard and Navarro [13] on a triangular grid with an upwind discretization of the source terms and led to good results. The implementation of (3.41) for 2-D cases has been described in LeVeque [19]. All the computations presented below are performed by using the minmod flux limiter as described in LeVeque [19]. 3.6. Stability and boundary conditions Due to the explicit discretization of the convective terms, a CFL number less than 1 is an upper limit for stability. This is the case when the non-linear effects are not dominant, like in the LeVequeÕs 1-D test case presented below using a CFL value of 0.99. However, for most existing schemes, the stability considerations

A. Mohamadian et al. / Advances in Water Resources 28 (2005) 523–539

are more restrictive when nonlinear effects are dominant, and this is also the case for the present scheme. Based on numerical experiments, there is no significant difference between the stability condition of the original RoeÕs scheme and the modified one proposed here, in the case of smooth flows. On the other hand, the proposed corrections in (3.29) and (3.30) induce a minor effect on the TVD property of the RoeÕs scheme in simulating discontinuous flows with a high CFL number. However, by choosing a slightly smaller time step (usually 10–20% less than the original RoeÕs scheme), the oscillations disappear. For example, in the dambreak test case 4.1 presented below, using a first order scheme, the oscillations are not seen when the CFL number is less than 0.83, and the scheme is stable up to a CFL of 0.93. When using high-resolution scheme, the oscillations are not observed when the CFL number is less than 0.75 and the scheme is stable up to a CFL of 0.9. In 2-D, the CFL condition is defined as " pffiffiffiffiffiffiffiffiffiffiffiffiffiffi # ð u2i þ v2i þ ci ÞDt CFL ¼ min ; ð3:42Þ d ij where dij is the distance between the centroid of each cell i and its neighbors j. The numerical treatment of the boundary conditions is performed by setting the variables at the boundary edges based on the theory of characteristics [4]. For subcritical flows, two external conditions must be specified at inflow boundaries, whereas only one is required at the outflow one. However, two-dimensional supercritical flows require the imposition of three inflow boundary conditions and none at the downstream side. In the particular case of solid wall, both velocity components are set equal to zero. The depth variables are calculated by using the information carried out by the outgoing bicharachteristics (the Riemann invariants). However, in most practical cases those may be simply set to the corresponding values of the adjacent inner cells.

4. Numerical results In order to study the performance of the numerical scheme, some tests have been performed herein. Tests 4.1 and 4.2 show the ability of the model for simulating dam-break type flows. The remaining tests show the performance of the proposed scheme for variable topographies and recirculating flows. In tests 4.1–4.4 and 4.9.1 we examine 1-D problems and the other tests are performed using the two dimensional model. In all tests the CFL number ranges between 0.35 and 1.0. In all figures and tables, the SI system has been used, i.e. the water depth and the water surface elevation are in

529

m and the water discharge is in m3/s. In all tests a high resolution scheme is employed except otherwise mentioned. 4.1. The dam break problem The dam break problem is a common test for evaluating shock capturing schemes in shallow flows. Here, a dam is located at the mid point of a channel of length 50 m. The water depth at left and right hand side of the dam is 0.5 m and 0.25 m respectively. The dam is instantaneously removed across its entire width and the flow conditions are computed up to time t = 8 s. A CFL number of 0.65 was chosen in order to compare the first and second order schemes introduced in (3.38) and (3.41), respectively. The results consist of a right sock and a left rarefaction. The results for the depth (h) and discharge (uh) are presented in Fig. 3a and b respectively. As it can be seen, both first and second order schemes are able to provide good results without numerical oscillations. Further, the second order scheme requires less grid points than the first order one, to capture the shock wave. In the dam break tests, the maximum value of hL/hR is a measure of the performance of the numerical scheme. In the next section, we show that the proposed method can simulate the dambreak problem with all depth ratios (even on a dry bed). 4.2. ToroÕs Riemann problems We now investigate the ability of the proposed scheme in performing some challenging test cases designed by Toro [27]. A summary of the initial conditions used in these test cases is presented in Table 1, where hL, uL, hR, uR, denote the initial depth and velocity in the left and right hand side of the discontinuity, x0 is the position of the discontinuity and tout is the output time. In all tests a wide horizontal, rectangular and frictionless channel of length 50 m is used. In all tests a second order scheme is used and the results are compared with those of the Roe scheme (second order) and the analytical solutions. 4.2.1. Left sonic rarefaction and right shock The initial condition has been chosen as to produce a strong shock wave and a sonic or transcritical left rarefaction wave. The results of this test are presented in Fig. 4, and they are compared with those of the original RoeÕs scheme. As it can be seen, both the RoeÕs scheme and the proposed method give good results. 4.2.2. Two rarefaction and nearly dry bed In this test, the initial conditions have been chosen in order to produce two rarefaction waves and a nearly dry bed. The results are shown in Fig. 5 for the original RoeÕs scheme and the proposed method. They confirm

530

A. Mohamadian et al. / Advances in Water Resources 28 (2005) 523–539

0.5

0.5

Exact Calculated

Exact Calculated

0.45

0.4

0.4

h

h

0.45

0.35

0.35

0.3

0.3

0.25

0.25 0

(a)

10

20

30

40

0

50

10

20

30

40

50

x

x

0.25

0.25 Exact Calculated

Exact Calculated

0.15

0.15 q

0.2

q

0.2

0.1

0.1

0.05

0.05

0

0

0

(b)

10

20

30

40

50

0

10

20

x

30

40

50

x

Fig. 3. Calculated depth (a) and discharge (b) at t = 1.6, 3.2, 4.8, 6.4 and 8 s with the proposed method; second order (left) and first order (right).

Table 1 Data for ToroÕs test problems Test

hL (m)

uL (m)

hR (m)

uR (m)

x0 (m)

tout (s)

4.2.1 4.2.2 4.2.3 4.2.4

1.0 1.0 1.0 0.1

2.5 5.0 0.0 3.0

0.1 1.0 0.0 0.1

0.0 5.0 0.0 3.0

10.0 25.0 20.0 25.0

7.0 2.5 4.0 5.0

that the proposed modifications do not change the performance of the RoeÕs method for this problem.

4.2.3. Dambreak on dry bed Jha et al. [15] reported that the RoeÕs method using the second order accurate correction fails to simulate the dambreak problem when hL/hR < 0.002. In order to overcome this problem, the flux limiter function is also set here to zero when h < 105, i.e. the model is second order accurate everywhere except close to the shock region and the wet/dry border. The results are presented in Fig. 6, and they are again compared with those of the original RoeÕs scheme (with the modified flux limiter function). They show that the proposed method is able

1.1

1.1

1

1 Exact Calculated

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1 0

10

20

30 x

40

Exact Calculated

0.9

h

h

0.9

50

0

10

30

20

40

50

x

Fig. 4. Left sonic rarefaction and right shock test case. Calculated water depth with the proposed method (left) and the RoeÕs scheme (right).

A. Mohamadian et al. / Advances in Water Resources 28 (2005) 523–539

1.1

1

1

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

h

h

1.1

531

0.5

0.5

0.4

0.4

0.3

0.3

Exact Calculated

0.2

Exact Calculated

0.2

0.1

0.1

0 0

10

20

30

40

0 0

50

10

20

x

30

40

50

x

Fig. 5. Two rarefactions and nearly dry bed test case. Calculated water depth with the proposed method (left) and the RoeÕs scheme (right).

0.8

0.8

0.6

0.6 h

1

h

1

0.4

0.4 Exact Calculated

0.2 0

Exact Calculated

0.2 0

0

10

20

30

40

50

0

10

x

20

30

40

50

x

Fig. 6. Dambreak on dry bed. Calculated water depth with the proposed method (left) and the RoeÕs scheme (right).

to simulate the dambreak type flows whatever the choice of depth ratio hL/hR.

where L = 14,000 m is the channel length. By defining H(x) = z(x), the initial conditions write

4.2.4. Generation of a dry bed The initial condition has been chosen in order to produce a solution consisting of two rarefaction waves with a portion of a dry bed between them. The results of this test are presented in Fig. 7, and they are again compared with those of the RoeÕs scheme. As it can bee seen, both schemes produce similar results, i.e. the proposed modification does not change the quality of the results of the RoeÕs scheme for this test. Hence, tests 4.2.3 and 4.2.4 show the ability of the proposed scheme in simulating dry bed problems.

hðx; 0Þ ¼ H ðxÞ;

ð4:2Þ

uðx; 0Þ ¼ 0

ð4:3Þ

4.3. Tidal wave flow 4.3.1. Tidal wave flow over regular topography This one dimensional test was proposed by Bermudez and Vazquez [2]. The topography is defined as    40x 4x 1 zðxÞ ¼ 50:5 þ  10 sin p  ; ð4:1Þ L L 2

with the boundary conditions    4t 1 þ hð0; tÞ ¼ H ð0Þ þ 4  4 sin p ; 86; 400 2 uðL; tÞ ¼ 0:

ð4:4Þ ð4:5Þ

For this test, Bermudez and Vazquez [2] derived the following analytical solution:    4t 1 hðx; tÞ ¼ H ðxÞ þ 4  4 sin p þ ; ð4:6Þ 86; 400 2    ðx  LÞp 4t 1 cos p þ uðx; tÞ ¼ : ð4:7Þ 5400hðx; tÞ 86; 400 2 At time t = 7552.13, the results are shown in Fig. 8, and are in a good agreement with the exact solution.

A. Mohamadian et al. / Advances in Water Resources 28 (2005) 523–539

0.1

0.1

0.075

0.075

0.05

0.05

h

h

532

0.025

0.025 Exact Calculated

Exact Calculated

0

0 0

10

20

30

40

50

x

10

0

20

30

40

50

x

Fig. 7. Generation of a dry bed. Calculated water depth with the proposed method (left) and the RoeÕs scheme (right).

0.2

0.06 Exact Calculated

0.16

Exact Calculated

0.04

u

η

0.12 0.08

0.02 0.04 0 0

5000

0

10000

x

0

500

1000

1500

x

Fig. 8. Tidal wave flow in a channel with variable topography. Fig. 9. Tidal wave flow over an irregular bed.

4.3.2. Tidal wave flow over an irregular bed Here, the topography is defined in Table 2. This irregular bed was proposed at a workshop on dambreak wave simulations [11]. The initial and boundary conditions are identical to those used in test 4.3.1, with

( zðxÞ ¼

0:2  0:05ðx  10Þ2

if 8 < x < 12;

0

otherwise:

ð4:8Þ

4.4. Steady flow over a bump

Depending on the initial and boundary conditions, the flow may be subcritical, transcritical (with or without a steady shock), or supercritical. The analytical solution of this problem is given in [11]. The global relative error R is defined as vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u uX hn  hn1 2 i i ; ð4:9Þ R¼t hni i

This problem was also proposed at the workshop on dam-break simulations [11]. The topography is defined as

where hni and hn1 are the local water depths at the curi rent and previous time levels. A mesh of 150 nodes was used in all computations.

L ¼ 1500 m;

H ðxÞ ¼ H ð0Þ  zðxÞ;

H ð0Þ ¼ 16 m:

The analytical solution is given by (4.6) and (4.7). At time t = 10,800 s the results, shown in Fig. 9, indicate the ability of the scheme in simulating irregular topographies.

Table 2 Topography for irregular bed X (m) Z (m)

0 0

50 20

100 2.5

150 5

250 5

300 3

350 5

400 5

425 7.5

435 8

450 9

475 9

500 9.1

505 9

X (m) Z (m)

530 9

550 6

565 5.5

575 5.5

600 5

650 4

700 3

750 3

800 2.3

820 2

900 1.2

950 0.4

1000 0

1500 0

A. Mohamadian et al. / Advances in Water Resources 28 (2005) 523–539 1.58

1.1 1

1.56

Exact Calculated

0.9

Exact Calculated

1.54

0.8

1.52

q

η

533

0.7 0.6

1.5 1.48

0.5

1.46

0.4

1.44 0.3 0

10

x

0

20

10

20

x

Fig. 10. Steady transcritical flow over a bump without a shock: water surface elevation (left) and discharge (right).

Three different cases are considered: (i) Transcritical flow without a shock The boundary conditions are defined as Downstream: The water level h = 0.66 m is imposed only when the flow is subcritical. Upstream: The discharge uh = 1.53 m3/h is imposed. The surface profile is plotted in Fig. 10, and it shows a good agreement with the analytical solution. The computed discharge is also compared with the theoretical one in Fig. 10. The numerical oscillations observed in Fig. 10 (right) are also present in most existing schemes. This behavior is related to the fact that the C property, although it is widely accepted in the literature as a good-enough measure to test the ability of the numerical schemes in the presence of

real topographies, does not guarantee that the steady state condition with non-zero discharge is well captured. (ii) Transcritical flow with a shock The boundary conditions are defined as Downstream: The water level h = 0.33 m is imposed. Upstream: The discharge uh = 0.18 m3/h is imposed. Fig. 11 (left) shows that a good agreement is obtained between the analytical solution and the computed water surface elevation. A comparison of the computed discharge with the theoretical results is shown in Fig. 11 (right). The oscillations have the same origin than previously mentioned in (i). (iii) Subcritical flow The boundary conditions are now defined as

0.5

0.26

0.45

0.24

Exact Calculated

0.4

Exact Calculated

0.22 0.2 0.18

0.3

q

η

0.35

0.16 0.14

0.25

0.12

0.2

0.1

0.15

0.08

0.1 0

10

0.06

20

0

10

x

20

x

Fig. 11. As for Fig. 10 but for the steady transcritical flow over a bump with a shock.

4.48 2

4.46

Exact Calculated

1.98

Exact Calculated

4.44

4.4

q

η

4.42 1.96

4.38 1.94

4.36 4.34

1.92

4.32 4.3

1.9 0

10

20

x

0

10

20

x

Fig. 12. As for Fig. 10 but for the steady subcritical flow over a bump.

A. Mohamadian et al. / Advances in Water Resources 28 (2005) 523–539 1.00E+00

1.E+00

1.E-01

1.00E-01

1.E-01

1.E-02 1.E-03 1.E-04 1.E-05 1.E-06

1.00E-02 1.00E-03 1.00E-04 1.00E-05

1.E-07

1.00E-06

1.E-08

1.00E-07 0

500

1000

1500

2000

2500

3000

3500

4000

Releative error R

1.E+00

Releative error R

Releative error R

534

1.E-02 1.E-03 1.E-04 1.E-05 1.E-06 1.E-07

0

2000

4000

6000

8000

10000

12000

0

1000

2000

Iteration number

Iteration number

3000

4000

5000

6000

7000

8000

Iteration number

Fig. 13. Convergence history for tests (i), (ii), (iii), from left to right respectively.

Fig. 14. Unstructured triangular grid used for the steady flow over a bump.

Fig. 15. Different cases of flow over a bump (the 2-D model): Water surface elevation.

Downstream: The water level h = 2 m is imposed. Upstream: The discharge uh = 4.42 m3/h is imposed. The numerical results are depicted in Fig. 12 (left), and again they show an excellent agreement between the analytical solution and the computed water surface elevation. A comparison of the computed discharge with the analytical results is shown in Fig. 12 (right). Again, the oscillations have the same origin than previously mentioned in (i). The convergence history performed for the tests (i), (ii), (iii) is shown in Fig. 13.

In order to examine the 2-D model in different flow regimes over variable topographies, the present test case is also performed with the 2-D model. The unstructured grid used here is shown in Fig. 14. The results corresponding to the boundary conditions (i), (ii), (iii) are compared with the analytical solution in Fig. 15, and they show that the current method can simulate all cases accurately. 4.5. Bore reflection 4.5.1. Alcrudo and NavarroÕs oblique wave problem An oblique wave is generated by the interaction between a supercritical flow and a converging wall, with

Fig. 16. Oblique wave problem: (a) 3-D view of calculated water level, and (b) grid and contours of water level.

A. Mohamadian et al. / Advances in Water Resources 28 (2005) 523–539

535

Fig. 17. (a) A schema of the shock-on-shock interaction problem (left) and the triangular grid used in the simulation of the shock-on-shock interaction problem (right). (b) Shock on shock interaction: contours of calculated water surfaces elevation.

an angle of h = 8.95 between the wall and the streamwise direction. The upstream condition is a uniform supercritical flow with a depth of 1 m, and the flow velocity is such that u = 8.57 m/s and v = 0.0 m/s. The analytical ffiffiffiffiffiffiffiffiffiffiffiffiffiffi solution for this problem is hd = 1.5 m and p u2 þ v2 ¼ 7:9525 m/s, respectively, and the numerical results are shown in Fig. 16a and b. The calculated depth ffiffiffiffiffiffiffiffiffiffiffiffiffiffiand flow velocity are hd = 1.4995 m and p u2 þ v2 ¼ 7:9556 m/s. Hence, it appears that the model is able to simulate this test case accurately, without exhibiting numerical oscillations.

The water depth is 1 m in the undisturbed region, the upstream Froude number is 2.7 and the wall deflection angle is 12. The simulation was pursued until a steady state solution was obtained using the high-resolution scheme. A two dimensional view of the water level is shown in Fig. 17b. A comparison between the computed and the exact solution is presented in Table 3. As it can be seen, the numerical model can accurately simulate the water surface and the shock-on-shock interaction.

4.5.2. Shock-on-shock interaction This test, including an exact solution, was proposed by Causon et al. [6]. A schematic representing the problem and the unstructured grid used in the simulation, are presented in Fig. 17a. Based on the interaction of the shock-waves, four regions (1), (2), (3) and (4) are formed, where bA is the angle between the shock wave and the original flow direction and bD represents the angle between the shock wave and the side wall. The numerical simulation has been performed by using an averaged value of 0.65 for the CFL number.

To test again the performance of the proposed method, we now consider the breaking of a cylindrical dam and the time evolution of the subsequent waves [4]. Two regions of still water are separated by a cylindrical wall of radius r = 11 m, such that the water depth is 10 m on the inner side and is 1 m outside the dam. The computed water surface profile using 6328 cells at time t = 0.69 s is shown in Fig. 18 (left) and it is compared with the HLL scheme results (Fig. 18, right). As it can be seen, the method performs well and smooth results are obtained. Further, the radial symmetry of the computed solution numerical method is respected.

Table 3 Comparison between calculated results and exact solution

4.7. Recirculating flow after a sudden expansion

State 2

The 2-D laminar flow past a sudden expansion in a side wall (horizontally analogous to the backward facing step) provides a well documented benchmark test for evaluating the level of numerical diffusion in conservation law schemes. Many researchers [7], have provided experimental and numerical data on the flow pattern immediately after a backward facing step, where a steady vortex appears with its length depending on the

Variable Fr Depth (m) Velocity (m/s) bA ()

State 4 Theory Computed Variable

Theory Computed

1.868 1.676

1.87 1.674

Fr Depth (m)

1.249 2.562

1.249 2.560

7.571

7.58

Velocity (m/s)

6.258

6.26

33.688

33.69

48.102

48.764

bD ()

4.6. Circular dam break

536

A. Mohamadian et al. / Advances in Water Resources 28 (2005) 523–539

Fig. 18. Computed water surface level at t = 0.69 s for the circular dam break problem; a 3-D view (top) and a comparison with the HLL scheme (bottom).

Fig. 19. Predicted velocities with the proposed method in the xdirection compared with the experimental data.

Reynolds number, defined as Re = bU1/m, where U1 is the mean inlet flow velocity, b is the dimension of the wall expansion and m is the fluid kinematic viscosity coefficient. In Denham and Patrick [7], U1 = qx/D where qx is the inlet x-direction inlet momentum flux and D is the flow depth. The channel width is 2b before, and 3b after the expansion, where b = 1 m. At the expansion, the mean flow velocity is U1 = 0.5 m/s, the flow depth at the downstream is 1 m and the eddy viscosity is

0.00685 m2/s, also considered by Denham and Patrick, corresponding to Re = bU1/m = 73. No bottom friction is applied. All numerical experiments are performed using the first order scheme, and an structured 300 · 100 grid is employed. Fig. 19 presents the predicted velocity vectors in the x direction and the experimental data. The predicted steady-state velocity profiles at different positions along the channel past the expansion are in very close agreement with Denham and PatrickÕs experimental results, as shown in Fig. 19. The streamtraces predicted by using the proposed scheme are shown in Fig. 20. The length of the recirculation region is approximately 3.89b which is near the value 3.95b given by Denham and Patrick [7] from experimental data. This recirculation zone is correctly modeled, thus demonstrating the ability of the proposed scheme to simulate circulating flows with a low level of numerical diffusion. The results of the NujicÕs method are also shown in Fig. 20. As it can be seen the length of the recirculating

Fig. 20. Streamtraces calculated with: (a) the proposed method, (b) the NujicÕs scheme, (c) the central upwind scheme, and (d) the HLL flux (all first order).

A. Mohamadian et al. / Advances in Water Resources 28 (2005) 523–539

zone is roughly half of the expected value, which shows a high level of numerical diffusion for the circulating flows. The HLL flux (used by [21]), and the KurganovÕs central upwind scheme [18] have even more numerical diffusion, and the recirculating zone is roughly 1.2b as shown in Fig. 20. 4.8. Partially closed channel A partially closed channel of length 3 km and width 1 km is considered under the steady state flow conditions. Although this is an academic case, a number of different models have been applied to it as reported in the literature [17,30,34]. A headland of length 0.4 km is located in the middle of the channel, as shown in Fig. 21. The bed roughness height is assumed to be 0.25 m and the horizontal mixing coefficient is 0.5 m2/ s. The upstream hydrodynamic boundary condition is assumed to be a discharge of 4000 m3/s and the downstream boundary is h = 6 m [17]. Fig. 22 shows the mesh and the streamlines predicted by the proposed model.

3000 m 400 m 1000 m

Q y Fig. 21. Schematic view of the channel.

Fig. 22. Streamline predicted by the proposed model.

537

The results of the depth-averaged velocities along the streamline B obtained with the present model are compared in Fig. 23 with three other models. Those are GEO-DIVAST model [17], which is a modified version of DIVAST [8] and uses a finite volume method; the finite element model SUTRENCH, developed by Delft Hydraulics [30]; and the ESMOR two dimensional depth averaged model, also developed by Delft Hydraulics [34]. Based on a comparison of the results in Fig. 23, the results of the proposed model are in a close agreement with other models. However, the proposed model has a wider applicability and is able to simulate super critical and transcritical flows on unstructured grids. Further, it should be mentioned that in this test case, the bed is movable and the current model successfully simulated the morphological processes around a groyne [22]. 4.9.1. One dimensional, small perturbation of a steady-state solution This test case is concerned with wave propagation and it has been proposed by LeVeque [20]. The channel has a length of 1.0 m and the topography is defined as 8 2 2 > < 0:8 expð5ðx  0:9Þ  50ðy  0:5Þ Þ zðx; yÞ ¼ for jx  0:5j 6 0:1; > : 0:0; otherwise: ð4:10Þ The flow velocity is zero and the surface profile is  1:0 þ e for 0:1 < x < 0:2; gðx; yÞ ¼ ð4:11Þ 1:0 otherwise: Two cases have been considered: e = 0.2 m and 0.01 m. The solutions obtained on a grid of 100 nodes and g = 1 m2/s are shown in Fig. 24. When e = 0.01 the experiment is more difficult to conduct. A CFL number of 0.85 and 0.99 is used for the cases of e = 0.2 m and 0.01 m, respectively. Contrary to most existing schemes, which have used much more grid points than here (600 in [13]), the present method can capture the quasi steady solutions quite well, with a very low level of numerical diffusion and without numerical oscillations. 4.9.2. Two dimensional, small perturbation of a steady-state solution In the test case presented by LeVeque [20], the SW equations are solved in the domain [0, 2] · [0, 1], and the bottom surface is an elliptical shaped hump   zðx; yÞ ¼ 0:8 exp 5ðx  0:9Þ2  50ðy  0:5Þ2 : ð4:12Þ

Fig. 23. Depth averaged velocities along the streamline B obtained by using different models.

The surface is initially flat with h(x, y) = 1 except for 0.05 < x < 0.15, where h(x, y) = 1.01 m. Fig. 25 displays the right-going disturbance as it propagates past the hump, using coarse (11,653 control volumes) and fine

538

A. Mohamadian et al. / Advances in Water Resources 28 (2005) 523–539

Fig. 24. Computed water surface level at t = 0.7 with e = 0.2 (left) and e = 0.01 (right).

5. Conclusion A mass conservative approach has been proposed to solve the shallow-water equations. It is able to accurately simulate mild flows such as recirculating and tidal flows with low numerical diffusion. Numerical results indicated that the proposed method can also accurately compute sub, super and transcritical flows with discontinuity over complicated topographies. In the proposed scheme, it is not necessary to perform an extra upwinding or Riemann solution for the source terms, which makes our method computationally efficient. Many test cases show that the proposed method can be efficiently used for a wide range of SW problems. Contrary to many available schemes, our approach can be easily implemented on unstructured grids and has the advantage of being flexible for irregular boundaries and local mesh refinement.

Acknowledgments A.M. is supported partially by the water research center of Iran, the HEMAT project. A.M. and D.Y.L. are supported by grants from the Natural Science and Engineering Council (NSERC) and the Fonds Que´be´cois de la Recherche sur la Nature et les Technologies (FQRNT).

References Fig. 25. Computed water surface level with a coarse grid (left) and a fine grid (right) for t = 0.6, 0.9, 1.2, 1.5 and 1.8 s.

(29,877 control volumes) unstructured grids. The current method is clearly able to simulate this problem without facing spurious oscillations due to variable bed topography.

[2] Bermudez A, Vazquez ME. Upwind methods for hyperbolic conservation laws with source terms. Comput Fluid 1994;23:1049. [3] Alcrudo F, Benkhaldoun F. Exact solutions to the Riemann problem of the shallow-water equations with a bottom step. Comput Fluid 2001;30:643. [4] Alcrudo F, Garcia-Navarro P. A high resolution Godunov-type scheme in finite volumes for the 2D shallow-water equations. Int J Numer Meth Fluid 1993;16:489.

A. Mohamadian et al. / Advances in Water Resources 28 (2005) 523–539 [6] Causon DM, Mingham CG, Ingram DM. Advances in calculation methods for supercritical flow in spillway channels. J Hydraul Eng ASCE 1999;125:1039. [7] Denham MK, Patrick MA. Laminar flow over a downstreamfacing step in a two-dimensional flow channel. Trans Inst Chem Eng 1974;52:361. [8] Falconer RA. Numerical modeling of tidal circulation in harbors. J Hydraul Eng ASCE 1980;106(1):31. [9] Fenemma R, Chaudhry M. Simulation of one-dimensional dambreak flows. J Hydraul Res 1987;25(1):41. [10] Glaister P. Approximate riemann solutions of the shallow-water equations. J Hydraul Res 1988;26(3):293. [11] Goutal N, Maurel F, editors. Proceedings of the 2nd workshop on dam-break wave simulation. HE-43/97/016/B, France; 1997. [12] Harten A, Hyman P. Self adjusting grid methods for onedimensional hyperbolic conservation laws. J Comput Phys 1983;50:235. [13] Hubbard ME, Garcia-Navarro P. Flux difference splitting and the balancing of source terms and flux gradients. J Comput Phys 2000;165:89. [14] Hubbard ME, Garcia-Navarro P. Balancing source terms and flux gradients in finite volume schemes. In: Toro EF, editor. Godunov methods: theory and applications. New York: Kluwer Academic Publishers; 2001. p. 447. [15] Jha AK, Akiyama J, Ura M. First and second-order flux difference splitting schemes for dam-break problem. J Hydraul Eng ASCE 1995;121(12):877. [16] Jin S. Steady-state capturing method for hyperbolic systems with geometrical source terms. Math Model Num Anal 2001;35: 631–46. [17] Kolahdoozan M. Numerical modeling of geo-morphological processes estuarine waters. PhD thesis, University of Bradford, UK, 1999. [18] Kurganov A, Levy D. Central-upwind schemes for the Saint– Venant system. Math Model Numer Anal 2002;36:397.

539

[19] LeVeque RJ. Finite volume methods for hyperbolic problems. UK: Cambridge University Press; 2002. [20] LeVeque RJ. Balancing source terms and flux gradients in highresolution Godunov methods. J Comput Phys 1998;146:346. [21] Mingham CG, Causon DM. High resolution finite-volume method for shallow-water equation flows. J Hydraul Eng ASCE 1998;124:605. [22] Mohammadian A, Azad FL, Tajrishi M. Two dimensional numerical simulation of flow and geo-morphological processes near headlands by using unstructured grid. In: XXX IAHR Congress. Greece; 2003. [24] Nujic M. Efficient implementation of non-oscillatory schemes of free surface flows. J Hydraul Res 1995;33(1):101. [26] Shu CW, Osher S. Efficient implementation of non-oscillatory shock-capturing schemes. J Comput Phys 1988;77:439. [27] Toro EF. Shock capturing methods for free surface shallow flows. John Wiley and Sons; 2000. [30] van Rijn LC. Principles of sediment transport in rivers, estuaries and coastal seas. The Netherlands: Aqua Publications; 1993. [31] Vazquez-Cendon ME. Improved treatment of source terms in upwind schemes for the shallow-water equations in channels with irregular geometry. J Comput Phys 1999;43(2):357. [32] Vukovic S, Sopta L. ENO and WENO schemes with the exact conservation property for one-dimensional shallow-water equations. J Comput Phys 2002;179:593. [34] Wang ZB. Mathematical modeling of morphological processes in estuaries. The Netherlands: Delft University of Technology; 1989. [35] Xu K. A well-balanced gas-kinetic scheme for the shallow-water equations with source terms. J Comput Phys 2002;178(2):533–62. [36] Zhao D, Shen HW, Tabios III GQ, Lai JS, Tan WY. Finitevolume two-dimensional unsteady-flow model for river basins. J Hydraul Eng ASCE 1994;120(7):863. [38] Zhou JG, Causon DM, Mingham CG, Ingram DM. The surface gradient method for the treatment of source terms in the shallowwater equations. J Comput Phys 2001;168:1.