Adaptive quadtree simulation of shallow flows with wet–dry fronts over complex topography

Computers & Fluids 38 (2009) 221–234

Contents lists available at 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

Adaptive quadtree simulation of shallow flows with wet–dry fronts over complex topography Qiuhua Liang a,*, Alistair G.L. Borthwick b a b

School of Civil Engineering and Geosciences, Newcastle University, Newcastle upon Tyne NE1 7RU, UK Department of Engineering Science, University of Oxford, Oxford OX1 3PJ, UK

a r t i c l e

i n f o

Article history: Received 13 February 2007 Received in revised form 18 June 2007 Accepted 16 February 2008 Available online 25 March 2008

a b s t r a c t Many natural terrains have complicated surface topography. The simulation of steep-fronted flows that occur after heavy rainfall flash floods or as inundation from dyke breaches is usually based on the nonlinear shallow water equations in hyperbolic conservation form. Particular challenges to numerical modellers are posed by the need to balance correctly the flux gradient and source terms in Godunov-type finite volume shock-capturing schemes and by the moving wet–dry boundary as the flood rises or falls. This paper presents a Godunov-type shallow flow solver on adaptive quadtree grids aimed at simulating flood flows as they travel over natural terrain. By choosing the stage and discharge as dependent variables in the hyperbolic non-linear shallow water equations, a new deviatoric formulation is derived that mathematically balances the flux gradient and source terms in cases where there are wet–dry fronts. The new formulation is more general in application than previous a priori approaches. Three benchmark tests are used to validate the solver, and include steady flow over a submerged hump, flow disturbances propagating over an elliptical-shaped hump, and free surface sloshing motions in a vessel with a parabolic bed. The model is also used to simulate the propagation of a flood due to a dam break over an initially dry floodplain containing three humps. Ó 2008 Elsevier Ltd. All rights reserved.

1. Introduction The non-linear shallow water equations are very useful for simulating long wave hydrodynamics when the vertical acceleration of water particles can be neglected and the flow can be reasonably assumed to be nearly horizontal. Given their hydrostatic nature, hyperbolic integral forms of the non-linear shallow water equations provide a surprisingly accurate representation of steepfronted flows, such as dam breaks, dyke breaches, flash floods, tidal bores, flood waves, and coastal inundation. Godunov-type finite volume solvers of the non-linear shallow water equations have a shock-capturing property that is essential to preserve the discontinuous or steeply varying gradients that occur in transcritical and sharp-fronted shallow flows. Typical examples of Godunovtype shallow flow models are given by Alcrudo and García-Navarro [1], Zhao et al. [2], Fraccarollo and Toro [3], Anastasiou and Chan [4], Sleigh et al. [5], Causon et al. [6], Rogers et al. [7], Toro [8], and Liang et al. [9]. Although numerical solvers of the shallow water equations have been in existence for more than 20 years, two persistent problems dog the schemes in practice. One is the treatment of source terms in the equation system built from conservation laws. * Corresponding author. Tel.: +44 191 2226413; fax: +44 191 2226502. E-mail address: [email protected] (Q. Liang). 0045-7930/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2008.02.008

The other is connected with the modelling of the moving wet–dry interface at a shallow water wave front. With regard to the first problem, considerable progress has been made in devising schemes for computing the source terms. Several numerical and mathematical treatments are proposed in the literature for balancing the flux gradient and source terms in Godunov finite volume schemes (see e.g. [10–16]). Of these, the mathematical conditioning method of Rogers et al. [16] leads to the shallow water equations in deviatoric form by subtracting an equilibrium solution. The approach is attractive in that it is generalised, and leads to a balanced set of hyperbolic equations that does not require specific numerical algorithms to have to be implemented when computing the source term. However, Rogers et al.’s formulation of the shallow water equations is expressed in terms of still water depth and surface elevation above the still water level, and so is difficult to apply to cases involving wetting and drying as the still water depth is undefined in dry areas. This brings us to the second problem: the modelling of a moving wet–dry interface when solving the hyperbolic shallow water equations. In Godunov-type finite volume schemes where the depth-averaged velocity components are normally determined by dividing the discharge per unit width by the local water depth, erroneously high velocities can be predicted at wet–dry fronts where the water depth is extremely small. These can lead in turn to physically meaningless predictions of negative water depth, and instability. Although a variety of wetting

222

Q. Liang, A.G.L. Borthwick / Computers & Fluids 38 (2009) 221–234

and drying algorithms have been proposed specifically for the hyperbolic shallow water equations (e.g. [17,18]), most of these methods are not general and further research is required. A further important issue is that of computational efficiency, especially when it is intended to use solvers based on the hyperbolic non-linear shallow water equations to simulate multiple-scale shallow flows over complicated bed topography, such as is the case with flood inundation. To reduce the computational overhead, simplified diffusion-wave models have been derived from the fully twodimensional non-linear shallow water equations by neglecting the dynamic convection terms (e.g. [19–21]). However, when a flood event involves complex flow situations (e.g. transcritical flows, sudden breaching or overtopping of a dam or dyke, or inundation over very complicated terrain), it is necessary to include the dynamic terms in the governing equations in order to evaluate accurately the front velocity and flood route. Efficient solvers of the hyperbolic non-linear shallow water equations are therefore required for complex large-scale flood simulations. Various approaches to adaptive mesh resolution for improving computational efficiency have been suggested for the shallow water equations in the literature (see e.g. [5,7,9,22–26]). Typical mesh generation techniques may be classified as block-structured (uniform, curvilinear, mapped), unstructured (including advancing front and Delauney–Voronoï methods), hierarchical, and hybrid. Extensive reviews of grid generation techniques are available in the literature; two of the most noteworthy are by George [27] and Thompson et al. [28]. An example of an adaptive shallow flow model based on boundary-fitted curvilinear grids is that described by Ivanenko and Muratova [22], who used the theory of harmonic maps to create grids that conformed to the bed topography of the Azov Sea. In Ivenenko and Muratova’s model the grid cells could change size according to local flow features without altering the total number of cells. An obvious advantage of adaptive boundary-fitting is the accurate description of curved shorelines. However, there is the drawback that the highly stretched curvilinear cells created by the adaptation process may adversely affect solution accuracy and stability. Examples of adaptive shallow flow models based on unstructured triangular grids are given by Sleigh et al. [5] and Skoula et al. [24]. Unstructured grids are widely utilised in computational fluid dynamics because of their robustness and ease of generation, especially for domains with highly complicated boundary configurations. A shortcoming is the awkward way in which grid connectivity is handled, which can slow down the speed of the solver when applied on an adaptive unstructured mesh. Hierarchical (quadtree, tritree) grids are created by domain decomposition, and although they appear unstructured when fitted to a complicated domain, their underlying tree structure is easy to interrogate in order to identify neighbouring cells. Typical adaptive quadtree solvers of the shallow water equations are presented by Rogers et al. [7,16], Liang et al. [9], and Krámer and Józsa [26]. Applications include steep-fronted dam breaks, dyke breaks, and wind-induced circulation in shallow lakes. Lamby et al. [29] use fully adaptive multiscale wavelet methods to simulate transcritical shallow flows, with grids that are effectively hierarchical. Recently, Liang et al. [25] demonstrated that a Godunov-type solver of the non-linear shallow water equations on adaptive quadtree grids can produce simulations of similar resolution but at nineteen times less cost than for the same solver on a uniform grid. It is therefore reasonable to anticipate that an adaptive grid based full shallow flow solver could be as efficient as a simplified diffusion-wave model. The paper presents details of a computationally efficient numerical approach to simulating complex shallow flows, which is potentially applicable to large-scale flood inundation over natural terrain. In Sections 2 and 3, a new deviatoric formulation of the hyperbolic non-linear shallow water equations is derived that is more generally applicable than previous versions in that it encap-

sulates wetting and drying processes. Section 4 briefly introduces the adaptive quadtree grid. Section 5 describes the numerical solution of the equations on dynamically adaptive quadtree grids using a second-order Godunov-type finite volume scheme. Section 6 presents the results of model validation tests and a simulation of a dam break flood wave as it travels over complicated topography. 2. Governing equations The two-dimensional non-linear shallow water equations may be derived by depth-integrating the three-dimensional Reynoldsaveraged Navier–Stokes equations, neglecting the vertical acceleration of water particles, and taking the pressure distribution to be hydrostatic [30]. In matrix form, a conservation law of the twodimensional non-linear shallow water equations may be written as ou of og þ þ ¼ s; ot ox oy

ð1Þ

where t denotes time, x and y are Cartesian coordinates, and u, f, g and s are vectors representing the conserved variables, fluxes in xand y-direction, and source terms, respectively. Ignoring Coriolis effects, viscous terms and surface stresses, these vectors are given by 2 3 2 3 f uh 6 7 6 7 u ¼ 4 uh 5; f ¼ 4 u2 h 5; 2

vh

uvh

3 vh 6 7 g ¼ 4 uvh 5

3 0 6  sbx  gh of 7 s¼4 q ox 5; s of  qby  gh oy

and

v2 h

2

ð2Þ

where f is the free surface elevation above still water depth ðhs Þ and so the total water depth h ¼ hs þ f; u and v are depth-averaged velocity components in the two Cartesian directions; g is the acceleration due to gravity; sbx and sby are bed friction stresses; and q is the density of water. The bed stress terms sbx and sby represent the energy dissipation influence of bed roughness on the flow and are estimated empirically from pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sbx ¼ qC f u u2 þ v2 and sby ¼ qC f v u2 þ v2 : ð3Þ 1=3

The bed roughness coefficient C f is evaluated using C f ¼ gn2 =h , where n is the Manning coefficient. In cases involving wetting and drying, any grid cell with h < 1:0  106 m is assumed to be dry and C f set to zero. In order to facilitate solving the shallow water equations using a of of and gh oy Godunov-type scheme, different ways of splitting the gh ox terms in (2) have been used. Taking the x-direction momentum equation as an instance, a popular way is such that 2

gh

of g oh oZ b þ gh : ¼ ox ox 2 ox

ð4Þ

After applying the same splitting technique to the y-direction momentum equation, the following hyperbolic matrix system of conservation laws for the shallow flow equations can be derived (see e.g. [17,18]): 2

h

3

6 7 u ¼ 4 uh 5;

2

6 g¼4

3

6 27 f ¼ 4 u2 h þ 12 gh 5;

vh 2

uh

uvh 3

vh uvh 2

v2 h þ 12 gh

7 5

and

2

0

3

ð5Þ

7 6 sbx   gh ozoxb 7; s¼6 5 4 q sby  q  gh ozoyb

where zb is the bed elevation above the datum and ozb =ox and ozb =oy represent bed slopes in the x- and y-directions, respectively.

223

Q. Liang, A.G.L. Borthwick / Computers & Fluids 38 (2009) 221–234

The hyperbolic conservation laws for the shallow flow equations constituted by (1) and (5) may be solved using a Godunov-type shock-capturing solver without numerical balancing of the flux gradient and source terms in cases where the bed is flat and horizontal. However, Rogers et al. [16] prove that a Godunov-type finite volume scheme based on Roe’s approximate Riemann solver cannot be used (in the absence of further treatment of the source term) to solve shallow water equations based on (5) for quiescent water conditions in a domain with non-uniform bed topography. The proof can be extended to other approximate Riemann solvers (e.g. HLLC). To overcome the balancing problem, Rogers et al. [7] proposed an alternative split of the free surface gradient terms. For example, in the x-direction, of g oðf2 þ 2fhs Þ ohs : gh ¼  gf ox ox 2 ox

ð6Þ

This, and its counterpart in the y-direction, can give a deviatoric form of the hyperbolic non-linear shallow water equations, which, as presented by Rogers et al. [16], is properly balanced when solving in a finite volume Godunov-type framework. In Rogers et al.’s equation set, the total water depth h is effectively separated into the still water depth hs and the free surface elevation f above hs . However, the free surface elevation f and still water depth hs are both difficult to define at wet–dry interfaces and in dry areas (especially where the bed topography is highly variable). Rogers et al.’s shallow water equations are therefore awkward to apply to problems involving wetting and drying. To overcome this drawback for flows over initially dry land, while retaining the mathematical balance between the flux gradient and source terms, we propose splitting the free surface gradient term in the x-direction momentum equation as follows: 2

gh

of g oðg  2zb gÞ ozb ; ¼ þ gg ox ox 2 ox

ð7Þ

and similarly in the y-direction. As shown in Fig. 1a, g is defined as the surface water level above the datum; and the water depth is then evaluated by h ¼ g  zb . As suggested by one of the anonymous reviewers, Eq. (7) can be also derived by considering the momentum balance for the vertical column of fluid shown in Fig. 1b. In the shallow water equations (1) and (2), the source terms of the momentum equations represent the external forces, which are the pressure difference and bed friction in the current formulation. In Fig. 1b, the x-direction net pressure force per unit width is 2 h2  F b ¼ 1 h1  p p

1 1 2 2  qgh1  qgh2  qg hDz b 2 2

ð8Þ

 is the mean pressure, F b is the projected component of in which p the force from the bed, Dzb is the incremental increase in bed eleva is the mean depth within the fluid coltion in the x-direction, and h umn. Dividing by the size of the fluid element, Dx, and assuming Dx ! 0, gives the net pressure force per unit plan area as

2

1 oðh Þ ozb ;  qg  qgh ox 2 ox

ð9Þ

which can be written as qgh

oðh þ zb Þ og of ¼ qgh ¼ qgh ; ox ox ox

ð10Þ

which is the x-direction pressure force component of the source terms in (1) and (2), apart from a factor of q. Substituting h ¼ g  zb , (9) can be also rearranged as 1 o ozb  qg ðg2  2zb g þ z2b Þ  qgðg  zb Þ ox 2 ox 1 o 2 ozb ¼  qg ðg  2zb gÞ  qgg : ox 2 ox

ð11Þ

Equating the expressions in (10) and (11) and dividing both sides by q, Eq. (7) is obtained. The vector terms for the shallow water equations now become 2 3 2 3 g uh 6 2 7 6 7 u ¼ 4 uh 5; f ¼ 4 u h þ 12 gðg2  2gzb Þ 5; vh 2 6 g¼4

uvh 3

vh uvh 2

v hþ

1 gðg2 2

 2gzb Þ

7 5

and

3 0 6 sbx ozb 7   gg ox 7: s¼6 5 4 q s  qby  gg ozoyb 2

ð12Þ

It is straightforward to prove that the present formulation of shallow water equations preserves a motionless steady state. Consider the x-direction momentum equation at the initial steady state. After eliminating terms with zero velocities, the equation reduces to: g oðg2  2gzb Þ ozb : ¼ gg ox 2 ox

ð13Þ

For an HLLC approximate Riemann solver based Godunov-type numerical scheme (used herein), the term on the left hand side of (13) is approximated by ðfE  fW Þ=Dx, where fE and fW are fluxes through the east and west cell interfaces of cell (i, j), respectively. In the HLLC solver, fE and fW are calculated by one of three possible formulae, depending on the local wave speed at the cell interface. However, because og=ox  0 and thus g = constant for a wet-bed application with h–0, it is trivial to prove that the fluxes are fE ¼ 2g ðg2  2gzbE Þ and fW ¼ 2g ðg2  2gzbW Þ no matter which calculation formula is used. Therefore ðfE  fW Þ=Dx ¼ ggðzbW  zbE Þ=Dx and the discretized form of (13) is balanced, provided the bed slope term at the right hand side of (13) is discretized as ggðzbE  zbW Þ=Dx. A similar analysis can be applied in the y-direction momentum equation. Therefore, Eqs. (1) and (12) are in flux gradient and source terms balanced form and can be directly applied to predict wet-bed cases of shallow flow hydrodynamics with uneven bed bathymetry without any need to use special numerical treatment for the source terms, provided that the same discretization method is applied to ozb =ox on both sides of (13). The conser-

Fig. 1. Shallow water equations: (a) definition sketch of shallow flow bed topography; (b) a vertical column of fluid.

224

Q. Liang, A.G.L. Borthwick / Computers & Fluids 38 (2009) 221–234

vation of still water state is discussed for dry-bed applications in Section 3. The flux Jacobian for Eqs. (1) and (12) is given by 2 3 ny 0 nx oF 6 2 7 uny ð14Þ ¼ 4 ðc  u2 Þnx  uvny 2unx þ vny J¼ 5; ou 2 2 uvnx þ ðc  v Þny vnx unx þ 2vny where the flux vector F ¼ fnx þ gny , with nx and ny being Cartesian pffiffiffiffiffiffi components of the unit vector in x- and y-directions; and cð¼ ghÞ is the wave celerity. The three eigenvalues of the Jacobian matrix are k1 ¼ unx þ vny þ c;

k2 ¼ unx þ vny ;

k3 ¼ unx þ vny  c:

ð15Þ

The eigenvalues are all real and they are distinct if h–0, confirming that the proposed system given by Eqs. (1) and (12) is strictly hyperbolic. 3. Wetting and drying The shallow water equations (1) and (12) have been derived for application to wetting and drying processes on a complex terrain. A simple method is used to modify the local bed slope so as not to induce spurious flow in dry areas. In discretized form, the surface elevation and bed elevation are actually constant throughout each grid cell. When a wet cell is adjacent to a dry cell, there are two cases that need to be considered. Fig. 2a illustrates a wet cell (left) next to a dry cell on the right, where hL –0, hR ¼ 0, and gR ¼ zbR > gL . In this case, there is no flux through the interface wall. Without any modification however, the fact that gR –gL means that a numerical model will predict a non-physical flux across the common cell interface. To overcome this problem, the surface elevation gR and bed level zbR for the dry cell are temporarily replaced by the following values of g0R and z0bR whenever the case in Fig. 2a is encountered g0R

¼ gR  Dg

and

z0bR

¼

g0R ;

ð16Þ

where Dg ¼ gR  gL . At the same time, the wet cell velocity component normal to the wet–dry cell interface is set to zero to ensure the no-flow condition. In Fig. 2b, the bed surface elevation of the right hand dry cell is smaller than (or equal to) that of the left hand wet cell, i.e. gR ¼ zbR 6 gL . The flow continues to flood the dry cell and there is usually no need to modify the bed slope. When the bed slope is steep, there is a possibility that more water than is actually contained in the wet cell could be computed as flowing into the dry cell, causing the water depth in the wet cell to become negative and the scheme to become unstable. To prevent instability, the

velocity components are set to zero in any cell that has dried out. Following Brufau et al. [31], water is subtracted from the adjacent cell containing most water in order to maintain mass conservation (provided the depth in the adjacent cell is at least 1:0  106 m afterwards). To maintain the front velocity components, the conserved variables uh and vh in the adjacent cell where water has been subtracted are also modified accordingly so that u and v remain the same as before. If no adjacent cell can provide sufficient water for subtraction, then the negative mass of water required for conservation is carried onto the next iteration with the cell treated as dry until the cell is sufficiently recharged or a neighbour having adequate water is identified. Mass conservation has been monitored for all of the test cases considered in this paper, and the wetting and drying method is found to ensure absolute mass conservation. It should be noted that a more sophisticated approach whereby the fluxes are modified instead of the depths would also deal with the negative depth problem and is worth further investigation. Fig. 2c depicts two adjacent dry cells with hL ¼ 0 and hR ¼ 0. In such a case, a local bed slope modification similar to that for the situation in Fig. 2a, i.e. to modify locally the bed and water surface elevation according to (16), is implemented in order to avoid the generation of a spurious flux through the dry–dry interface. By incorporating this local bed modification method, the shallow flow model automatically captures the wet–dry front, within the framework of a Godunov-type scheme. After modifying the bed and water surface elevation at the wet–dry front, og=ox  0 and thus g = constant remain true across any cell interface, and so the previous analysis for still water state preservation in a wet bed case is still valid for a case involving wet–dry fronts. 4. Dynamically adaptive quadtree grid The hyperbolic system of conservation laws given by (1) and (12) is discretized using finite volumes on dynamically adaptive quadtree grids. The generation of quadtree grids is fast and automatic, and they are easy to adapt during the computation process. The grid is created by recursive spatial decomposition of a square root cell according to the following procedure: (1) Input seeding point data to describe the boundaries of the physical flow domain and any internal areas of interest (e.g. where the bed gradients are large). (2) Rescale the physical flow domain so that it fits into an initial square cell and specify the maximum and minimum levels for the grid. (3) Divide the initial square cell into four equal quadrant cells.

Fig. 2. Local modification of bed slope for dry-bed applications: (a) hL –0, hR ¼ 0 and gR ¼ zb

R

> gL ; (b) hL –0, hR ¼ 0 and gR ¼ zbR 6 gL ; (c) hL ¼ 0, hR ¼ 0.

Q. Liang, A.G.L. Borthwick / Computers & Fluids 38 (2009) 221–234

(4) Check each cell (including newly subdivided cells) in turn, and subdivide if two or more seeding points are found inside the cell and its subdivision level is less than the maximum prescribed. (5) Carry out further cell subdivision to ensure that no cell has a side length more than twice that of its neighbours. The first step towards generating a satisfactory quadtree grid is to define a set of seeding points that adequately describe the domain boundaries. Predetermined sets of seeding points can be produced analytically for simple regular boundaries (e.g. a circle) or else read in from external data files for more complicated domains. Seeding points can also be directly inserted in areas of interest, such as regions where the initial water surface and/or bed profile have steep gradients. A primary quadtree grid is then generated by recursively subdividing the initial square cell according to the criterion that a cell is subdivided when it contains two or more seeding points and its subdivision level is less than the maximum. On the primary mesh, boundary identification is performed so that grid regularisation is only carried out with regard to cells inside or on the boundary. This process minimises the total number of cells that need to be generated and hence reduces storage requirement and overall computational cost. Boundary identification involves checking each cell in turn by drawing two lines from the cell-centre in both directions in the vertical axis (or horizontal axis) until reaching the edge of the unit square. The intersections between each of the two lines and the boundaries are recorded. Boundary perimeters are described by closed sets of seeding points, which can be transformed into continuous form by linear interpolation. For cases where the computational domain is located inside the boundary, if the numbers of intersections are odd for both lines, the cell lies within the flow domain; otherwise the cell is outside the boundary. Each cell is then given a flag to indicate whether it lies entirely outside the computational domain, in which case the cell does not have to be regularised. Rogers [32] and Liang et al. [9] give a full description of the quadtree grid generation procedure. During grid generation, grid information (e.g. neighbour references, parent and child cell pointers) is automatically stored using a hierarchical data structure composed of a linked list, which can be interrogated to locate neighbour cells during grid regularisation and adaptation. The data structure permits easy addition and removal of grid cells according to user-specified adaptation criteria at intervals during the simulation, thus providing a locally highresolution, dynamically adaptive evolving grid. The side length (ds) of a cell can be calculated via ds ¼

L 2lev

;

ð17Þ

where L is the dimension of the square domain, and lev is the subdivision level of the cell.

225

The numerical procedure for solving the discretized governing equations on a quadtree grid is locally similar to that on a uniform grid. When a cell is adjacent to neighbours of different size, simple formulae based on natural neighbour interpolation are used to allocate flow variable values on a uniform cell template arranged around the cell, for a limited number of mesh configurations [9]. Similar natural neighbour interpolation formulae are also invoked when new cells are created by grid refinement or coarsening [9]. The grid adaptation criterion used in this work is based on the root mean square value of the free surface elevation gradients, sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2  2 og og : ð18Þ þ H¼ ox oy A wet cell is subdivided whenever H is greater than a prescribed value and the subdivision level is less than a set maximum. Cell coarsening for wet cells is as follows. Each parent cell that is subdivided into four wet cells is identified, and the child cells removed if the H values of all the child cells are less than a prescribed value and the subdivision level of the child cells is greater than a set minimum. In completely dry areas, all cells having a subdivision level greater than a prescribed minimum are removed to minimise the total number of grid cells. Cells adjacent to the wet–dry interface are refined to the prescribed highest level to capture accurately the wet– dry front. The threshold values for H are case-dependent and determined by trial and error. A more general adaptation scheme, such as the one presented by Krámer and Józsa [26] is worth further investigation. 5. Numerical model Eqs. (1) and (12) are discretized spatially on the adaptive quadtree grid using a finite volume Godunov-type approach, with the HLLC approximate Riemann solver chosen to evaluate fluxes across volume interfaces. The MUSCL–Hancock method is adopted to achieve overall second-order accuracy. A non-linear slope limiter is used to damp physically meaningless oscillations at sharp fronts. 5.1. Finite volume method

o ot

The integral form of Eq. (1) is  Z  Z Z of og u dX þ s dX; þ dX ¼ ox oy X X X

ð19Þ

with all the vector terms given by (12). Applying Green’s theorem, (19) becomes I Z Z o u dX þ F dS ¼ s dX; ð20Þ ot X S X where S is the boundary of the volume X. In a Cartesian system, the surface integral term is evaluated as

Fig. 3. Fluxes through the boundary of a control volume.

226

I

Q. Liang, A.G.L. Borthwick / Computers & Fluids 38 (2009) 221–234

F dS ¼ ðFE  FW ÞDy þ ðFN  FS ÞDx;

ð21Þ

S

where Dx and Dy are the side lengths of a grid cell in the two Cartesian directions; FE , FW , FN and FS are flux function vectors directed outwards across the corresponding cell interfaces as illustrated in Fig. 3a. In a quadtree grid, any arbitrary cell under consideration may be adjacent to cells of different size. An example is shown in Fig. 3b where the cell of interest is adjacent to two smaller cells on its right. In this case, the flux across the east interface of the cell under consideration is calculated from FE ¼ FE1 þ FE2 , where FE1 and FE2 are the fluxes through the mid-point of the west interfaces of the two smaller cells. Consider the smaller cell centred at i. The flux calculation may involve flow variables at locations w, i and e. As the conserved flow variables g, uh and vh are stored at the cell-centre, the flow variables at i are available directly. The values at w and e are obtained via natural neighbour interpolation. From (20) and (21), the explicit time-marching conservative finite volume formula is uinþ1 ¼ uni 

Dt Dt ðf E  f W Þ  ðg  gS Þ þ Dtsi ; Dx Dy N

ð22Þ

where superscript n represents the time level; subscript i is the cell index; Dt is the time step; f W and f E are the fluxes through the west and east cell interfaces; and gS and gN are the fluxes through the south and north cell interfaces. 5.2. The HLLC approximate Riemann solver The HLLC approximate Riemann solver [33] is utilised because it is suited to cases involving a wet–dry interface (see e.g. [8]) and is easy to implement. Referring to the solution structure of the HLLC approximate Riemann solver given in Fig. 4, the interface flux, for example f E , is computed from 8 if 0 6 SL ; fL > > > if SM 6 0 6 SR ; f R > > : f R if 0 P SR ; where f L ¼ fðuL Þ and f R ¼ fðuR Þ are calculated from the left and right Riemann states uL and uR for a local Riemann problem. The face values are obtained from the central values of the flow variables using a piece-wise linear reconstruction (leading to a spatially second-order accurate finite volume scheme). The fluxes f L and f R correspond to the left and right sides of the middle (contact) wave; and SL , SM and SR are estimates of the speeds of the left, middle (contact) and right waves. When evaluating fluxes in the x-direction, it should be noted that the middle wave arising from the presence of the y-direction momentum equation is always a shear wave, across which the tangential velocity component v changes discontinuously

while normal velocity component u and water depth h remain constant [8]. Expressing f ¼ ½f1 f 2 f 3 T in (12), the fluxes in the middle region may be denoted f  ¼ ½f1 f 2 f 3 T . Either side of the middle wave, the middle region fluxes f L and f R are different, and can be evaluated from 2 3 2 3 f1 f1 6 7 6 7 ð24Þ f L ¼ 4 f2 5 and f R ¼ 4 f2 5; vL f1

vR f1

where vL and vR are the left and right tangential velocity components of the Riemann states, which remain unchanged across the left and right waves, respectively. The fluxes f  in the middle region are calculated from the HLL formula proposed by Harten et al. [34] f ¼

SR f L  SL f R þ SL SR ðuR  uL Þ : SR  SL

ð25Þ

From (24), it is obvious that only f1 and f2 are calculated from the above HLL formula. The third flux component in the middle flux vector f  is calculated by f3 ¼ vL f1 or vR f1 depending on which part of the middle region the flux component is evaluated. Different choices of wave speed estimates are possible. Fraccarollo and Toro [3] suggest the following approximations including dry-bed options from the two-rarefaction approximate Riemann solver, whereby the left and right wave speeds are determined from ( pffiffiffiffiffiffiffiffi uR  2 ghR if hL ¼ 0; ð26Þ SL ¼ pffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffi minðuL  ghL ; u  gh Þ if hL > 0; and ( SR ¼

pffiffiffiffiffiffiffiffi uL þ 2 ghL pffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffi maxðuR þ ghR ; u þ gh Þ

if hR ¼ 0; if hR > 0;

where uL , uR , hL and hR are the components of the left and right Riemann states for a local Riemann problem, and are calculated from qffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffi 1 ð27Þ u ¼ ðuL þ uR Þ þ ghL  ghR ; 2 and h ¼

 qffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffi 2 1 1 1 ð ghL þ ghR Þ þ ðuL  uR Þ : g 2 4

ð28Þ

For general problems including dry-bed cases, Toro [8] recommends calculating the middle wave speed from: SM ¼

SL hR ðuR  SR Þ  SR hL ðuL  SL Þ : hR ðuR  SR Þ  hL ðuL  SL Þ

ð29Þ

Similar formulae are used to calculate f W , gS and gN . 5.3. Discretisation of source terms Eqs. (1) and (12) are mathematically conditioned so that the discretised flux gradient and source terms balance, and so there is no need for any special numerical treatment of the source terms. Herein, the source terms are cell-centred with the bed slope terms ozb =ox and ozb =oy approximated by second-order accurate central differences. The x-direction source terms are calculated from  n z  z  sbx ozb sbx bE bW  gg ¼  ggni  ; ð30Þ q ox q i Dx

Fig. 4. HLLC solution structure of a Riemann problem.

where zbE and zbW are the bed elevations at the mid-point of the eastern and western interfaces of cell i. It should be noted that the cell-centred flow variables obtained from the predictor step are used when computing the source terms in the corrector step within the MUSCL-Hancock framework described as follows.

227

Q. Liang, A.G.L. Borthwick / Computers & Fluids 38 (2009) 221–234

5.4. The MUSCL-Hancock method

In the predictor step, for the west and east face values of cell i, (33) becomes

Overall second-order accuracy in space and time is achieved using a two-step unsplit MUSCL-Hancock method. The predictor step comprises  inþ1=2 ¼ uni  u

Dt Dt Dt ðf E  f W Þ  ðg  gS Þ þ si : 2Dx 2Dy N 2

ð31Þ

The interface fluxes f W , f E , gS and gN are calculated from the face values at mid-points of cell faces, which are linearly interpolated from the corresponding cell-centre values of the flow variables. No Riemann solutions are required at the predictor step and the interface fluxes are essentially evaluated within each cell using interpolated values at the extremities of the inner interfaces. In the finite volume discretization of the governing equations, the flow states are usually different either side of an interface (interpolated from different cell-centre values), which implies that the resulting fluxes are distinct either side of the interface. Hence, the use of inner fluxes to represent the overall fluxes through the cell face in the predictor step is non-conservative. However, as this is merely an intermediate step it does not affect the conservative character of the overall numerical scheme. In the corrector step, the HLLC approximate Riemann solver is used to calculate interface fluxes and the flow variables are updated over a full time step using Eq. (22), which is fully conservative [35]. The Riemann states are calculated by a piece-wise linear reconstruction method based on intermediate cell-centre values obtained in the predictor step. 5.5. Slope limiter Face values of the flow variables, required in both the predictor and corrector steps of the MUSCL-Hancock scheme, are calculated using linear interpolation: uðx; yÞ ¼ ui þ rrui ;

uðx; yÞ ¼ ui þ WðrÞrrui ;

ð33Þ

in which the slope limiter WðrÞ is a function of the ratio of gradients defined as ðuhÞe  ðuhÞi r¼ ðuhÞi  ðuhÞw

and

ðvhÞe  ðvhÞi r¼ ; ðvhÞi  ðvhÞw

WðrÞ ¼ max½0; minðbr; 1Þ; minðr; bÞ;

ð35Þ

in which bð1 6 b 6 2Þ is the limiter parameter. b ¼ 1 gives the minmod limiter and b ¼ 2 leads to Roe’s superbee limiter. A minmod limiter, used throughout all the test cases in this work, tends to be more dissipative but gives smoother results than a superbee limiter. Although not considered in the present paper, there is scope for a parameter study to determine the value for b that achieves a sensible compromise between the minmod and superbee limiters. This is recommended for future work.

1 WðrÞðuni  unw Þ: 2

ð36Þ

1 þ WðrÞðuni  unw Þ; 2 1  enþ1=2  WðrÞðune  uni Þ: uRE ¼ u 2 nþ1=2

i uLE ¼ u

ð37Þ

 nþ1=2  enþ1=2 are the predicted cell-centre values at cell i where u and u i and its eastern neighbour obtained in the predictor step. When calculating the Riemann states in the corrector step, Hu et al. [37] suggest that use of the original flow variable values un for the limiter and gradient vector terms WðrÞðuni  unw Þ=2 and WðrÞðune  uni Þ=2 gives better results than the use of intermediate flow variables  nþ1=2 . The Riemann states at the other three cell faces are calcuu lated in a similar way. 5.6. Stability criteria The numerical scheme is explicit, and its stability is governed by the Courant–Friedrichs–Lewy (CFL) criterion. For a two-dimensional Cartesian grid, the CFL criterion for choosing an appropriate time step, Dt, may be expressed Dt ¼ C minðDt x ; Dt y Þ;

ð38Þ

with Dtx ¼ min i

Dxi pffiffiffiffiffiffiffi j ui j þ ghi

and

Dty ¼ min i

Dyi pffiffiffiffiffiffiffi ; j vi j þ ghi

where C is the Courant number specified in the range 0 < C 6 1. In the present scheme, an adaptive time step is used based on (38) with C set to 0.8. The overall time step at each iteration is generally constrained by the smallest grid cells, but this is offset by the improved local cell distribution in the dynamically adaptive quadtree grid. Higher computational efficiency may be achievable by implementing a more sophisticated local time step method (see e.g. [26,38]). 5.7. Boundary conditions

ð34Þ

for the three conserved flow variables in the x-direction. Herein, the subscripts ‘e’ and ‘w’ denote the eastern and western neighbours of cell i. A similar approach is used in the y-direction. When the cell under consideration has neighbours of different size, as illustrated in Fig. 3b, a locally uniform cell template is used in conjunction with natural neighbour interpolation. The slope limiter does not therefore require local modification. The slope limiter [36] is defined as

unE ¼ uni þ

These formulae are for discretization on a uniform cell template. The south and north face values are obtained in a similar fashion. In the corrector step, the fluxes through the cell interfaces are evaluated by the HLLC approximate Riemann solver. At the eastern interface of cell i, the left and right Riemann states are evaluated from

ð32Þ

where (x, y) is the point where the interpolants are required (i.e. the mid-point of the cell faces); r is the distance vector from the cellcentre to point (x, y) with east and north being the positive in the x- and y-direction, respectively; and rui is the gradient vector. A slope limiter is used to prevent the numerical oscillations of the solution in the vicinity of steep gradients or discontinuities. Eq. (32) becomes

g  gi r¼ e ; gi  gw

1 unW ¼ uni  WðrÞðuni  unw Þ; 2

At solid (slip) boundaries: hB ¼ hI ;

^ B ¼ 0; u

^B ¼ v ^I : v

ð39Þ

If the boundary is open involving inflow and outflow, the boundary conditions are implemented based on Riemann invariants according to the local Froude number: (1) Fr < 1 (subcritical flow) pffiffiffi pffiffiffiffi pffiffiffiffiffi ^B ¼ v ^ I  2 g ð hI  hB Þ; v ^I ; where hB is prescribed ^B ¼ u u ð40Þ or pffiffiffiffi  1 ^B Þ ; ^I  u hI  pffiffiffi ðu 2 g ^I ; where u ^ B is prescribed; ^B ¼ v v

^B ¼ h

ð41Þ

with the sign dependent on the flow direction: outflow +; inflow .

228

Q. Liang, A.G.L. Borthwick / Computers & Fluids 38 (2009) 221–234

(2) Fr > 1 (supercritical flow) ^ B , and v ^B are prescribed; and for For inflow, the variables hB , u outflow, hB ¼ hI ;

^B ¼ u ^I ; u

^B ¼ v ^I : v

ð42Þ

^ and v ^ are the depth-averaged velocity In the above equations, u components normal and tangential to the boundary. Subscripts B and I denote boundary values and inner values next to the boundary. 6. Results and discussion The numerical scheme is validated against several benchmark tests and results are compared with analytical solutions and alternative numerical predictions where available. The test cases considered include 1D steady flow over a hump, the interaction of a bore-like disturbance with an elliptical-shaped hump, liquid sloshing in a container with a parabolic bottom and a dam break flood wave travelling over an otherwise dry bed topography containing three large mounds. In all cases, g = 9.81 m/s2 and q = 1000 kg/m3. 6.1. One-dimensional steady flow over a hump This benchmark test proposed by Goutal and Maurel [39] validates the treatment of source terms (see e.g. [15,16,40]) as it involves spatially varying bed topography. The flow occurs in a 25 m long frictionless channel whose bed elevation is defined by ( 0:2  0:05ðx  10Þ2 if 8 < x < 12; ð43Þ zb ðxÞ ¼ 0 otherwise:

Depending on the initial and boundary conditions, the steady flow over the hump can be subcritical, transcritical with or without a hydraulic jump, or supercritical. Analytical solutions are suggested by Goutal and Maurel [39]. All the cases considered by Goutal and Mourel have been simulated, but results are only presented for steady transcritical flow with a jump and subcritical flow, which are most challenging to a numerical model. The computational channel is 5 m wide and uniform in the y-direction, with slip boundary conditions imposed at the north and south walls of the channel. However, the sidewall boundary conditions are irrelevant as the flow is essentially one-dimensional. Following Zhou et al. [15], convergence of the numerical solution is indicated by the global relative error, vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! u uX hn  hn1 2 i i ; ð44Þ R¼t n hi i n

n1

where hi and hi are water depths at the current and previous time steps at cell i. The solution is taken as having converged to steady state if R < 1  106 . For the case of transcritical flow with a jump, we impose a constant discharge per unit width q = 0.18 m2/s at the upstream boundary and a fixed water depth of h = 0.33 m at the outflow boundary. The numerical solution is obtained on a quadtree grid with maximum and minimum subdivision levels of 8 and 6, commencing from an initial grid created about seeding points located along the channel walls and inflow and outflow boundaries, with the adaptation criteria set to H ¼ 0:01 for refinement and H ¼ 0:005 for coarsening. At the start of the simulation, the surface elevation and discharge per unit width in the x-direction are set to g = 0.33 m and q = 0.18 m2/s at each grid point. Fig. 5 shows the

Water Surface Elevation

Discharge 0.26

0.4 0.22

q (m2/s)

η (m)

0.3 0.2 Bed profile Analytical Numerical

0.1 0 0

5

10

15

20

0.18 0.14 Theoretical discharge Numerical discharge

0.1 0.06 0

25

5

10

15

20

25

x (m)

x (m)

Fig. 5. Steady flow over a hump – transcritical flow with a shock: (a) water surface elevation; (b) discharge.

14000

1 0.8 Relative error

Cell number

12000 10000 8000 6000

0.4 0.2

4000 2000

0.6

0

1000

2000 3000 4000 Iteration number

5000

6000

0 0

2000

4000 6000 8000 10000 12000 Iteration number

Fig. 6. Steady flow over a hump – transcritical flow with a shock: (a) temporal change of cell number; (b) temporal evolution of relative error.

229

Q. Liang, A.G.L. Borthwick / Computers & Fluids 38 (2009) 221–234

general agreement between the predicted and analytical solutions at steady state in terms of the water surface elevation and discharge per unit width along the central line of the channel in the x-direction. For the discharge, there is a maximum difference of about 15% between the numerical prediction and analytical solution at a grid point near the middle of the domain. Similar deviations between the predicted and analytical solutions have also been observed by other researchers. For example, Valiani and Begnudelli [40], Zhou et al. [15] and Rogers et al. [16], report errors in discharge of about 23%, 10% and 15%. Fig. 6 presents the changes in number of quadtree leaf cells and relative error (44) against number of time steps. The number of leaf cells grows from an initial value of 2380 to a maximum of 13,210 after about 800 iterations, then settles to a steady value of 4252 after about 5200 iterations. The solution converges to steady state after about 12,000 iterations. Fig. 7 illustrates the effect of grid size, showing that the

Water Surface Elevation

0.3

η (m)

6.2. Shallow water disturbances propagating over an elliptical-shaped hump Extending the validation to two-dimensions, we next consider the case of shallow flow over an isolated elliptical-shaped hump [11] in a 2 m long, 1 m wide rectangular domain, with the hump bed elevation given by 2 2 zb ðx; yÞ ¼ 0:8 expð5ðx  0:9Þ  50ðy  0:5Þ Þ;

0.4

Bed profile Analytical lev = 8 lev = 7 lev = 6

0.2 0.1 0

numerical predictions from three uniform quadtree grids of different levels (level 6, 7 and 8) converge to the analytical solution. In the case of subcritical flow, the upstream inflow discharge per unit width is set to a constant value of q = 4.42 m2/s and the downstream water depth is fixed at 2 m. These values are also used to provide initial values for the flow variables throughout the domain. The same grid generation parameters are used as for the previous case. Fig. 8 shows that the numerical predictions match the analytical solutions of water surface elevation and discharge along the centreline of the channel. Fig. 9 presents the temporal changes in the number of cells and relative error against iteration number.

0

5

10

15

20

25

x (m) Fig. 7. Steady flow over a hump – transcritical flow with a shock: grid convergence.

where x and y are horizontal Cartesian distances with origin at the southwest corner of the domain. Initially the free surface of the water throughout the domain is flat and 1 m above datum, with h ¼ 1  zb . In the region 0:05 < x < 0:15, the free surface elevation is then disturbed so that it is raised to 1.01 m but remains otherwise still at t = 0 s. The bed is frictionless and the lateral wall boundaries are open. An initial quadtree grid, with maximum and minimum subdivision levels of 9 and 6, is generated with the finest mesh at the open boundary walls and along the line with discontinuous water surface. Fig. 10 shows the initial grid and bed topography contours for the problem domain. For grid adaptation, the free

Water Surface Elevation

Discharge 4.48

2 4.45

q (m2/s)

η (m)

1.6 1.2

Bed profile Analytical Numerical

0.8

4.42 4.39 4.36

0.4

Theoretical discharge Numerical discharge

4.33

0 0

5

10

15

20

4.3 0

25

x (m)

5

10

15

20

25

x (m)

Fig. 8. Steady flow over a hump – subcritical flow: (a) water surface elevation; (b) discharge.

7000 0.015 Relative error

Cell number

6000 5000 4000

0.01

0.005

3000 2000

0

500

1000 1500 Iteration number

2000

ð45Þ

0 0

1000 2000 3000 4000 5000 6000 7000 Iteration number

Fig. 9. Steady flow over a hump – subcritical flow: (a) temporal change of cell number; (b) temporal evolution of relative error.

230

Q. Liang, A.G.L. Borthwick / Computers & Fluids 38 (2009) 221–234

Fig. 10. Shallow water wave propagating over an elliptical-shaped hump: initial quadtree grid and bed topography contours.

surface enrichment and coarsening thresholds are H ¼ 0:01 and H ¼ 0:005. Fig. 11 presents the predicted wave patterns at different non-dimensional times, in terms of three-dimensional free surface elevation, water surface contours and the associated adapted grids

as the free surface disturbance interacts the submerged hump. pffiffiffiffiffiffiffiwith ffi The non-dimensional time is ~t ¼ t= gh0 (consistent with e.g. [11]) where h0 ¼ 1m denotes the undisturbed initial water depth away from the hump. The step-like perturbation of the free surface immediately releases two soliton-like waves, one directed westward, the other eastward. At ~t ¼ 0:6, the westward wave has already exited the domain. The eastward wave distorts as it approaches the hump, and there are spatial variations in the speed of the wave whose free surface moves more slowly in water of shallower depth over the hump (Fig. 11a). As the wave crest passes over the top of the hump at ~t ¼ 0:9, the local free surface elevation rises to a peak where energy is concentrated (Fig. 11b). Wave reflection and scattering due to energy transfer then occur at the hump (see e.g. Fig. 11c at ~t ¼ 1:2). The front of the initial eastward-directed wave gradually straightens out as it travels downstream (Fig. 11d and e at ~t ¼ 1:5 and 1:8). The wave pattern is complicated due to wave– wave interactions between the scattered and incident waves. The

Fig. 11. Shallow water wave propagating over an elliptical-shaped hump: 3D surface, surface elevation contours and adapted quadtree grids at different output times.

231

Q. Liang, A.G.L. Borthwick / Computers & Fluids 38 (2009) 221–234

results agree closely with those obtained by LeVeque [11] who used a high resolution Godunov-type method. 6.3. Moving shorelines on parabolic bottom topography Sampson et al. [41] derived analytical solutions of the non-linear shallow water equations for perturbed flow in a container with parabolic bed topography with bed friction, neglecting the Coriolis effect. This benchmark test is useful for validating the wetting and drying procedure in the numerical model, as well as providing a further check that the source terms (i.e. bed friction and bed slopes) are being evaluated correctly. The bed profile of the domain is defined by 2

zb ðxÞ ¼ h0 ðx=aÞ ;

ð46Þ

where h0 and a are constants and zb is uniform with zero slope in the y-direction. Sampson et al. present analytical solutions depending on a bed friction parameter s (related to the bed friction coeffipffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2 2 cientpffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi by C f ¼ ffi hs= u þ v Þ and a hump amplitude parameter p ¼ 8gh0 =a2 . We consider a case where s < p, and the analytical solution of the water surface time history above a given horizontal datum is

a2 B2 est ðss sin 2st þ ðs2 =4  s2 Þ cos 2stÞ 8g 2 h0   B2 est est=2 sB   Bs cos st þ sin st x; ð47Þ 4g g 2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where B is a constant and s ¼ p2  s2 =2. As t ! 1, gðtÞ ! h0 , defining the still water elevation above the datum. The projection of the moving shorelines (two parallel straight lines on the x—y plane) is given by   a2 est=2 sB x¼ Bs cos st  sin st  a: ð48Þ 2gh0 2 gðx; tÞ ¼ h0 þ

As t ! 1, x ! a, indicating that the oscillatory flow is damped by bed friction. The computational domain has plan dimensions of 10,000 m  1000 m. Slip conditions are invoked at the lateral walls. The coefficients are a = 3 km, h0 ¼ 10 m, s = 0.001 s1, and B = 5 m/s. The surface gradient is spatially smooth, and so the simulation lasting 6000 s is carried out on a uniform level 7 quadtree grid. Fig. 12 shows the excellent agreement between the numerical predictions and analytical solutions of the free surface profile at intervals throughout the simulation. The moving wet–dry inter-

t = 1000 (s)

t = 500 (s) 30 bed profile analytical solution numerical prediction

25

Surface Elevation (m)

Surface Elevation (m)

30

20 15 10 5 0

4000

2000

0 x (m)

2000

20 15 10 5 0

4000

bed profile analytical solution numerical prediction

25

4000

2000

t = 1500 (s)

Surface Elevation (m)

Surface Elevation (m)

bed profile analytical solution numerical prediction

20 15 10 5 4000

2000

0 x (m)

2000

bed profile analytical solution numerical prediction

25 20 15 10 5 0

4000

4000

2000

t = 3000 (s)

0 x (m)

2000

4000

t = 6000 (s)

30

30 bed profile analytical solution numerical prediction

25

Surface Elevation (m)

Surface Elevation (m)

4000

30

25

20 15 10 5 0

2000

t = 2000 (s)

30

0

0 x (m)

4000

2000

0 x (m)

2000

4000

bed profile analytical solution numerical prediction

25 20 15 10 5 0

4000

2000

0 x (m)

Fig. 12. Sloshing motions in a vessel with parabolic bottom topography.

2000

4000

232

Q. Liang, A.G.L. Borthwick / Computers & Fluids 38 (2009) 221–234

faces are correctly reproduced, thus validating the wetting and drying algorithm and the discretization of the source terms. 6.4. Dam-break wave propagating over three humps We finally consider the case of a dam-break wave travelling over an initially dry floodplain with three humps, which was proposed by Kawahara and Umetsu [42] and has been reconsidered by other researchers (e.g. [17]). The dam break occurs in a 75 m long and 30 m wide closed container, with the bed topography defined by  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 ðx  30Þ2 þ ðy  6Þ2 ; zb ðx; yÞ ¼ max 0; 1  8 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 ðx  30Þ2 þ ðy  24Þ2 ; 1 8 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 ð49Þ 3 ðx  47:5Þ2 þ ðy  15Þ2 : 10

Fig. 13. Dam-break wave propagating over three humps: initial state.

Initially, there is a dam located at x = 16 m that retains still water with surface elevation g = 1.875 m. Fig. 13 illustrates the bed topography within the computational domain and the upstream still water retained by the dam. At t ¼ 0, the dam collapses instantaneously. The bed roughness coefficient is evaluated with the Manning coefficient set to n ¼ 0:018. The initial quadtree grid has maximum and minimum levels of 8 and 6 and is generated about seeding points that delineate the lateral walls of the container. During adaptation, grid enrichment and coarsening are undertaken according to H ¼ 0:08 and H ¼ 0:05. Slip conditions are applied at the lateral walls of the container. Fig. 14 illustrates the inundation of the floodplain at different times after the dam break has occurred until steady state is reached. 3D visualisations of the water surface, plan views of water depth plot, and corresponding adapted quadtree grids are presented. Immediately after the dam failure, a flood wave led by a wet–dry front starts to inundate the floodplain. By t = 2 s, the wet–dry front has reached the pair of small hills and begun to rise over them. Curved bore-like reflection waves are created at the small hills and start to propagate upstream, while part of the wet–dry front continues downstream between the hills. At t = 6 s, the small hills are entirely submerged, and the wet–dry flood front has reached the large hill. Due to the momentum of the flood wave, part of the front runs more than halfway up the large hill, while water cascades around the side of the hill. The curved reflection bores continue to move upstream and begin to interact with each other and the side walls of the container. At t = 12 s, the flood water that is passing either side of the large hill starts to flood

Fig. 14. Dam-break wave propagating over three humps: 3D water surface elevation, plan views of water depth plot and adapted quadtree grids at different output times.

Q. Liang, A.G.L. Borthwick / Computers & Fluids 38 (2009) 221–234

istic flood inundation cases involving locally complex flows, such as trans-critical flows over complicated terrain. It should also be noted that the quadtree grid is generated by recursive domain decomposition and so is inherently adept at representing real topographic data supplied from a digital terrain model in G.I.S.

12500

Number of cells

10000

7500

Acknowledgements

5000

2500 0

233

50

100

150 t (s)

200

250

300

Fig. 15. Dam-break wave propagating over three humps: time history of total number of cells in the dynamically adaptive quadtree grid.

the lee of the hill. The upstream directed reflection shock-waves from the small hills have coalesced into a nearly straight bore. The main flood wave has reflected against the large hill, releasing a further upstream directed bore. The maximum water level at the large hill reduces as the reflection bore moves away. After further complicated wave–wave, wave–wall and wave–topography interactions and energy dissipation due to bed friction, the flow starts to approach steady state. By t = 300 s, steady state is achieved, with the peaks of the small hills no longer submerged. The numerical model properly simulates complicated wetting and drying processes and produces results that are very similar to those of Brufau et al. [17]. Absolute conservation of mass is also guaranteed for this test case. The adaptive quadtree grid evolves dynamically according to the development of the dam-break wave, with the finest meshing at the wet–dry front and where the flow pattern is locally rapidly changing. Fig. 15 shows the temporal evolution of the number of cells in the adaptive quadtree grid, which initially increases from 3042 to a peak of 12,216 and then decreases to a steady value of about 3500 as the flow becomes less dynamic. This highlights the substantial efficiency gain to be had by using dynamically adaptive grids for flood inundation problems. 7. Conclusions This paper presents details of a Godunov-type numerical model that solves the two-dimensional deviatoric shallow water equations derived as conservation laws with the flux and source terms conveniently balanced so that there is no need for further numerical treatment. The hyperbolic non-linear shallow water equations have been formulated so that flood inundation on otherwise dry terrain can be readily simulated using a Godunov-type finite volume solver that balances the flux gradient and source terms due to bed friction and bed slope. A second-order Godunov-type finite volume scheme is used to solve the governing equations on dynamically adaptive quadtree grids. Numerical wetting and drying is based on the local modification of bed slope, and so is easy to implement. The equation set together with the approach dealing with wetting and drying conserves still water conditions for both wet-bed and dry-bed applications. The model gives accurate results in agreement with analytical solutions and alternative numerical predictions for flow over humps in 1D and 2D, and liquid sloshing in a parabolic-bottomed container. The use of dynamic grid adaptation gives enhanced computational efficiency while retaining accuracy at sharp fronts, as is evident in the case of the flood inundation over an initially dry floodplain containing three humps. The present results indicate that the scheme works well for the range of laboratory scale tests considered. It is intended that the numerical model will be applied in future to real-

The authors thank Dr. J.G. Zhou from the University of Liverpool for his useful discussion on the test case of steady flow over a bump. The authors are grateful to the anonymous reviewers for many useful suggestions, including the momentum balance argument used as an alternative derivation of the net pressure forces. References [1] Alcrudo F, García-Navarro P. A high-resolution Godunov-type scheme in finite volumes for the 2D shallow-water equations. Int J Numer Method Fluids 1993;16:489–505. [2] Zhao DH, Shen HW, Tabious III GQ, Lai SJ, Tan WY. Finite-volume twodimensional unsteady-flow model for river basins. J Hydraul Eng ASCE 1994;120(7):864–83. [3] 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(6):843–63. [4] Anastasiou K, Chan CT. Solution of the 2D shallow water equations using the finite volume method on unstructured triangular meshes. Int J Numer Method Fluids 1997;24:1225–45. [5] Sleigh PA, Gaskell PH, Berzina M, Wright NG. An unstructured finite-volume algorithm for predicting flow in rivers and estuaries. Comput Fluids 1998;27(4):479–508. [6] Causon DM, Ingram DM, Mingham CG, Yang G, Pearson RV. Calculation of shallow water flows using a Cartesian cut cell approach. Adv Water Resour 2000;23:545–62. [7] Rogers BD, Fujihara M, Borthwick AGL. Adaptive Q-tree Godunov-type scheme for shallow water equations. Int J Numer Method Fluids 2001;35:247–80. [8] Toro EF. Shock-capturing methods for free-surface shallow flows. Chichester: John Wiley & Sons; 2001. [9] Liang Q, Borthwick AGL, Stelling G. Simulation of dam and dyke-break hydrodynamics on dynamically adaptive quadtree grids. Int J Numer Method Fluids 2004;46:127–62. [10] Bermudez A, Vázquez ME. Upwind methods for hyperbolic conservation laws with source terms. Comput Fluids 1994;23(8):1049–71. [11] LeVeque RJ. Balancing source terms and flux gradients in high-resolution Godunov methods: the quasi-steady wave-propagation algorithm. J Comput Phys 1998;146(1):346–65. [12] Vázquez-Cendón ME. Improved treatment of source terms in upwind schemes for the shallow water equations in channels with irregular geometry. J Comput Phys 1999;148(2):497–526. [13] Hubbard ME, García-Navarro P. Flux difference splitting and the balancing of source terms and flux gradients. J Comput Phys 2000;165(1):89–125. [14] García-Navarro P, Vázquez-Cendón ME. On numerical treatment of the source terms in the shallow water equations. Comput Fluids 2000;29:951–79. [15] 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):1–25. [16] Rogers BD, Borthwick AGL, Taylor PH. Mathematical balancing of flux gradient and source terms prior to using Roe’s approximate Riemann solver. J Comput Phys 2003;192(2):422–51. [17] Brufau P, García-Navarro P, Vázquez-Cendón ME. A numerical model for the flooding and drying of irregular domains. Int J Numer Method Fluids 2002;39:247–75. [18] Begnudelli L, Sanders BF. Unstructured grid finite-volume algorithm for shallow-water flow and scalar transport with wetting and drying. J Hydraul Eng ASCE 2006;132(4):371–84. [19] Bates PD, de Roo APJ. A simple raster-based model for flood inundation simulation. J Hydrol 2000;236:54–77. [20] Yu D, Lane SN. Urban fluvial flood modelling using a two-dimensional diffusion-wave treatment, part 1: mesh resolution effects. Hydrol Process 2006;20:1541–65. [21] Yu D, Lane SN. Urban fluvial flood modelling using a two-dimensional diffusion-wave treatment, part 2: development of a sub-grid-scale treatment. Hydrol Process 2006;20:1567–83. [22] Ivanenko SA, Muratova GV. Adaptive grid shallow water modelling. Appl Numer Math 2000;32:447–82. [23] Hubbard ME, Dodd N. A 2D numerical model of wave run-up and overtopping. Coast Eng 2002;47:1–26. [24] Skoula ZD, Borthwick AGL, Moutzouris CI. Godunov-type solution of the shallow water equations on adaptive unstructured triangular grids. Int J Comput Fluid Dynam 2006;20(9):621–36.

234

Q. Liang, A.G.L. Borthwick / Computers & Fluids 38 (2009) 221–234

[25] Liang Q, Zang J, Borthwick AGL, Taylor PH. Shallow flow simulation on dynamically adaptive cut-cell quadtree grids. Int J Numer Method Fluids 2007;53(12):1777–99. [26] Krámer T, Józsa J. Solution-adaptivity in modelling complex shallow flows. Comput Fluids 2007;36(3):562–77. [27] George PL. Automatic mesh generation: application to finite element methods. Chichester, UK: Wiley; 1991. [28] Thompson JF, Soni BK, Weatherill NP, editors. Handbook of grid generation. Boca Raton, USA: CRC Press; 1999. [29] Lamby P, Müller S, Stiriba Y. Solution of shallow water equations using fully adaptive multiscale schemes. Int J Numer Method Fluids 2005;49:417–37. [30] Falconer RA. An introduction to nearly horizontal flows. In: Abbott MB, Price WA, editors. Coastal, estuarial and harbour engineer’s reference book. London: Chapman and Hall; 1993. p. 29–36. [31] Brufau P, García-Navarro P, Vázquez-Cendón ME. Zero mass error using unsteady wetting–drying conditions in shallow flows over dry irregular topography. Int J Numer Method Fluids 2004;45:1047–82. [32] Rogers BD. Refined localised modelling of coastal flow features using adaptive quadtree grids. D.Phil. thesis, University of Oxford; 2001. [33] Toro EF, Spruce M, Speares W. Restoration of the contact surface in the HLL Riemann solver. Shock Waves 1994;4(1):25–34.

[34] Harten A, Lax PD, van Leer B. On upstream differencing and Godunov-type schemes for hyperbolic conservation-laws. SIAM Rev 1983;25(1):35–61. [35] Toro EF. Riemann solvers and numerical methods for fluid dynamics: a practical introduction. Berlin/New York: Springer; 1997. [36] Hirsch C. Numerical computation of internal and external flows. Computational methods for inviscid and viscous flows, vol. 2. New York, USA: Wiley; 1990. [37] Hu K, Mingham CG, Causon DM. Numerical simulation of wave overtopping of coastal structures using the non-linear shallow water equations. Coast Eng 2000;41(4):433–65. [38] Crossley AJ, Wright NG. Time accurate local time stepping for the unsteady shallow water equations. Int J Numer Method Fluids 2005;48(7):775–99. [39] Goutal N, Maurel F, editors, Proceedings of the 2nd workshop on dam-break wave simulation. HE 43/97/016/B, Départment Laboratoire National d’Hydraulique, Groupe Hydraulique Fluviale Electricité de France, France; 1997. [40] Valiani A, Begnudelli L. Divergence form for bed slope source term in shallow water equations. J Hydraul Eng ASCE 2006;132(7):652–65. [41] Sampson J, Easton A, Singh M. Moving boundary shallow water flow above parabolic bottom topography. ANZIAM J 2006;47(EMAC2005):C373–87. [42] Kawahara M, Umetsu T. Finite element method for moving boundary problems in river flow. Int J Numer Method Fluids 1986;6:365–86.