A TWO-DIMENSIONAL UNSTRUCTURED SOLUTION ALGORITHM APPLIED TO THE SIMULATION OF DAM-BREAK WAVES

 

P Andrew Sleigh^? and Philip H Gaskell^?

 

^? School of Civil Engineering

^? School of Mechanical Engineering

University of Leeds, Leeds, LS2 9JT, UK

Correspondence address: Dr Andrew Sleigh

School of Civil Engineering, University of Leeds, Leeds, LS2 9JT, UK

Tel: +44 (0)113 2332398, Fax: +44 (0)113 2332265, Email: P.A.Sleigh@leeds.ac.uk

 

 

ABSTRACT

A numerical algorithm is presented for the solution of geometrically challenging two-dimensional river flows. It is capable of simulating smooth and discontinuous transient flow conditions on complex shaped domains with complex topography.

The solution algorithm possesses the following properties:

·        Is centred on the conservative form of the shallow water equations.

·        Employs a finite volume formulation on unstructured triangular meshes;

·        Uses an approximate Riemann solver to calculate the convective fluxes together with an associated flux limiter that is monotonicity preserving.

·        Makes use of second order spatial and temporal discretisation.

·        Embodies the scope for automatic spatial and temporal error control.

·        Has the flexibility to specify general boundary conditions

The algorithm is applied to the problem of modelling a dam-break wave down a laboratory channel. Results are presented and compared with experimental measurements. The numerical solutions are shown to simulate the passages of the wave with good accuracy.

 

INTRODUCTION

The numerical algorithm presented here was developed for the solution of geometrically challenging two-dimensional river and estuary flows. In particular it was designed to: handle both steady and transient flow conditions (smooth or otherwise); describe and incorporate complex topography; simulate both sub- and super-critical conditions and to accommodate flow around and through structures such as weirs and gates. More exacting design criteria involve correct embodiment of the underlying physics; quantifiable accuracy and speed of computation.

 

The governing equations are the shallow water equations which are discretised using a finite volume approach embodying variable step time integrators to yield a method that is second order accurate in both space and time. The approach is capable of handling complex flow domains and yielding solutions for which errors are controlled automatically by the use of spatial re-gridding and time stepping based on local error estimates. Its range of applicability is thoroughly demonstrated in [1] through considering several problems involving super/sub-critical flow, wetting/drying, culminating in the solution of a complete estuary problem. An approximate Riemann solver is used to determine flux across both cell and physical boundaries enabling a straightforward and consistent inclusion of any boundary condition.

 

In summary the solution algorithm possesses the following properties:

·        Is centred on the conservative form of the shallow water equations.

·        Employs a finite volume formulation, the unstructured computational grid for which is generated automatically by triangular decomposition of the region of interest; any associated digital bathymetry data is similarly interpolated automatically.

·        Uses an approximate Riemann solver to calculate the convective fluxes together with an associated flux limiter that is monotonicity preserving.

·        Makes use of second order spatial and temporal discretisation.

·        Embodies the scope for automatic spatial and temporal error control.

·        Has the flexibility to specify general boundary conditions, including the capacity to handle wetting/drying characteristics.

 

In this paper the method is applied to the problem of dam-break wave modelling. The method is particularly suited to the solution of this problem due to its use of Riemann solvers to model the steep gradients and shocks which frequently occur.

 

GOVERNING EQUATIONS AND DOMAIN DECOMPOSITION

The model is based on the two-dimensional shallow water equations, containing source terms representing frictional stress and momentum change due to a sloping bed, which when written in conservative form can be expressed as

 

 

where

 

Equation 1

 

and u, and v are the velocities in the x and y directions respectively. , where h is the depth and g the acceleration due to gravity. K is a conveyance term that can be calculated from one of the common friction formulae, for example Manning of Chezy.

 

MESH GENERATION

The type of mesh employed is an unstructured triangulation of the solution domain enabling arbitrary shaped geometries to be accommodated more easily than with a square grid system. Many algorithms are available that generate automatically an initial mesh with the required properties. The property most desired is that the mesh is smooth, that is, that adjacent triangles do not have greatly differing qualities and that none of the triangles have a quality less than a prescribed value. Use was made of the GEOMPACK [2] mesh generator. GEOMPACK allows the mesh to be distributed as required via a function defining a weight between 0 and 1 for any point within the solution domain. The higher the function value the greater the mesh density. The above was interfaced to software requiring only an outline of the full solution domain and a minimum quality of triangle to be supplied.

 

METHOD OF SOLUTION

The numerical algorithm presented here is centred around SPRINT2D [3] software - a general package for solution of equations discretised on triangular domain decomposition. Its methods enable accurate solutions to be found for both smooth and discontinuous flow problems. It allows flexibility of definition of problems and boundary conditions by passing control of the flux calculation, initial conditions, source term and boundary specification to external routines while performing the integration routine within.

 

A key feature of SPRINT2D is its ability to quantify the solution error and to control it as necessary. It contains methods that control temporal and spatial errors inter-dependently - while the mesh is allowed to refine/de-refine, the time step is governed by the temporal integration procedure - see refs in [1][3].

 

DISCRETISATION PROCEDURE

A cell centred finite volume method is formulated for equations (1) over a triangular shaped control volume with the dependent variables of the system represented as piecewise constants. The association of these variables with particular points enables the use of a high-order interpolation scheme. Integrating equation (1) over the i th triangle and using a one-point quadrature rule and the divergence theorem , leads to the semi-discrete equation

 

Equation 2

 

where A_i is the area of the triangle and subscripts indicate sides of adjacent triangles, for example, refering to figure 2, side ij is common to the triangles associated with U_i and U_j and F_n(U) is the outward flux vector.

 

NORMAL FLUX CALCULATION: THE RIEMANN APPROACH

Evaluation of the normal flux in equation (2) is made by a series of solutions local to the lines which make up the triangular mesh. The Riemann problem is defined by the solutions on the left and right of the cell face (or internal and external to the finite volume), and the order of the numerical scheme is determined by the definition of these two data states.

 

 

Figure 2: Spatial discretisation procedure

 

If is a local orthogonal coordinate system centred at the mid-point of a cell face, with in the normal outward direction equation (1) can be transformed to this local coordinate system as the one-dimensional problem

 

 

This can be viewed as a Riemann problem and solved using an approximate Riemann solver and the solution transformed back again to give the desired flux F_n(U). Many Riemann solvers exist for application to the shallow water equations, Toro [4] identifies several, including an exact solution. In this work the solver due to Roe used.

 

SPATIAL INTEGRATION

For a first order scheme the values of the variables to the left and right of side ij are , , the solution at the centres of the cells. By employing an upwind weighted linear interpolation function second order accuracy can be achieved [3]. Limiting the gradients eliminates the possibility of under and over shoots that second order schemes can exhibit in regions of steep gradients. The limited values on the left and right of face ij, see figure (2), are given by

 

 

Here is a suitable limiter, such as the modified Van Leer, see [3], and the r's are the internal and external upwind bias ratios of gradients and and are internal and external linear upwind values obtained from interpolation between the quantities which can be most clearly seen by examining figure (2).

 

TIME INTEGRATION

Equation (2) may be rewritten as a system of ordinary differential equations

 

 

and integrated numerically to give an approximation to the exact solution values. The solution of this system of equations constitutes the major computational task. The difference between these exact and approximate solution is the global time integration error due to the numerical integration method applied. SPRINT2D contains a number of time integration modules- both explicit and implicit. The Theta method, used in this work, has automatic local error control, which when applied with local mesh refinement/de-refinement enables a balance between spatial and temporal errors to be found, allowing optimum accuracy [3].

 

BOUNDARY CONDITIONS

Boundary conditions are applied in similar way to the flux calculation, via a Riemann solver. A combination of the Riemann invariant and an equation for the particular boundary provides sufficient information for an explicit solution to be calculated. Further details can be found in [1]

 

RESULTS: DAM-BREAK WAVE MODELLING

The method has been applied to the problem of modelling of dam-break waves. These flows typically posses regions of highly complex transitional flow. The particular problem investigated was that posed by the European Union group know as CADAM (Concerted Action on Dam-Break Modelling). The test consisted of a physical model that combined a square-shaped upstream reservoir and a channel with a 45° bend. The flow was essentially two-dimensional in the reservoir and in the bend region joining the two channel reaches. Experimental tests were carried out by the CADAM members the data obtained has been used for comparison with results here.

 

 

Figure 3: Channel Dimensions and Gauge Locations

 

Details of the channel dimensions can be seen in figure 3. The channel had a slope of zero and experimentally measured Manning coefficients of 0.0095 s/m^1/3 for the bed and 0.0195 s/m^1/3 for the wall. A guillotine-type gate simulated the dam and connected the channel the reservoir. The downstream end of the channel was a free outflow. Depth/pressure gauges were positioned at the locations shown.

 

The reservoir initially contained water with the free surface 0.25m above the channel bed level, i.e. the water depth in the reservoir is 0.58 m. The initial water depth in the channel was zero. To simulate an instantaneous dam-break the gate was raised (by hand) as quickly as possible.

 

NUMERICAL SOLUTIONS TO THE 45° CHANNEL

Solutions were computed using several different meshes. A 2008 uniformly distributed triangle mesh, was fine enough to give good resolution. To investigate the flow in more detail at the bend the mesh shown in figure 4 consisting of 2004 triangles weighted near to the bend was used.

 

 

Figure 4: The computational mesh of 2004 triangles with weighted distribution

 

Figure 5 shows detail of the flow solution near the bend with the above 2004 triangle weighted mesh. On the left are shown contours of the Froude number, the lighter indicating high values greater than 1 and the dark low, subcritical, values. To the right are shown the velocity vectors, note how the flow changes direction rapidly and how the regions of slow flow correspond to red regions of the left figure. The two figures very clearly show several transitional shocks. These are very sharply defined due to

 

 

Figure 5: Froude number and velocity vectors at the bend

 

the combination of Riemann solvers and high order spatial discretisation used in this algorithm.

 

COMPARISON OF NUMERICAL AND EXPERIMENTAL RESULTS

Experimental measurements of depth were made at all gauge points and compared to the results from numerical simulation. Figure 6 shows results from three gauges - G4, G6 and G8. Apart from the 'noise' of the experimental measurements the numerical solution matches very well at all gauges. The numerical peak depths are very slightly lower at G4 and higher at G6 and the reflected wave shown by the second depth increase is slightly late and 0.02m too low. This difference is probably due to the highly 3-D nature of the reflected wave. The time of arrival of the main wave at all gauges is almost exact to the second.

 

 

Figure 6: Depth histories at gauges G4, G6 and G8

 

CONCLUSION

A solution algorithm has been described that is capable of modelling highly transitional two -dimensional open channel flows. It is capable of high resolution due to its use of Riemann solvers and both second order spatial and temporal discretisation. The algorithm has been applied in a computer simulation of dam-break waves and shown to compare well to experimental measurements.

 

REFERENCES

[1] Sleigh PA, Gaskell PH, Berzins M and Wright N. An Unstructured Finite-Volume Algorithm for Predicting Flow in Rivers and Estuaries. Computers and Fluids Vol 27. No 4. pp479-508, 1998.

[2] Joe B and Simpson RB. Triangular meshes for regions of complicated shape. International Journal for Numerical Methods in Engineering. Vol 23, pp 987-997, 1991.

[3] Berzins M. temporal error control in the method of lines for convection dominated equations. SIAM Journal of Scientific Computing, May 1995.

[4] Toro EF. Riemann problems and the WAF method for solving the 2-dimensional shallow-water equations. Philosophical Transactions of the Royal Society of London, Series A, Physical Sciences and Engineering 1649, 1992.