F.L.B. Ribeiro,' A.C. Galeao, L. Landau* Programa de ......^Laboratorio Nacional de Computaqao Cientifica, Rua Lauro Muller 455, Rio de Janeiro, RJ 22 290- 160, Brasil Abstract This

A space-time finite element formulation for the

shallow water equations

F.L.B. Ribeiro,' A.C. Galeao," L. Landau*

"Programa de Engenharia Civil, COPPE / Universidade Federal

do Rio de Janeiro, Caixa Postal 68506, Rio de Janeiro, RJ

^Laboratorio Nacional de Computaqao Cientifica, Rua Lauro

Muller 455, Rio de Janeiro, RJ 22 290- 160, Brasil

Abstract

This paper presents a space-time finite element formulation for problemsgoverned by the shallow water equations. A constant time-discontinuousapproximation is adopted, while linear three node triangles are used for thespatial discretization. The streamline upwind Petrov-Galerkin (SUPG) methodis applied in its equivalent variational form to fit the time discretization. Also,the correspondent semi-discrete SUPG version is established, and somenumerical results are presented in order to compare the performance of thesemethods.

1 Introduction

As it is well known, the use of the classical Galerkin method to approximateconvection-dominated phenomena leads to numerical solutions contaminated byspurious oscillations that are spread over the entire computational domain. Inthe context of weighted residual methods, a remarkable improvement on thenumerical solution of such problems was provided by the consistent variationalSUPG method proposed in by Brooks and Hughesfl]. Since then, many SUPGbased methods have been used for multi variable systems of equations.. InSharkib[2] and Almeida and Galeao[3], space-time Petrov-Galerkin (STPG)formulations were derived for the compressible Euler and Navier-Stokesequations, showing their inherent control of derivatives along thecharacteristics. This fact turns out to be an improvement over the classical semi-discrete SUPG formulation, where only streamline derivatives are controlled.The discontinuity capturing approach proposed in those references lead to

Transactions on Ecology and the Environment vol 10, © 1996 WIT Press, www.witpress.com, ISSN 1743-3541

404 Computer Techniques in Environmental Studies

stable and accurate methods to solve all details of sharp layers and/or shockdiscontinuities. If we realize that the mathematical structure of the abovementioned equations is identical to those that govern the shallow waterproblem, it can be immediately concluded that for such problems, P-G weightedresidual methods will also perform well. This was done in Bova and Carey[4]and Saleri[5], where semi-discrete P-G (SDPG) finite elements were employed,and in Carbonel, Galeao and Loula[6], where a space-time P-G formulation wasderived, using linear interpolation for both spatial and temporal discretization.In this paper, the continuous linear interpolation is retained for the spatialdiscretization, but constant time-discontinuous interpolation is adopted. For thischoice of interpolation, the weighting function gives no contribution to theterms involving time derivatives. Even so, the resulting discrete space-timeGalerkin (STG) equations coincide with the semi-discrete Galerkinapproximation with an Euler backward difference scheme to approximate thefirst time derivatives. Nevertheless, this will not occur with the P-G weightingresidual methods. The numerical examples that will be presented later will showthat, under these circumstances, the STPG formulation will give more accurateresults than the SDPG model.

2 Problem statement

Let (x, y) G Q e W* define a set of points on an horizontal plane and let z e [-

A, 77] denote the vertical direction, where h(x> y) represents the water depth andrj(x,y,z) is the water surface elevation, both measured from the undisturbed

water surface. We start from the 3-D incompressible Navier-Stokes equations,after turbulent time-averaging, integrating these equations along the z directionusing depth-averaged horizontal velocities. Under the simplifying assumption ofa hydrostatic pressure distribution (negligible vertical acceleration), we arrive atthe shallow water equations:

1 * - lKA + (qw|vt' + g|u|fO+ //(*/,*,+f ) (la)

- Ww

In these equations, H = A» 77 is the total water depth and u is the depth-averaged velocity, with components w and v in x and y directions respectively.The gravitational acceleration is given by g and/is the Coriolis parameter. Thewind velocity is w, with components v/ and w/ a and C are, respectively, thesurface and Chezy friction coefficients, and // is the eddy viscosity.


Computer Techniques in Environmental Studies 405

Multiplying the third equation by g and observing that,

- cClc} O\- K ;>, W

where c = (gH) , and considering similar expressions for (g//),* and (gff),y , weobtain the shallow water equations in the velocity-celerity variables (see Saleri[5] ) which, in matrix form, can be written as:

U,,+A.VU = F ; A.VU = A*VU

u

0

c

0

w

0

c

0 , /i2 0

0

0

c

o"

c

(3 a)

(3b)

(3c)

—

C(3d)

Once an initial state Ua(x, y) is specified at / = 0 and appropriate boundaryconditions are prescribed, the system of equations above can be solved to givethe unknown column vector U.To obtain the space-time description of (3a-d) we introduce the variable change:s = (I)/, where (1) has units of velocity. Then if we define:

(4a)

where I is the (3x3) identity matrix, we will say that the space-time solution ofthe original problem (3a-d) is the (3x1) column vector U that satisfies thetransformed equation:

in (4b)



3 Petrov-Galerkin finite element model

In order to construct the space-time finite element subspace, let us considerpartitions 0 = < t\ <...(„< *„+, of 5R* and denote by /„ = ((„ , ;„+,) the n* timeinterval. For each n the space-time integration domain is the "slab" S» = Q x /„,

with boundary F = DC/,, . If we define S* as the e'** element in £„, e = 1, ...

n , where (TVy „ is the total number of elements in 5%, then for n = 0, 1,2, ...

(/) the space-time finite element partition H^ is such that:

,= '; ,%=0,,x/,,; Q=UQ,;Q,nQ.=0 ybrz#y (5a)g=l e=l

(ii) the space-time finite element subspace consists of continuous piecewisepolynomials on the slab $„, and may be discontinuous in time across thetime levels t^ that is:

tf* EC"(S.); U" ep*(5;); f/"|-=O (5b)

where P^ is the set of polynomials of degree less than or equal to k.

According to the above definitions, the variational STPG formulation for theproblem (4a-b) reads:

Find t/* e »* such that for n = 0, 1 ,2,. . .

W. _ _X J(T . V L/" ). A*d&dt + (6a)

where,

= A.VU* - F" = Uj+AW*).VU'' - F* (6c)



• Remarks

(1) If in (6a) the integrals are taken over Q and fig instead of & and 5^,

respectively; the P-G weighting function (rZ.Vf)'') is replaced by

(rA.VU**)', and in (5b) the finite element subspace is defined for all /

(making the jumping term in (6a) to disappear); and finally, if weapproximate the partial time-derivative by a time-differencing operator, wegenerate the semi-discrete SUPG method.

(2) If in (6a) we do not consider the added SUPG contribution represented bythe second term under the summation symbol, we reproduce the time-discontinuous Galerkin method. If in this case, for instance, we use constanttime-discontinuous interpolation, the resulting space-time finite elementmethod will be identical to the backward Euler semi-discrete finite element

procedure. Since constant time interpolation is used, C/,J' = 0, and therefore

it is clear that the jumping condition, represented by the third integral in(6a), is the term responsible for this equivalence. Nevertheless, the STPGmethod and the correspondent SDPG method will be different.

(3) The definition of the matrix T will be also different in these twoformulations. This point is focused in the next section.

4 Purely hyperbolic problems

To simplify our analysis let us assume F = 0. Using the T matrix definition foundin Sharkib[2], we have:

• for the space-time formulation,

_ s ^ -1-1/2tfVaf 1

(7a)

where x<, = t : *i = x ; X2 = y ; & (k = 0,1,2) are the local coordinates of theparent element S^ ; and A\, AI are the Jacobian flux matrices defined in (3c).

• for the semi-discrete formulation,

- -1/2

(7b)



where £/ (/ = 1,2) are the local coordinates of the parent element i% If, onceagain, for the sake of simplicity , we consider unidimensional problems and

11*' partitions of equal (hexAt) S* elements, (7a-b) simplify to:

(8a)

(8b)

where,

4 4(9a)

(9b)

are the eigenvalues and correspondent eigenvetors of (z

With these definitions we have,

)" and /4, respectively.

*. 17 - | 2

11 + C

«-c

+

-

u-c

u + c

H-Cl-lU

U 4- C \ + \ U - C

c (lOb)

Now let us introduce the non-dimensional factor a = -~- , or what is the same,

CFL = — =>a>\> where CFL denotes the well known Courant-Friedrichs-a

Levy number. Using this factor, the definition of T*t becomes:



A28

(lla)

(lib)

»Remarks

=aV (lie)

(i) Notice the intrinsic dependence of r,, on the used time step A/, which is notaccounted for in %</. In the limit, as At -» 0 ; T,, -» 0, in a consistent way.This does not occur with the semi-discrete formulation because r^ isindependent on At.

(ii) Although not realizable from the practical point of view, %,= r^ if and onlyif a = 0. Remind that a must be greater than one, in order to the CFLcondition be attained. We restate this comment saying that r,, approaches %</as the time step At becomes larger, or, in other words, when accuracydecreases.

(in) If we assume that u (( c, then r,, and r«/ become almost diagonal matrices,

and can be replaced by,

I * Ii; r,U=-—i

Because a > 1, the above ratio between brackets is always less than one.Even in the most unfavorable situation, CFL = 1 = a, r>, introduces lessdissipation than r,/.

In order to get a deeper insight about the performance of the STPG andSDPG methods, some numerical experiments will be performed in the nextsection. For these examples, constant discontinuous time interpolation andpiecewise linear continuous spatial interpolation will be adopted for the STPGmethod. For the correspondent SDPG formulation this same spatialdiscretization will be used, combined, with an implicit backward finite differencescheme for time discretization.



5 Numerical results

Our first example is the well known dam break problem, which consists of awall separating two undisturbed water levels that is suddenly removed (Figurela). Friction effects are neglected and the spatial discretization is given by a2x100 triangular elements mesh, as illustrated in Figure Ib. Figures 2-3 showthe results for / = 2.50, respectively, comparing the solutions obtained with theGalerkin, the space-time and the semi-discrete formulations. For a time stepAf = 0.10 (Figure 2), Galerkin solutions exhibit oscillations in the entirecomputational domain. This does not occur for both, the semi-discrete and thespace-time P-G solutions, which accurately approximate the high gradientsbetween the three horizontal water levels. For this time step, the semi-discretesolution is sharper than that obtained with the space-time formulation. Theeffect of reducing the time-step is shown in Figure 3, where the resultscorresponding to Af = 0.05 are plotted for / = 2.50. For the Galerkin solutions,the oscillations grow up. For the semi-discrete solutions some localizedoscillations appear near the sharpest layer, while a sharper solution withoutoscillations is obtained with the space-time method.

(a)

h= 1

g = 10L= 100

L/24

L/2

L/ 100

I Jl ~~ ML/ 100

Figure 1: Dam break problem.



2.20-

2.00

1.80

1.60

1.40

1.20

1.00

0.80-

= 2.50

GalerkinSpace-timeSemi-discrete

20 40 60 80 100

Figure 2: Solution for time / = 2.50, A/= 0.1.

2.20-

2.00

1.80

1.60

1.40

1.20

1.00

0.80-

= 2.50


20 40 60 80 100

Figure 3: Solution for time / = 2.50, A/ = 0.05.



The second example, illustrated by Figure 4a, is the problem of areflecting wave in a frictionless horizontal channel of length L = 5000,discretized with 2x10 elements, as shown in Figure 4b. The channel is open atthe inflow boundary and closed at the opposed boundary. The system issubjected to a boundary condition at point A, raising the water level suddenlyfrom the initial state of rest (H = 10) to H = 10.1, within one time step. Theresults can be seen in Figures 5-6. In these figures, the time-history responsesfor the water surface elevation at point B are depicted. For A/ = 10 (Figure 5),the curves of both, Galerkin and semi-discrete solutions, almost coincide. TheSTPG method presents an overdiffusive behavior for this time step, leading to asolution that progressively damps along time. However, a completely differentbehavior occurs when A/ = 1. For this time step, the STPG solution reaches therectangular pulse form, while the Galerkin and SDPG solutions present someoscillations, as can be observed in Figure 6.

(a) h= 10

L = 5000

(b) ~>-J ___ __zz L__JZz J___L :— __ _Lc=rc2 Lf L/25

Figure 4: Reflecting wave in a frictionless channel.



0.25-

0.15

0.05

-0.05 J

At = 10

6000

Galerkin-Space-timeSemi-discrete

8000

Figure 5: Solution at point B, with A/ = 10.

0.10

0.05-

o.oo

-0.05


6000 8000

Figure 6: Solution at point B, with A/ = I.



6 Conclusions

In this paper, a STPG finite element model was derived for problems governedby the shallow water equations. A piecewise linear continuous interpolation (inthe space variables) was used, and piecewise constant discontinuous functionswere adopted for time interpolation. In addition, the correspondent semi-discrete P-G model was also presented. The numerical results showed in thiswork indicate that the STPG method performs better than the SDPGformulation. Finally we should comment about the necessity of using aconsistent time-space shock-capturing PG weighting function such that, in thelimit, when Af-> 0, the approximated numerical characteristics approaches thetrue characteristics. This can not be actually provided by the SUPG method.The CAU generalized method proposed in Almeida and Galeao [3] forcompressible flows, fulfill these requirements and, can therefore be successfullyapplied to shallow water problems.

References

1. Brooks, A. N., Hughes, T. J., Streamline Upwind Petrov-GalerkinFormulation for Convection-Dominated Flows with Particular Emphasis onthe Incompressible Navier-Stokes Equations, Compiit. Meth. Appl. Mech.Engrg, Vol. 32, pp 199-259, 1982.

2. Shakib, F, Finite Element Analysis of the Compressible Euler and Navier-Stokes Equations, Ph.D. Thesis, Stanford University, 1988.

3. Almeida, R. C , Galeao, A. C , An Adaptive Petrov-Galerkin Formulationfor the Compressible Euler and Navier-Stokes Equations. Comput. Meth.Appl Mech. Engrg, Vol. 129, pp 157-176, 1996.

4. Bova, S. W., Carey, G. F., An entropy Variable Formulation and Petrov-Galerkin Methods for the Shallow Water Equations, in: Finite ElementModeling of Environmental Problems-Surface and Subsurface Flow andTransport, ed. G Carey, John Wiley, London, England, 1995.

5. Saleri, F., Some Stabilization Techniques in Computational Fluid Dynamics,Proceedings of the 9"' International Conference on Finite Elements inFluids, Venezia, 1995.

6. Carbonel, C , Galeao, A. C , Loula, A. D., A Two-dimensional FiniteElement Model for Shallow Water Waves, Proceedings of the 73"' BrazilianCongress and 2"** Iberian American Congress of Mechanical Engineering,Belo Horizonte, 1995.


Documents

F.L.B. Ribeiro,' A.C. Galeao, L. Landau* Programa de ......^Laboratorio Nacional de Computaqao Cientifica, Rua Lauro Muller 455, Rio de Janeiro, RJ 22 290- 160, Brasil Abstract This