Fox J.A. Transient Flow In Pipes, Open Channels And Sewers

176

Unsteady

flow in

open channels

lch.

7

Fig.7.23

-

Surge initiation.

known.

Then

q(t)

-

u(t)d(t)

assuming that the reach

is

of

rectangular cross-section so

that

q is

the

flow

per

unit

width.

So

g(t):'@ry

but

Sec.7.12l

The

unsteady

travelling

surge

177

For

a specified

value

of

q(t') a

value of

c(r) can

be

obtained

by solving

the

cubic

equation

above.

Then

dq(t)tdt

{6lc(t)12

+ 2(u

s

-

?*dc(t)

|

s)

Having

calculated

c(r)

and

dc(t)tdt

the

value

of x"

can

be calculated

from

the

"qouti6n.

It

is

possible

that

the

value

of x"

so calculated

may

be

greater

than

the

length

of the

channel,

in

which case

no surge

will occur-

ittong

the

r axis,

q rvill be

changing

and

may

be increasing

at an increasing

rate.

In

this

case,

surges

starting

from

iater

times

on

the

i axis

may'oe

initiateriand

will

r-un

after

surges

ittiti"ted

from

earlier

times,

eventually

overtaking

them

and

incorporat-

ing

them

into

themselves.

Once

the surge

is

initiated,

it must

be tracked

and

this

should

be

done

using

the

complete

charaiieristic

equations

so

that

frictional

effects

and the

attenuatuation

that

they

cause,

can

be

taken

into

account.

In stlady

state,

the

surge

is simple

to

calculate.

The

upstream

florv

must

have

a

Froude

number

gSeater

than

unity

and

the downstream

flow must

have

a Froude

number

less

that

o"ity.

When

the

surge

moves,

i.e.

it becomes

a travelling

surge,

it

is

possible that

both

Frbude

numbers

are

less

than

unity

or

they

may

span

unity.

If a case

in

which

the

Froude

numbers

on

either

side

of the

moving

jump

were

both

less than

unity

and

analysed

using

characteristics,

and

ignoring

the

jump,

the

results

would

produce

a

profile in

which

a

smooth

wave-like

surface

will

occur

where

the

jump

is

located.

ThG

is

a consequen@

of

the

diffusive

nature of

characteristics.

However,

this

will

not

work if

the

flow

is such

that

Froude

numbers

span

unity

across

the

jump

unless

the

supercritical

flow

is correctly

analysed

using

the

trro

forward-

stoping

characteristics

ippropriate

to

supercritical

flow.

This

could

be

considered

a

t".itoi!o"

of analysis

bul

theie

is an

accurate

way of

tracking

the surge

and

this

must

be

considered

the

technique

of

choice.

The

possible types

of

zurge

are

shown

in

Fig. 7

.24.T\emathematical

conventions

of sign

will be observed

throughout

the

following

development.

For

instance'

a surge

travJfling

downstream

traveli

at

a

velocity

Dnand

this

is considered

positive;

when

the

wavelravels

upstream,

its

velocity

is

considered

negative.

Th9 flow

velocities

are

also

considered

positive

if

they

are

directed

in the

direction

of x increasing.

The surge

trivels

along

a

path in the

r

versus,

space

TP.

Characteristics

can

be

drawn

from

P: a

positive

characteristic

is

r

SrPz.Point

Thas

twovalues

Tland

T2,a;rrr

surge

located

at these

points is a discontit

the

two

possible cases.

If a travelling

surl

magnitude

u,,

but

of

opposite

sign

to

that

it

So

dc(t)

dt

u(t)

-

Zc(t):

Do-ko .

Therefore

u(t)-2c(t)+us-2c,

and

hence

q(t):QrQ)* uo- zrdry

For

a specified

value ot

q(t)

a

value

of c(t) can be obtained as follows:

q(t)

-?k4)l'*(uo

-?ti

[r(r)]'

.

8'"' I

zone

of

quiet

178

Unsteady

flow

in

open channels

this

case

is impossible

Fig. ?.24

-

Possible types of surge.

upstream

supercritical'flow

(a)

Ax

,l!

Ax

Tt Tz

upstream

subcritical

flow

line of

the

wave

Path

line

downstream

/

subcritical flow

downstream

subcritical flow

surge

wave

lch.

7

1l'2

At

,z/i

A-

Ns2

C

I

(b)

p"r(

Fig.7.25

-

A surge located at

point

T.

Sec.7.12l

The

unsteady

travelling

surge

t79

surge

is

brought

to

rest.

The

velocity

of flow

upstream

of the

now

stationery

wave

will

thei

be ur-u,.,

and

downstream

the

velocity

of

flow

will

be

ur-u,,

and,

using

these

values,

thl

Froude

numbers

relative

to

the surge

can

be calculated-

These

Froude

numbers

Fr1

and

Fr2wlllthen

span

unity

for

all surges'

Fig.7.Zi

shows

ihut

thr

flows

relative

to the

stationary

jump

are

from

relative

Fig.7.26.

supercritical

flow

to

relative

subcritical

flow

based

on

the

relative

velocities

ur-u,,.,

and

,-u*.The

relative

velocity

u2-u,,will

be of

opposite

sign

to

the

velocity

ur'

Note

that,

iti fi!.

7.25(b),the

travellinglurge

has

a smooth,

almost

glassy surface

and

does

not

displiy

the

iough

turbulent

nature

that is

usually

associated

with

a

jump.

Surges

may

travelin

either

direction

but

it should

berealized

that

there

is one

case

which

cannot

happen.

This

forbidden

case

is

illustrated

in

Fig.

7.24.

|t is

a negative

surge

and

cannot

occur

because

it

would

require

a

decrease

in entropy

which

cannot

occur

in

this

universe.

7.12.1

The

analysis

of

the

case

in

Fig.

7.25(a)

point

p,

can

be solved

directly

from

the

two

super

critical

characteristic

equations.

* upr

un,

*

E^,

dt:0

,

_

CSr

+

upr

D5,

*

E5,

dr=

0

The

values

of

.Ap,and

ur, can

be obtained

by solving

these

two

equations

simulta-

neously.

once

A",

has

bden

determined

the

value

of

d",can

be

evaluated.

+v1

180

Unsteady

flow

in

open channels

There

remains the

negative characteristic coming from

.Sr,

tion and

the surge equation.

The

variables zr, and

z

rrare

the depth

of the centroids

of the cross-sections upstream

and

downstream

respectively.

These three equations

must next

be solved

simultaneously.

7

.12.2

The analysis of

the case

in Fig. 7.25(b)

In this

case the two

characteristic

equations downstream of

the surge

give

an

immediate solution

for

Aprand u";:

/A, -\

.

cRz

lt'

-

I

|

* ur,- DRz *

Enzdt:0

,

'\an,

/

-

/r \

-cs,

(*-1f

+u",-us,* Er,dt-0

.

\nSz

/

Upstream

of the surge the

positive

characteristic equation

is available:

/A"

\

.",

|,fr

-

r) *

up,-uR, * E^, dr=0,

Ar,(urr- u*):

Apr(upr-

Dn)

,

t)n

:

u

P,!fut

(A

P'1.P'

-

A''z'')

-

/rr

Arr- Ar,

There

are five equations

the first

two of which

are simultaneous

and can be solved

directly.

The remaining

three remain to

be

solved simultaneously.

By combining

the

first

two of

these

three,

up, can

be

eliminated

from the second rwo.

This

leads to an

expression

for u,.,from both

the second and the third equations and

these two

values

must

be

equal. z is a function

of

A and, as up,andAp, have already been

found,

it only

rernains to solve for Apr. Using

the

value

of.

Arras an approximation

for Ap,

an

iterative

process

can

be

slarted

to solve for Ap,. This is

simple

if the cross-section

is

Sec.

7.121

The unsteady travelling

surge

rectangular but,

when it is

not, difficulties

will be encountered.

A method

based on a

modification of the

Newton-Raphson

method is therefore

to be

preferred.

The

expressions for

the

two values of

u,.,, namely u*, and

uwz,can be

used

to form a

function of the

depth dp"

which must be

zero, i.e.

f(d):

Dn,-

un,

u,.,,

is the above

speed from

the fourth eqn.

and D,,,, is the

wave

speed

from

the fifth

eqn.

u,,,,r must equal

u,,, of

course.

Now

6f(d)

6d

0.0Ld

and

then

d:d'-W

where

dr is

the

earlier

estimate

of d and

d is the improved

estimate

of d.

This method

is

sometimes called

the secant

method.

The

process

should

be iterated

until there

is

an insignificant

difference

between d' and

d.

The choice of

which

of the

rwo signs in

the t

sign

in the

wave

equation

should

be

used is decided according

toFig.

T .27

Fi1.7.27

-

(a)

Use

of

the

positive

sign;

(b)

use the

negative

sign.

Once the

value

of

the unknown

depth

is determined,

the

value of u,., should

be

calculated and this can be used

to calculate

the

location of

the surge

at the end of

the

next Ar

period.

The entire

process

can

be repeated

until

the entire

time

period

of the

simulation time has elapsed.

Normal

analysis

of the

two sections

of

the

reach into

which the surge divides

it can

be

performed so that

the

entire

reach

is analysed.

The

'

[ch.7

the

continuity equa-

181

-

,,,

(*,-

1) + up,

,", * 8",

dr: o

,

Ar,(u"r- u*)

:

Apr(Dpr-

D-)

and

l)n:

Dpr-

(a)

gA",

(Ar,z",-

Ar.zr.\

Ap,

\

Ar,- Ar,

/

1,82

Unsteady

flow in open channels

lch.

7

programming

to do

this

is not so much

difficult as tricky.

This is

primarily due

to

the

need

to

perform all

interpolations

between the

upstream

end of

the

reach

and

the

point T, and the

point

T2

and the

downstream

end

of the

reach.

No interpolation

process

may cross

the surge

because

of

the discontinuity

that it

presents.

The depth

and

velocity

values immediately

upstream

and

downstream

of

the surge

have

little

interest and

are

only

calculated

so

that u,.,,

the surge

velocity, can

be

calculated

accurately.

There

is therefore

no reason

to retain

or store

them

once

they

have

served

their

purpose.

7.I3 APPLICATIONS

OF

TIIE FRICTIONLESS

ANALYSIS

From the equations

for frictionless

analysis

developed

previously in this

chapter

(p.

15a)

it is

possible to solve

almost

any circumstance

in

which the

assumption

of

frictionless or

near-uniform

conditions

is applicable.

7.13.1

Negative

waves

A - - --r!--- :- ^ ---^--^ --,L:^L +L^ l^-+L +^ l^^-^^^- ^-:+ ^rnftraccac in fha

A nggauvc) wave

ls

a

wavc wtllull ualuD(iD

urg LrgPtu

L(,

trgLrf/.lDe

(rD

rL

l,r\.rSrvrrvD

rrr Lrrv

direction

of

r

increasing.

The equation

(p.

157) that

is then

relevant

is

, dx

^,

'

*-3c(r)

*uo-Zcs

.

This equation

is

that

of a

positive

characteristic

line. For

a

given

depth

dthe

value of

c(r) has'a

fixed

value, so

dxldt

is also

defined for this

depth. Integrating

gives

v

+:3fi|,*uo_

z\fgh

.

t-tr

d, defines the

depth

at

time /.

Note that

/, is the time at

x

:0.

us and

co

are the

velocity

and celerity

respectively

in the zone of

quiet.

This

equation

defines

the

wave

profile

and its changes

with time

and

distance.

For example

at

a time

t and for

a

depth

d(t)

the value of

r

appropriate

can

be evaluated.

Alternatively

at a

position

x

the

variation

in d

or c

with time

can

be obtained.

There

is one

situation

in which the

method

is

extremely

useful;

this

is

the dam

break problem.

The solution

obtained

is admittedly only

two

dimensional

whereas

the real

problem is

three

dimensional

but

the result

would apply

to the

rupture

of a

barrage and

the

analysis

applies

to the

flow downstream

of a

bursting

dam.

When

a

dam

collapses,

a negative

wave

is created

and the situation

can

be

regarded

as

frictionless. The

dam can

be

regarded

as a

plate

across

the

channel

and,

to

represent

the collapse of

the

dam,

the

plate

can

be imagined

as being

moved

downstream

faster

than

the

water can follow

it.

To develop

the

theory

of

this, imagine

such a

plate representing

the

dam,

which,

at the moment

of collapse,

starts

to accelerate

downstream

until

it

reaches

a

velocity

u6.

This

velocity can

be increased

to a large

value and so

will represent

the collapse

of

the dam

(Fig. 7.28).

Sec.7.13l

Applications

of the

frictionless

analysis

183

The

characteristics

drawn

in

Fig. 7.28(a)

are all

negative

characteristics

except

the C*

line drawn

back from

A to its intersection

with OF(shown

as a broken

line).

Along

this broken

line,

on*Zca:

Ds*zcs

but

u, is the

velocity

of the

plate and so

cA caneasily

be calculated

given

that

the

plate

velocity is

known.

Line

AA' is

a C- characteristic;

so for

this

line

-Cn

Earlier

in this

chapter

it

was shown

that, for

a

positive

characteristic

d-rldr

-

Do-ko+

3c(r)

-

1.5u(t)

-

0.5u6

* cs.

In an exactly similar

manner it

can

be shown

that,

for a negative

characteristic,

du/dt

:

uo

* Zco- 3c(t)

:1.5u(r)

-

0.5uo

-

cs

from

unt?rn=

uo * ?*s;

so c3

:0.5u6

*

co

-

0.5ur. Substituting

this

result into

dxld/:

na- ce

gves

the

required

result.

However,

the flrst

statement

has u(t)

and

c(r)

as

variables

and

the second

has

u, and ca.

It can therefore

be concluded

that

a

disturbance

may

be initiated

from

any location

in

the flow

field

and need

not

be

specified

only

on

the

f axis.

In Fig. 7

.28(a)

the line OABC

is not

a characteristic

line;

ii

is the

trajectory

of

the

plate in the

.r-t space.

The curved

section

of

it

is

the

acceleration

phase

of its

motion

and

the

straight

portion is

the constant-velocity

section.

Along

line

AA'

,

dxldt:

Do * Zco

-

3ca

from

the

result

just

obtained

and

so.

A'

is

a

straight

line becaus

e dxldt

is constant;

thus

all the other

C

-

lines

are straight

also.

firereiore

u * 2c

is constant

everywhere

and

the slopes

of

the C-

lines coming

off

the

BC

portion

of

the OABC

line

become

the same;

hence

the lines

become

parallel.

Along

all these

line u

-

2c

-

un

-

Zco. The

space

between

Ox

and OABC

can

be

split

into three

zones.

7,onel.

This

is the

zone of

quiet in

which the

velocity is

zero

and the

depth

is

do.

7,onell.

This is

the

zone in

which

C- characteristics

diverge

owing

to

the

accele-

ration

of

the

plate.

7,one

ltr. This

is the

zone

in which

the C-

characteristics

are all

parallel. The

depth

and

velocity

are therefore

constant;

so

all the

water is moving

at

the

velocity of

the

plate.

The

depth of

water in

zone

III adjacent

to the

plate

is

obtained

as follows:

ut*2c6:

Do*Zco

.

This

comes

from

? C

*

characteristic

starting

from OF

up

to B

(not

shown

on

the

diagram);

so

d.r

&:'n

184

Unsteady

flow in open channels

lch.

7

Fig.

7 .?A

-

(a)

the

profile of the moving

dam

(physical plane); (b)

the x

versus r diagram

of the

movingdam.

cs

:0.5uo

-

0.5vu

* co -

Noting

that

us:Da

ard

uo=0

where uo is the

velocity

of

the dam;

then

cr:

-

0.5uo

* c6 and so

)

_r'"

:(co

-

0.5u0)2

ou:T.

g

Now dr

cannot be

less

than

zero; so

ud cannot exceed

2ro.

ff uo

--

ksthen

c

n:

d

a:

0,

i.e. there is

zero depth at the

plate.

Thus,

if the

plate

moves

faster

that Zcs,

the

water

loses contact

with the

plate

and the

water then

feathers

out to

zero depth,

moving

at a velocity ko. The segment

OAB

of

the

line

OABC

can be

shrunk

to a

point

by imposing an

infinite acceleration

on the

plate

up

to

a speed

of.Zcoand

then

zone

III

disappears

(Fig.

7.29). The

field of flow

is contained

in the fan-shaped

area

FOC.The

depth

and velocity at any

point within the

fan

can be

obtained

by using

the equation

dxldr:

uo*Zcs-3cn

which is

the

result for

a C-

line. us

:

0 and

co--

ffiand

so

Sec.7.13l

Applications

of the

frictionless

analysis

185

Fig.7

.29

-

(a)

r

versus

I

plot of complete

dam

removal;

(b) physical

plane.

where

A is

a

point anywhere

in the

flow.

When

dxldt:0,

i.e.

the

vertical

characteristic

through

O,

A

must

be on

this

vertical.

Then

o:2@o-3\M

.

So

Using

the

C*

equation,

us*Zcs:

ue*Zcn

.

Now

co-

tfsd.n:?€O:'rro-

Also

Yor:lrlgdo-

g![sdA

186

Unsteady

flow

in

open

channels

lch.7

Dn:2cs-

2cn

and

so

Thus on

the

vertical

characteristic

at the original

dam site

the

local

velocity equals the

wave

celerity.

7.13.2

Dam break

wave

with a downstream

submerged

river bed

The analysis

already

developed

can

be further

developed

to describe

the

situation of

a

dam

burst

into

a channel

in

which there

already exists

a

depth of water.

The

flow

profile

shown

in Fig.

7.30 cannot occur.

It cannot

happen

because,

at a

Fig.

7.30

-

The

profile

of the flow.

point

B

where

the

downstream

flow

profile

meets

the dam

break

flow already

described,

the

flow

to

the

left of the

point

has a

quite

considerable

velocity but

immediately

to

the

right of

B the

flow

is

stationary.

This is an

impossible situation

because a

change

in momentum

flux,

without

a

force

to cause

it, is a

direct

contravention

of

Newton's

second

law.

To solve

the

surge

the

continuity

equation

and the

momentum

equation

are

available

and together

with

the C*

line to the

point

C

(shown

as

a broken

line)

the

surge can

be solved

(Fig.

7.31).

It is important

to remember

that

the

depth

from

the

point

C

to

the surge

is constant

because

there is no

friction

and so

the

velocity

upstream of the

surge

is the

same

as that

at

C.

For C

to be located

on

the

vertical

through

O

the

depth

do downstream

of the

ju*p

must

be

equal

to or

less

than

Sec.

7.131

Applications

of

the

frictionless

analysis

r87

/4\z

,":\t-j)

^:r'o

I

horizontal

surface

because

this

is a

frictionless

flow analysis

I

?

da

zone

ll

f

.---_._...--

zone

lll

d3

d_

d^

zone

lV

Fig.7.31_(a)Profileofflowwithjump;(b)thexversusrdiagram.

o.t3}do.If

the

depth

is

less

than

this,

the

point

C

moves

downstream

and

the outflow

.8

stays

steady

atLncfis.If

the

depth

dais

greater than 0.138do,

the

flow

remains

steady

but at

a

value

less

than

*^0"

7.13.2

The

efrect

of

the

resistance

of

a dry

bed

With

a

dry bed

the

feathered

edge

travelling

at

a speed

of

?rsis

going

to

be heavily

slowed

down

by

friction.

This

effect

cannot

be

analytically

describgd

but can be

represented

if

working

in a

full characteristic

representation

with

friction.

The finite

difference

forms

alreidy

developed

must

then

be used.

The

moving

edge

must

be

represented

by

treatiog

"

wedfe-shaped

block

of

fluid

as

if

it

were

a solid

and

esiimating

all

ihe

forceJ

acting

upon

ii.

fne

resultant

force

can

then

be

substituted

into

Newton's

second

law to

oUtiin

the

acceleration

of

the

block.

This

can

then

be

used

to

locate

the

block

at the

end

of

the

period and

to find

its

velocity.

The

force

per

unit

width of

the

channel

for

the

wedge

is

whz

wlui

olu

du

z+2-rs:2dt'

h is

the height

of

the

wedge

s is its

length.

i

is

the channel

bed

slope.

188

Unsteady

flow

in

open channels

r can

be calculated

from

the

Chezy

C. It is

given

by

wD2

-c2

and

from this the

value

of duldt can

be

obtained.

Then

lch.7

x: Dtar +

o.s

ff

toO'

x

is the distance moved

by the

wedge

in

time A/ and

t)z: t)r+

ff

nr .

The

value

of s is the

length

of

the

wave

remaining after the integer number of Ax

lengths

have been deducted as shown in Fig. 7.32. s is a

variable

distance

and this

'-

ail'+o;-.>l--A"---*l<;--->l

Fig.7 .32-

The

basis of

calculating the s

value.

technique

is

much

simpler than

tracking the

leading

edge of the

wave

and calculating

the

duldr value

using

a fixed s length.

In the case

of

a downstream

flow the feathered

edge

will

be lubricated and then

there is no need to

make

special

allowances for friction. If desired,

such a case can

be

analysed

using

the frictional

characteristic finite difference method

also.

Finite

difference

methods

8.T

INTRODUCTION

Pressure

transients are easily dealt

with

by use

of

the

method

of

characteristics

and

there is no doubt that

this is the method of choice.

Unsteady open channel hydraulics is only infrequently

solved by the method of

characteristics. Most analysts

prefer

to use

various

techniques

based upon the finite

difference integration

of

the basic differential

equations. The

reasons

for

their

preference

is historical. It has only been

possible

to solve the

differential equations

of

either transient bounded

flows

or free surface flows

since the computer

became

available.

In the days before the small microcomputer,

the

cost

of operating

main

frame computers was

a controlling factor

in the

choice

of method.

If the

run

time

of a

program

was

more than a few minutes, the cost

of the analysis was regarded

by many

as

unacceptable. Finite difference methods were

found to be

very

much faster

than

the method of characteristics and were therefore

preferred.

Now

that micro-

computers are

cheaply and easily available and have large

memories which

operate

at extremely high speeds, this

point

has lost a lot

of its force.

At

the time

of writing

this book, 1.989, it is

not

yet

true that

one can

have

a

supercomputer on one's

desk

but

it is

possible

to have a machine

possessing

l Mbyte

of memorywith a

40Mbyte hard

disc

and a

processor

that

will

perform

at a

rate

of 7 million

instructions

per

second

-

All

this for

sums

of money of the order of a few

thousands of

pounds.

What was an

excessively

expensive

run of a

program

now costs

trivial

sums and

perhaps

should

no

longer

be considered as a factor

in

the choice of method. A run

of.2h is not

thought

exceptional and the author

has run

programs

overnight

on his own machine. The

only cost is that of electricity which is negligible.

8.2

CONSERVATION

FORMS

OF

TIIE

UNSTEADY FLOW EQUATTONS

It is

generally

accepted that the forms of the unsteady

flow equations to be used when

performing

either

explicit

or implicit finite difference

integrations

should be

of

characteristics to

be

used

190

Finite

difference

methods

lch.8

conservative

type.

Some

of them

only

work

in a stable

and

non-diffusive

manner

when

the equations

are

of

this

form.

The

lawi

of

fluid

mechanics

are statements

that

certain

properties

are

conserved.

sions

of

these

principles

often

expressed

in the

form

of a

relationship

between

partial

derivatives.

E.2.1

Conservative

forms

of

the St

Venant

equations

These

are

where

A is

the

cross-sectional

area,

p

is the

flow,

i is

the

bed slope

of the

channel

and

j

isthe

energy

loss

per

unit

weight of

fluid

per

unit

distance

along

the

channel.

Z

is the

dejrthof

the

centroid

below

the

water

surface.

E.3

TIfi

LEAP

FROG

METHOD

This

is

the

method

which

anyone

meeting

finite

difference

methods

for

the

first

time

would

invent.

First

substitute

any

paftial differential

by

its finite

difference

form,

i.e.

Sec.

8.31

The leap

frog

method

similar

expressions

can be

written

for 0Al0x

and0Al0t

In

the above

expressions,

i

describes

ihe

location

in time

andi

describes

the

location

in

space

of the

node

which

the

subscripts

are

defining.

Having

substituted

expressions

such

as

these

into the

basic

equations,

i.e. the

continuity

and

dynamic

equations

it

will be

possible

to solve

for

Q,*r,,

and

A,*r,y.

If

these

calculations

are

carried

out

for each

node

in the

ne$;k,';esults

for

a'new

time

level

wilt have

been

found

for

all

node

points- It is of

approximately

second-order

accuracy

and

mass

conservation

is

very

good- The

sotution

that

results

is

a saw-toothed

line

and this

is due

to the

fact

that

the

integration

is

done

in an

overlapping

fashion.

Node

2

is calculated

using

nodes

1',

2

and 3 at

the earlier

time

level;

node 3

is calculated

using

nodes

2,

3 and

4 at

the earlier

time

level;

and

so

on.

So

essentially

there

are two

independent

grids being

used

as

can be seen

from

Fig. 8.1.

This suggests

that

the number

of

Ar

intervals

in

the reach

should

be even.

191

aA,aQ_n

:TF-

-v

r

dt

ox

,

X.

*.Go'.f)

+

,A(i-

i)

:

o

aQ

:LQ

and

y:L9-

Ox

Lx

0t

Lt

Then

represent

the

finite

difference

form

by some

formula

which

description

of

it.

The

formula

chosen

for the

leap

frog method

is:

LQ

_Q,,,i*,,-

Q,,i-,

Al

ZLx

and

LQ

-Q,*r,-

Q,*r,

Lt

ZLt

Fig. 8.1-

The

two

independent

grids

of

the leap

frog

method'

Itisveryimportantto

observe

the

requirementof

the Courant

[16]

condition,

i.e.

Ltl

-s-

At-lrttl

'

If this condition

is

complied

with,

the scheme

is non-diffusive

and

there

is

no

domain

1 of

finite

difference

domain

2

of

finite

difference

is

a

good

i+? i+?

i+2

;+t

;+1

41 ;J-',

tta

1, 5

2,5

i-1

I

i+1

j+t

i+1

i+z

i+1

1,4

2,4

j-1

I

..

:

t.

I

r

-

---i

l1

lp

I

t,

i+l

I

i,

i+z

I

1,3

2,3

j-1

L-- --I

i-1,

I

i-1

j+l

i-1

i+2

1,2

2,2

W

i-1,i-2

i,i-z

i+1,i-2

i+2,i-2

x

r92

Finite difference

methods

lch.8

likelihood

of

instability.

This method

is in fairly common

use and

is a satisfactory

method

provided that the

points mentioned

above are observed.

Lax

[17]

produced

a first-order scheme

which

was of

first-order

accuracy

and

diffusive.

His student,

Wendroff,

produced a

second-order

scheme

and

the combi-

nation

of the

two

-

the Lax-Wendroff

[9]

scheme

-

is

generally regarded

as being

the

best

finite

difference scheme of

all

the many available.

8.4

THE

LAX-WENDROFF

[9]

METHOD

It

is necessary

to express

the equations

to be integrated in conservative

form

(see the

leap

frog method

above).

Requoting

the conservative

forms

of the St

Venant

equations

from

before:

r*,tw denote

the column

vector

l: l,

O,

equal to the

vector

lr'

*3''"1=

lll

"",

r

be equar

to the

vector

l*t;

to'l:

I

;;l

*" two

st

venant

equations

become

aw

aE(w)

ar

+?*z:o

The following

finite

difference

form

can

be used to

obtain

the

required

values

at the

later time

level

t

* Lt:

wi+t,j-

wi,i

-

fi{oL,r,.j+

TLjLx)

#.Nu,:o'

X.U*(*'.f)+sA(i-i):o

+ 0.s

(t*)'*

i*o.s(LF;,i+o.s

* T,,i*o.rM)

-

R,r-o.t(

LF,,i

-o.,

* T,.i

-o.rAr)

+ Pt,i Mtl

Sec.

8.51

Implicit

schemes

where

the difference

operators

have

the

following

definitions:

Ai,i*o.s:(

),,i*t-(

)t.i

Aj,;

:

0.5[(

),.r*,

-

(

),,;-

r]

,

(

),,i*o.r:0'5[(

),,i*,

*

(

),,r]

'

.r

-lsei

-

sAil

'i'j

I

o

I'

T,,

1

*

o.s

:

0.5(T,.,

* T

t,i

* r)

and

D

--ar

t

i,i_

il

The

matrix

R

is

defined

as

lzete

ga(AilaA

-

Q'tA'l

R:11

o

I

The

Lax-Wendroff

scheme

is

capable

of

passing a shock

-

in

the

case

of

open

channel

flow

this

corresponds

to

passing

a travetling

surge

or

jump

-

and

is

capable

of

using

vOry

much

laiger

Ar

vilues

than

is

possible

when using

the

method

of

charactiristiis.

The

schJme

is of

second-order

accuracy

and

is not

dissipative;

so

it

can

pass

a shock-

a

travelling

surge

or

jump

in the

context

of

an

open

channel-

but

it will

not smooth

out

a

pertuibati-on

generated

by

initial

conditions

which

do

not

fit

the

basic

flow

equationi

weil. The

boundary

conditions,

in

almost

all

explicit

finite

difference

schemes,

must

be

computed

using

the

techniques

of

the

method

of

characteristics.

8.5

MPLICIT

SCIIEMES

The

manner

in which

a

method

can

become

implicit

instead

of

explicit

can

be seen

from

Fig.

8.2.

It can

be seen

that,

by interpolating

between

the

values of

variables

at

the current

time

level and

the

futuri

time

ievel,

values of

variables

at the

base

of

characteristics

can be obtained.

Using

these

values

in

conjunction

with

the

characteristic

equations

gives

expressions

for

ihe

variables

at the

level.

The

interpolations

for

values .i the

base

of

the characteristics

will be

in

terms of

the

values

at

the

current

time

level

and

at

the

future

time

level.

This results

in equations

containing

the

variables at

the

required

future

time

for

positions

j-

t,i

and

i+

l; so

there

will six

variables

ur,j-r,ut,i,ur,i+t,dt,i-t,d,randd,,i*ttosolvebutonlytwoequationsare

193

194

Finite

difference

methods

lch.8

Fig.

8.2

-

Timewise interpolation.

available

to

solve them. However, if the

set

of equations

for the next

7

value

are

obtained, four

of the

variables

that occurred

in the first

set

will

be repeated,

i.e. u,,,,

Di.i+r, d,,randdi,i*iso

it can

be seen that

only two newvariables

are introduced into

each

set

of

equations for each node.

As there

are n

-

2 internal nodes

in the reach

length

(where

n is

the number of nodes

in the

reach counting

the first node as 1)

and

there

are two

boundary nodes

which

each

give

two equations from

the two

characteristics

and

the boundary conditions

themselves, it follows that

there will

be

sufficient

equations

to soive the entire

future

time ievei. As

each

variabie

appears in

at

least

six

equations

(except

for the

variables

at the

boundary nodes which

appear

only

in

four)

the

variables

are in implicit

form and the

solution is to

solve all the

equations

simultaneously

for theZnvariables

in the

problem.

The matrix

to be

solved

is

a banded

matrix

and is relatively

easy to

solve and it is

possible

to use

a standard

procedure.

The

solution

is

rapid and the

process

does not require

much storage.

The

foregoing

uses i

to denote spatial location

and

7

to denote

timewise location

of the

node.

8.6 THE

STRELKOV METHOD

The

former method is based

on characteristics

but there

is no need

to employ them.

Strelkov suggested two methods:

one suitable

for waterhammer

analysis

and the

other

for open channel

problems.

The author

is not

persuaded

that finite

difference

methods

are well

suited

to waterhammer problems

and believes that

characteristics

methods

are

superior in this context.

The

reason for this is

that

events in

pipe

networks,

such as

valve

closures

or

pump

trips, occur in

small times which

are small

multiples

of the

pipeperiod2Llc

and, if large

Atvalueswere

to be usedin

an analysis,

which

is

the main

reason for using explicit

or implicit

finite difference

methods, very

poor

modelling

of such

valve

or

pump

behaviour

would

result.

For

completeness,

both of the

Strelkov methods

are

given

below

(Fig.

8.3).

E.6.1

The first

Strelkov

method

In

the first method,

Strelkov

suggested

the following

finite difference

forms for the

partial

diffential coefficients that

occur in the waterhammer

equations:

Ah

:(hr*r,,*r+

ht*r,,)

-

(h,,,*r+

hr,,)

dx

2Lx

Sec.

8.61

The Strelkov

method

195

(b)

(a)

Fig.

8.3

-

The

Strelkov

methods:

(a)

the first

method;

(b)

the

second

method.

ah_

0t

ZLt

AD_(u,*r,i*t+u,*t

0x

2Lx

Ou

_

(u,*r,r*, +

u,*r,r)

-

(u,,r*, * ur,r)

Ot

ZLt

Also

,

-

4*r,L*

ui,i

These expressions

must

into

the

waterhammer

equations:

Ah

czdu

ar

+'ar*?ar:u

'

Ah

uOu

10u

zfulul

*+;a'.;i*i:s

The

result

will be

two equations

of

the

type

ut,i+t*

u;*r,i*r

-

ahi,j*r*

Fht*r,;*r

*

T

:0,

Dt,1+t*

ur*

t,j+r

-

6hr,i*r*

ehr*r,i*

r

*

1': 0

where o,

F,

6, e,

y

and

l,

are constants

involving

ui,i,,

ui*r,i,,

ht,l

and

h,*r,,

These

last

two

equations

contain

four

unknowns.

As

with

the

when each

node

is eiamined

in turn,

it will

be seen

that

only

two

new

variables

are

introduced

for

each

node;

so there

are

only

2(n-

2) equations

available

from

the

internal

nodes.

As

with

the

four

additional

equations

can

be

Fox J.A. Transient Flow In Pipes, Open Channels And Sewers

Подождите немного. Документ загружается.