A MODEL EQUATION FOR SHALLOW WATER FLOWS

 

GUILLERMO HAUKE

 

Correspondence address: Departamento de Mecánica de Fluidos,

Centro Politécnico Superior, C/María de Luna 3,

50015 Zaragoza, Spain.

Tel.: +34 976 761881, Fax.: +34 976 761882,

email: ghauke@ideafix.cps.unizar.es

 

EUGENIO OÑATE

Centro Internacional de Métodos Numéricos en Ingeniería

Universidad Politécnica de Cataluña

 

 

ABSTRACT

The shallow water equations present important source terms which can become singular as the water level tends to zero. Therefore, as the source terms grow larger and larger, they pose a challenge to numerical methods, specially to finite element methods. In order to study the numerical problems associated to source terms and improve upon finite element methods, a simple model equation is presented. As a result an improved finite element method is proposed which, while simple, behaves correctly with large source terms.

 

INTRODUCTION

In the source terms of the shallow water equations, it is customary to include many transport phenomena, such as bed slope, laminar/turbulent friction, wind shear forces, Coriolis forces, infiltration laws, etc. Many of these terms are singular, and therefore pose a challenge to the stability of numerical methods. In particular, as it will be made apparent below, finite element methods (e.g. Hauke 1998) are at a disadvantage with respect to finite difference methods.

 

Therefore, a one-dimensional model equation is proposed to analyze the behavior of several finite difference methods and finite element methods under large source terms. The cause of the lack of stability will be determined and a simple finite element method will be proposed which has good stability and is free from oscillations for large source terms.

 

THE MODEL EQUATION

The shallow water equations form a nonlinear hyperbolic system of equations with source terms. Hyperbolic equations represent wave phenomena, which can be modeled with advection, , where is the wave speed. The source terms can be modeled according to , where c>0 for production and c<0 for dissipation or destruction. Therefore, an appropriate one-dimensional model for steady shallow water flows can be given as follows. Find such that:

 

(1)

 

The exact solution for (1) can be readily obtained as

 

(2)

 

Equation (1) can also be solved numerically. For that purpose, the domain is discretized with N+1 nodes of coordinates x_i, i=0,....N, which define N equally sized elements of area and length h=1/N. The discrete solution consists of nodal values of the searched function, , which approximate the exact solution at the node x_i, i.e., .

Numerical methods applied to (1) give rise to the corresponding stencil,

 

(3)

 

where the parameters s_i take on different values according to the numerical method employed. A solution of (3) can be found by assuming that the discrete solution obeys the recursive law,

 

(4)

 

where is the discrete amplification factor and represents the growth of the solution from one node to the next. Operating, it can be shown that there are two amplification factors,

 

(5)

 

From (2), the exact amplification factor can be seen to be

 

(6)

 

where

 

(7)

 

is an dimensionless number, ratio of source to advection.

If the discrete amplification factor equals the exact one, the numerical solution will be nodally exact. Note that since there are two discrete amplification factors, one of them can be spurious.

 

FINITE DIFFERENCE METHODS

Two methods will be studied here: the central difference scheme and the upwind scheme for positive wave speed . The stencils and amplification factors are tabulated in Table 1. Figure 1 shows the amplification factors and and a function of .

 

Table 1. Stencil and amplification factors for finite difference schemes.

 

 

As can be seen in Figure 1, the upwind scheme has only one amplification factor which behaves correctly for , that is,

 

	(dissipation)
	(8)
	(production)

 

As for the central finite difference scheme, there are two amplification factors, from which is a spurious root. The root behaves correctly.

Therefore, both methods should perform adequately, oscillation-free, with source terms.

 

 

Figure 1. Amplification factors for the central difference and upwind schemes.

 

FINITE ELEMENT METHODS

Finite element methods for (1) can be formulated as follows. Given the usual solution finite element space S and weighting space V, find such that for all :

 

(9)

 

The first integral is the Galerkin contribution. The second integral is the stabilizing term, where L denotes a differential operator which is applied to the weighting function w.

This weak form encompasses various finite element methods. Different definitions of the stabilizing operator L give rise to different stabilized methods, in particular SUPG (Brooks & Hughes 1982), Galerkin/least-squares-GLS (Hughes, Franca & Hulbert 1985) and the variational multiscale method-VMS (Hughes 1995). The Galerkin finite element method is recuperated by taking .

 

Hereto forth, it will assumed that the finite element spaces are formed by linear elements. Also, expressions will be written as a function of the dimensionless parameter

 

(10)

 

Table 2 shows the stencil for the methods analyzed and Figure 2 displays the amplification factors as a function of . These have been calculated with , the correct value for the advective limit in the advective-diffusive equation.

 

Table 2. Stabilizing operator and stencil for serveral finite element methods.

 

 

Note that for all finite element methods, in the regime of dissipation for , the amplification factor is negative, . Therefore, the numerical solution presents spurious oscillations, even for the Galerkin method. Recall that this problem is not present with finite difference methods, and therefore, the cause for this unstable behavior can be traced back to averaging the source term within an element.

 

STABILIZED FINITE ELEMENT METHOD

Following the work in Onate 1998, the following time scale parameter can be approximately derived:

 

 

With this, the amplification factors for the SUPG method are displayed in Figure 3. Now, the numerical solution is converging with x free from oscillations because for . Even in the absence of nodal exactness, we do obtain good stability.

 

 

Figure 2. Amplification factors for several finite element methods.

 

CONCLUSION

An analysis of the model equation revealed that usual finite difference methods treat correctly source terms. On the other side, to avoid stability problems, finite element methods require more care in the presence of large source terms. The lack of stability can be traced back to averaging the source term within each element. To improve the stability, a simple stabilized finite element method is proposed which treats correctly source terms of any size.

 

 

Figure 3. Amplification factors for the proposed stabilized methods

 

REFERENCES

Brooks, A.N. & T.J.R. Hughes (1982). "Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations", Computer Methods in Applied Mechanics and Engineering, 32, 199-259.

 

Hauke, G. 1998. "A symmetric formulation for computing transient shallow water flows", Computer Methods in Applied Mechanics and Engineering, to appear. 163, 111-122.

 

Hughes, T.J.R. 1995. "Multiscale phenomena: Green's functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins Computer Methods in Applied Mechanics and Engineering, 127, 387-401.

 

Hughes, T.J.R., L.P. Franca & G.M. Hulbert 1989. "A new finite element formulation for computational for computational fluid dynamics: VIII. The Galerkin/least-squares method for advective-diffusive equations", Computer Methods in Applied Mechanics and Engineering, 73, 173-189.

 

Oqate, E. 1998. "Derivation of stabilized equations for numerical solution of advective-diffusive transport in fluid flow problems", Computer Methods in Applied Mechanics and Engineering, 151, 233-265.