An efficient shock capturing scheme for two-dimensional, open channel, unsteady flows in a generalised coordinate system

Computers Math. Applic. Vol. 29, No. 1, pp. 83-88, 1995

Pergamon

Copyright©1995 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0898.1221(94)00208-8

0898-1221/95 $7.00 + 0.00

An Efficient Shock Capturing Scheme for Two-Dimensional, Open Channel, Unsteady Flows in a Generalised Coordinate System P. GLAISTER Department of Mathematics, University of Reading P.O. Box 220, Reading R G 6 2AX, United Kingdom

(Received and accepted February 1994) A b s t r a c t - - A n efficientalgorithm based on flux differencesplittingis presented for the solution of two-dimensional, open channel flows. A transformation maps a nonrectangular, physical domain into a rectangular one. The governing equations are then the shallow water equations, including terms of slope and friction,in a generalised coordinate system. A regular mesh on a rectangular computational domain can then be employed. The resulting scheme is efficient,has good jump capturing properties and the advantage of using boundary/body-fitted meshes. Keywords--Shock

capturing, Shallow water equations, Generalized coordinates.

1. I N T R O D U C T I O N In a recent paper [1], an efficient numerical scheme was proposed for the solution of the twodimensional shallow water equations in a rectangular domain. For two-dimensional, open channel flows in complex geometries, it is convenient to make predictions using a boundary/body-fitted, computational mesh. The purpose of this paper is two-fold. First, it is to extend the scheme in [1] to apply to such flows. A transformation maps a nonrectangular, physical domain into a rectangular one, on which a regular mesh is placed. The resulting equations apply in a generalised coordinate system and include the terms of nonzero bottom slope and friction. The scheme incorporates a treatment of these terms. Second, it is intended to generalise this work, which involves the use of the arithmetic mean for the average of the flow variables, to the Euler equations, and this task was stated in [1]. Thus, this paper is also to be seen as providing an important intermediate link in this development. Time-marching solutions of the governing flow equations are frequently used for open channel flow analysis. Any numerical scheme must have the ability to compute mixed subcriticalsupercritical flows with automatic capturing of hydraulic jumps, and can be utilised to eliminate the most common cause of spillway failures, namely the improper design of steep shutes. Several authors [2-9] have proposed schemes for such flows; however, there remains the difficulty of capturing of bores. The scheme presented here has good jump-capturing properties by design, and this was successfully demonstrated in [1] for the special case of a rectangular physical domain. 2. G O V E R N I N G

EQUATIONS

The two-dimensional St. Venant equations governing open channel flows can be written in conservation form as Typeset by AA4S-TEX 83

84

P. GLAISTER

wt + F x + G z = f + g ,

(2.1)

where W

-~- (•, ¢?~, CW) T,

F(w) =

(2.2a)

Cu, Cu2 + T ' Cuw

,

(2.2b)

G(w) = (¢w, ¢uw, ¢w2 + ¢---~],

(2.2c)

f(w) = (0, g¢ (ht - st), 0) T

(2.3a)

g(w) = (0, 0, g¢ (hz - sz)) T ,

(2.3b)

and

¢ = g (7 + h).

(2.a)

The quantities ¢ = ¢(x, z, t), and u = u(x, z, t) and w = w(x, z, t) represent g multiplied by the total height above the bottom of the channel, and the components of the fluid velocity in the x and z directions, respectively, at a general position x, z and at time t. The gravitational constant is represented by g and the undisturbed depth of the water is given by h(x, z). The elevation 7 = 7( x, z,t) above the plane y = 0 is measured in the vertical y direction. The quantities sx, Sz are the slopes of the energy grade lines in the x-, z-directions, respectively, and are determined from the steady-state friction formulae (in SI units) sz=

n 2 . ~

(¢/g)4/3

,

Sz=

n2w v ~ r T ~ (¢/g)4/3

'

(2.5a, b)

where n represents Manning's roughness coefficient. Alternatively, the following form could be used:

sx =

¢ C2

,

Sz =

¢ C2

,

(2.6a,b)

where C is Chdzy's flow friction coefficient. If we introduce a transformation x = x(~, 7), z = z(~, y) from the physical (x, z) space to a computational (~, y) space, equations (2.1) become ( J w ) t + F¢ + ( ~ = f + ~,

(2.7)

where

w = (¢ ¢u, ¢w) v,

(2.8a)

~'(w) =

Cu, 51 z~¢ 2 + CuU, - ~1x ~ ¢2 + ¢wU

(2.8b)

~2(w) =

1 x¢ ¢2 + 4noW CW, - ~1 z¢ ¢2 + Cuw, ~

(2.8c)

f(w) = (0, g¢ (z~ he - s x~ W ) , g ¢ ( - x ~ he - s z~ W)) T ,

(2.9a)

$(w) = (0, g ¢ ( - z ¢ h, - sz¢ V), g¢(z¢ h , - s z ¢ ~]))T,

(2.9b)

and U=z~u-xnw,

W=-zcu+xcw,

(2.10a,b)

An Efficient Shock Capturing Scheme

where

n2 ~ s=

85

g or

(¢/g)4/3

¢C 2

,

(2.11)

depending on whether (2.5a, b) or (2.6a,b) has been employed. The terms f, ~ are associated with derivatives in the (, ~? directions, respectively. The Jacobian of the grid transformation is given by J = x( Y,7 xn Y¢. (2.12) -

NOTE. We also have

1

(¢~, ¢~, n~, ~ ) = ~ (y,, - x , , -y¢, x¢). 3. O P E R A T O R

SPLITTING

We propose solving equation (2.7) using similar ideas to those given in [1], in particular, the technique of operator splitting [10], i.e., solve successively wt + F ¢ -- f,

and

(3.1a)

wt + (~, = ~,

(3.1b)

along ~- and y-coordinate lines, respectively. We give the scheme for the solution of equation (3.1a), and the scheme for the solution of equation (3.1b) will then follow by symmetry.

4. N U M E R I C A L

SCHEME

Following the approach outlined in [1] for the shallow water equations in Cartesian coordinates, a first order, explicit numerical scheme for equation (3.1a) along a ~-coordinate line 77 = ~0 can be written be written as as At %vre'b1 ~ ~i ~i e~, (4.1) k = w~ Jj-(1/~) A¢ i=1 where At, A~ represent mesh spacings in the t and ~ directions, respectively, w~n represents the numerical approximation to the solution w(~j, W0, m At), after m time steps, and at a point ~j on this coordinate line. The quantity J j - 1 / 2 represents an approximation to the grid Jacobian (2.12) in ( ( j - l , ~j), which can be given analytically, or constructed numerically, once the grid has been generated [11]. The value of k could depend on the sign of each ~i (see below) as the technique of upwinding is employed, i.e., k could be j or j - 1. Thus, the practical implementation of (4.1) is to a d d - J~_l/2a¢ At At )-~ ' -y i e-i to w~"_l, if~i < 0, to obtain w~n+l, or to add ~i 7i ei to w[*, if ~i > 0, to obtain w~n+l, for each computational cell (~j-1, ~j) on the coordinate line y = ~0. The quantities ~i, ~ represent approximations to the eigenvalues and eigenvectors, respectively, of the Jacobian aF/Ow in the cell (~1-1, ~j). These, together with the appropriate wavestrengths, ~ , represent an extension of those contained in [1], and which we now derive.

5. D E R I V A T I O N

OF AVERAGES

We begin by considering the case when the source term f = 0. 5.1. C o n t i n u o u s C a s e

We first note that the Jacobian of the flux function in the continuous case is 0

A = 0~" = z~ (¢ - u 2) + u w z , 0w \ _ x ~ ( ¢ _ w2)_ u w z ~

z, 7 2 u z n - w z,7 wz,

--U X~

u z,~ - 2 w x, 7 /

(5.1)

86

P. GLAISTER

whose eigenvalues, A~, are given by

A~= U+v~,

U,

i = 1,2, 3,

(5.2a-c)

with corresponding eigenvectors

e~,2 --

1, u+ v ~ ' w~: v ~ J

'

and

(5.3a,b)

e3 = (0, x,, z,7)T,

(5.3c)

where U = z, u - x , w,

d=

and

(5.4)

+

(5.5)

5.2. Discrete Case Here we seek to determine an approximate Jacobian .& which is the average of the approximate Jacobian evaluated at w}nl and w ~ , such that (5.6)

A ( W R -- WL) ---~F ( W R ) -- F ( W L ) ,

where we have denoted w~_ 1 = WL and w~n = WR for simplicity, and w and F are given by (2.8a,b). We achieve this by determining matrices I3 and (~ such that

l~(qR -- qL) ----W R -- WL,

(5.7)

and

(~(qR -- QL) = F(WR) -- F(WL),

(5.8)

where q = (¢, u, w) T, so that k = C ]~-I.

(5.9)

First, since A(¢u) ----fi A¢ ÷ ~ Au,

and

(5.10a)

A(¢w) = ~ A¢ ÷ CAw,

(5.10b)

where the overbar denotes the arithmetic mean of left and right states, we have

f3 =

~

.

(5.11)

0 For (5.8), it is necessary to treat terms of the form A(¢u) = A (¢ [z,7u - x v to]) = A(z n Cu) - A ( z ~ Cw). For this, we assume that the grid derivatives z~ and x n are constant in the interval (~j-1, Cj) -- (~L, ~R), and denote these by Z and X, respectively. Thus A ( ¢ u ) = Z A ( ¢ u ) - X ACCw) = Z (~ Z~u + ~ A ¢ ) -- X (~ A w + ~ = ~A(Zu

- X w ) + (Zf~ - X ~ , ) A ¢ = Z ~ A u

---- C A U + U A ¢ ,

~¢)

-- X ~ : , ~ . + (Zr~ - X , r , )

a¢ (5.12)

for example, where 0 = Zfi - XtD.

(5.13)

An Efficient Shock Capturing Scheme

87

Similarly, (5.14a,b) A(¢uU) = A (¢u [z, u - x , w]) = Z A(¢u ~) - X A(¢uw).

(5.15)

For the first expressions, we take A(¢u ~) = ~ A ¢

(5.16)

+ 25~ Au,

where (5.1T)

.~ = ~ (~t + . ~ ) , and for the second a(¢UW) = $ ~(UW) + ~-~ ~ ¢ = ~ W

where

+ ~ U

+ W~a¢,

(5.1s)

1

(5.19)

u"~ = fi (u,~ WL + UR Wa),

so that

a(¢.u) = ( z ~ -

x~)

a,+ (2z~- x~)~a.-

x~a~

= u-UA¢ + ~ 6" a u + ~ AU.

(5.20)

Similarly, A( CwU)

~ ) ,',~ + z $,~ ,',,, + (z ~ -

(z ~ - x

2x ,~) $ a~

= w U A ¢ + ~ O A w + ~,~AU.

(5.21)

Combining equations (5.12), (5.14a,b), (5.20), and (5.21) yields

[ z~-x,~ C= [ z ~ + z - ~ - x ~

z~ 2z~,~-x~,~

-x~ -x(~

] |.

L - x ~ +z w ~ - x - ~ z };~ z ~ - 2x ~ j Combining(5.21) and (5.22), the requiredJacobian is obtainedfrom(5.9) as Z -X ] A.= [ Z (~ - fi~O)+ Xh-'~ 2Z~-X~ -X~ , k - x (¢ - ~2) - z Z ff~ Z ~ - 2X ~

wherewe havesimplifiedexpressions 2 ~ - ~ = fi2, 2 if,2 _ ~-~ = ~b2,

(5.22)

(5.23)

(5.24a,b)

with, say, ~ = v ~ UR,

~ = V ~ WR,

and

(5.25a,b)

1

(5.26)

5.3. Eigenvalues and Eigenvectors The eigenvalues and eigenvectors of the matrix A. can then be shown to be ~=O+v'~,

[

~,,2= 1,~,+ (Z~+~AuAV)

(X~-~AwAV)] •

and (5.28a,b) (5.2Sc)

~3 = (0, x , z ) w,

where

(5.27a-c)

~,

1

a = ~2 ( x 2 + z ~) + ~ (av)2.

(5.29)

88 5.4.

P. GLAISTER

Wavestrengths

Finally, to complete the scheme, it is necessary to give the wavestrengths ~ , satisfying

WR--WL

= E ~ei'

(5.30)

i=1 and these are

"~I,2 ---- ~

A• 4- (~

,

(5.31a,b)

~2 (X Au + Z Aw) '~3 =

d

(5.31C)

For the case when f ~ 0, a projection of this evaluated at average states within the cell can be projected onto the eigenvectors, resulting in modified wavestrengths ~ as outlined in [1].

6. C O N C L U S I O N S A conservative finite difference scheme has been proposed for the solution of two-dimensional o p e n channel flows, based on flux difference splitting. T h e scheme treats flows where friction and b o t t o m slope terms are included, and represents a generalisation of an earlier scheme [1]. I m p o r t a n t l y , the j u m p - c a p t u r i n g properties of the original scheme remain. However, a n y physical d o m a i n can now be treated, including nonrectangular bodies/boundaries. This has been achieved using a grid transformation. T h e resulting scheme is c o m p u t a t i o n a l l y efficient, and the techniques discussed can be employed in future work on the Euler equations.

REFERENCES i. P. Glalster, A n efficientnumerical scheme for the two-dimensional shallow water equations using arithmetic averaging, Computers Math. Applic. 2'7 (1), 97-117 (1994). 2. J.V. Soulis, A numerical method for subcritical and supercritical open channel flow calculations, Int. J. Numer. Meth. Fl. 13, 437-464 (1991). 3. J.A. Liggett and S.U. Vasudev, Slope and friction effects in two-dimensional high speed channel flow, International Association for Hydraulic Research 11th Int. Con9 r. Leningrad (1965).

4. J.J. McGuirk and W. Rodi, A depth-averaged mathematical model for the near-field of side discharges into open channel flow, J. Fluid. Mech. 86, 761-781 (1978). 5. R.S. Chapman and C.Y. Kuo, Application of high accuracy finite-difference technique to steady free surface flow problems, Int. J. Numer. Meth. Fl. 3, 583-590 (1983). 6. R.S. Chapman and C.Y. Kuo, Application of the two equation k - e turbulence model to a two-dimensional, steady free-surface flow problem with separation, Int. J. Numer. Meth. Fl. 5, 257-268 (1985). 7. A.O. Demuren, Numerical computation of supercritical and subcritical flows in open channels with varying cross-sections, Proe. 4th Int. Conf. on Finite Elements in Water Resources, Hannover, June 1982, Computational Mechanics Centre, Southampton, 4.13-4.23. 8. J.B. Herbich and P. Walsh, Supercritical flow in rectangular expansions, J. Hydraul. Res. Div., Proc. ASCE, (100), 553--568 (1972). 9. R.J. Fennema and M.H. Chaudry, Implicit methods for two-dimensional unsteady free-surface flows, Jour. Hyd. Res. 27 (3), 321-332 (1989). I0. N.N. Yanenko, The Method of Fractional Steps, Springer-Verlag, Berlin, (1971). 11. J.F. Thompson, Z.U.A. Warsi and C.W. Mastin, Numerical Grid Generation, Foundations and Applications, North-Holland, (1985).