Solution of the St Venant Equations / Shallow

Download Report

Transcript Solution of the St Venant Equations / Shallow

Solution of the St Venant
Equations / Shallow-Water
equations of open channel flow
Dr Andrew Sleigh
School of Civil Engineering
University of Leeds, UK
www.efm.leeds.ac.uk/CIVE/UChile
St. Venant Equations
d
 d u
g
u

 g S o  S f 
x
 x t
A
u
d
u
A
b
0
x
x
t
St Venant Assumptions of 1-D Flow

Flow is one-dimensional i.e. the velocity is uniform over the cross
section and the water level across the section is horizontal.

The streamline curvature is small and vertical accelerations are
negligible, hence pressure is hydrostatic.

The effects of boundary friction and turbulence can be accounted for
through simple resistance laws analogous to those for steady flow.

The average channel bed slope is small so that the cosine of the
angle it makes with the horizontal is approximately 1.

Cunge J A : Practical Aspect of Computational River Hydraulics
Control Volume
y
c
re
di
n
tio
A1
V1
f
of
w
lo
h1
free surface
h
2
A2
V2

z1
x1
s
z2
x
x2
bed
Continuity

CV in x,t plane



between cross sections x=x1 and x=x2
between times t=t1 and t=t2
conservation of mass
 A   A  dx   Q
x2
x1
t2
t2
t1
t1
x2
 Qx1  dt  0
Momentum

Conservation of momentum
2
2










uA

uA
dx

u
A

u
Ax 2  dt
t2
t1
x1
x1
t1
x2
t2
g
t2
g 
t2
g
t2
t1
t1
t1
I 
1 x1

x2
x1
 I1 x 2  dt
I 2 dx dt
 AS
x2
x1
o
 S f  dx dt
Geometric Change Terms

I1  
0

()
vertical change
in cross-section
h x 
b
h
h-
d

hx    x,  d
change in width along the length of the
channel.
h x 
  
I 2   hx    
 d
0
 x  hho
Integral / Differential Forms



Integral form do not require that any flow
variable is continuous
We will see later finite difference methods
based on this integral form.
Can derive differential form … but


Must assume variables are continuous and
diferentable
Replace variable with Taylor’s series
 At 2
A
 2 A t 2
  At1 
t  2
 ...
t
t 2
Differential form

-
A Q

0
t x

Q   Q 2
 
 gI1   gASo  S f   gI2
t x  A

gI 1 
h
 gAx   gI 2
x
x
Q 
 h

 uQ  gA  So   gASf  0
t x
 x

In terms of Q(x,t) and h(x,t):

Where b = b(h), A=A(h)
Ah  A h
h

b
t
h t
t
h 1 Q

0
t b x
Q   Q 2 
h
    gA  gAS f  So   0
t x  A 
x
remember b = b(h), A=A(h)
Each of these forms are a set of non-linear differential equations which do not
have any analytical solution.
The only way to solve them is by numerical integration.
In term of u(x,t) and h(x,t):

using
 A h  A 

Q
u
A
u
 A u
 A  u
 

x
x
x
x

h

x
 x  hcons tant 

h A u
h u  A 

u   
0
t b x
x b  x  hcons tant
u
u
h
u
g
 g S f  S o   0
t
x
x
Characteristic Form




The St Venant equations may be written in a
quite different form know as the Charateristic
Form.
Writing the equations in this form enables
some properties and behaviour of the St
Venant equations become clearer.
It will also help identify some stability criteria
for numerical integration
will help with the definition of boundary
conditions.
Characteristic form

Consider a prismatic channel
h A u
h

u
0
t b x
x
u
u
h
u
g
 g S f  S o   0
t
x
x
Wave speed

Consider the speed, c, of a wave travelling in
the fluid.
c g

A
b
with respect to x and t gives:
c
h
2c  g
x
x
2
c
c
u
 2u  c
0
t
x
x
2c
c
h
g
t
t
u
c
u
 2c  u
 g S f  S o   0
t
x
x
Combining

Adding equations (3) and (4) gives
v
v
c
c
 v  c   2  2v  c   g So  S f 
t
x
t
x

Subtracting equations (3) and (4) gives
v
v
c
c
 v  c   2  2v  c   g S o  S f 
t
x
t
x
Characteristics

Equations (5) and (6) can be rearranged to
give respectively
v  c  v  2c   v  2c   g So  S f 
x
t
v  c  v  2c   v  2c   g So  S f 
x
t
Total differential


For some function f of x and t
df sense.
f dx f
In a physical





dt
fx is
dt some
t property of the flow e.g.
If the variable
surface level,
if an observer is moving at velocity v,
the observer will see the surface level change
only with time relative to the observers' position.
The characteristic form of the St Venant
equations

If we take f = (v + 2c)

Total differential is
d v  2c   v  2c  dx  v  2c 


dt
x
dt
t

Compare with
v  c  v  2c   v  2c   g So  S f 
x

Clearly

and
dx
 v  c 
dt
d v  2c 
 g S o  S f 
dt
t
Characteristic form of the St Venant
Equations

These pairs are known as the Characteristic
form of the St Venant Equations
d v  2c 
 g S o  S f  for
dt
dx
 v  c 
dt
d v  2c 
 g S o  S f 
dt
dx
 v  c 
dt
d v  2c 
 g S o  S f 
dt
for
for
dx
 v  c 
dt
Meaning of the characteristics
The paths of these observers that we have talked about
can be represented by lines on this graph.
c2
dt
1

dx v  c 
t
dt
1

dx v  c 
8
c1
6
c4
7
c0
5
c5
4
3
2
Zone of quiet
1
x
Information paths






it takes time for information to travel
E.g. a flood wave at u/s end
The channel downstream will not receive the
flood for some time.
For how long?
The line C0 represents the velocity of flood
wave
Everything below C0 is zone of quite
Zones of Dependence and Influence

The idea that characteristics carry information
at a certain speed gives two important
concepts
c2
c1
t
c2
P
Domain of dependence of P
c1
Zone of influence of Q
Q
x
Stability

These zones imply a significance to
numerical methods and stability of any
solver

The numerical method must take only
information from within the domain of
dependence of P

this limits the size of time step
Calculation with characteristics


If we know the solution at points 3 and 5
Can determine the solution at point 6.

For C5-6 then
d v  2c 
 g i  j 
dt
1

dx v5  c5 
v6  2c6   v5  2c5   t g i  j 
dt

For C3-6 then
d v  2c 
 g i  j 
dt
dt
1

dx v3  c3 
v
6
 2c6   v3  2c3   t g i  j 
Characteristics solution

Adding equations
v6 

v5  v3
 c5  c3   t g i  j 
2
Subtracting
c6 
v5  v3 c5  c3

4
2
MOC on Rectangular Grid
t
P
t+t
t
t
W
L
O
R
E
x
x
P
t+t
t
t
W
L
O
x L
x
R
x R
x
E
Midstream discretisation

Away from boundaries
dx
 v L  c L 
dt
x L
 v O  c O 
t
P
t+t
dx
 v P  c P 
dt
t
t
W
L
O
x L
x
x L
vO  vW 
v L  vO 
x
R
x R
E
x R
 v O  c O 
t
x
x R
vO  v E 
v R  vO 
x
Solution
d v  2c 
 g So  Sf
dt

vP  2cP   vL  2cL   t g So  Sf L
d v  2c 
 g So  Sf
dt

vP  2cP   vR  2cR   t g So  Sf R
vP 
vL  vR
t g
So  Sf L  So  Sf R 
 cL  cR  
2
2
cP 
vL  vR cL  cR t g
So  Sf L  So  Sf


4
2
4
R 
Stability

Considering characteristics
C2
t+t
t
t
P
dt
1
dx = v+c
L
C1
dx
 v  c 
dt
dt
1
dx = v-c
M
R
x
t 
v  c 
x
x
t  0.9
v  c max
Boundary conditions
Upstream boundary:

Downstream boundary:

forward characteristic
t+t
t
P
P
t+t
boundary

backward characteristic
boundary

t
t
t
O
R
E
W
L
O
x R
x
x L
x
Boundary Conditions

The second equation is a



boundary condition equation
Upstream depth boundary
c p  ghp
hP = H(t)
vP  2cP   vR  2cR   t g So  Sf R
vP  vR  2cR  cP   t g So  Sf R
Spuer-critical - mid

Right characteristic moves to left
P
t+t
t
t
W
L
R
O
E
x R
x L

x
x
solution method is exactly the same
Super-critical - upstream
No characteristics
P
boundary

L
R
t
O
E
Super-critical downstream
2 characteristics

No boundary condition
P
t+t
boundary

t
t
W
L
R
O
x R
x L
x
Finite Difference Schemes

Two basics classes




Implicit
Explicit
Commercial packages use implicit
Explict



for high accuracy (sometimes!)
Testing / understanding behaviour
Class examples!
Which Equations ?



Not always clear what equations a being used!
What are the shallow water equation?
We will look at schemes based on the integral
equations:
x1  At 2   At1  dx  t1 Qx2  Qx1  dt  0
x2
t2
2
2
x1 uAt 2  uAt1  dx  t1 u Ax1  u Ax 2  dt
x2
t2
g
t2
g 
t2
g
t2
t1
t1
t1
I 
1 x1

x2
x1
 I1 x 2  dt
I 2 dx dt
 AS
x2
x1
o
 S f  dx dt
Homogeneous Integral Equations

Without the gravity / frictions terms
 A
x2
x1
t2
 At1  dx   Qx 2  Qx1  dt  0
t2
t1
2
2










Q

Q
dx

u
A

gI

u
A  I1 x1  dt  0
1 x2
t1
x1 t 2
t1
x2
t2
Grid based

Consider the grid …
B
n+1
t
t
n
C

f

A
D
x
i-1
x
i
i+1
Integrate around the cell

Considering the cell ABCD,

Integral can be written in the general vector form :
  fdx  G  f dt   0
ABCD
Q
f  
 A
 Q 

G f    Q 2
 gI1 

 A

Gridpoints / Variables

Variables are all known or will be calculated
at the grid points



xi represents the value of x at position i
tn represents the value of t at position n
Derivation approximates values:
Gxi , t   Gin1  1   Gin
0  1
f x, tn   fi n1  1   fi n
0   1
Substitute in approximations

And the equation become
 f
xi 1
xi
n 1
i 1
 

t n1
tn
 f
n 1
i 1

 1   f i n 1  f i n1  1   f i n dx
G
n 1
i 1
 1   f i

 

 1   Gin1  Gin 1  1   Gin dt  0
n 1
 f
n
i 1
 1   f i
 
n
x

 Gin11  1   Gin1  Gin1  1   Gin t  0
Difference equation

f
n 1
i 1

Divide through byΔx Δt
 
 
 
And can see that it is an approximation of
f G  f 

0
t
x


 1   f i n1  f i n1  1   f i n Gin11  1   Gin1  Gin1  1   Gin

0
t
x
i.e. starting with integral form
discretisation is also valid for diferential form
Several schemes

This is a general discretisation scheme



Vary the parameters ψ and θ
Get a family of different finite difference schemes
Features are



They are implicit for values of ψ > 0. else explicit.
They link together only adjacent nodes.
Space interval can vary – no loss of accuracy.
They are first order, except for the special case of
ψ=θ=0.5 when they are second order.
Preissmann Scheme



ψ = 0.5 gives Preissmann 4-point scheme
Time derivative
f in11  f i n1  f in1  f i n
f

t
2t

 

Space derivative

 
G f  Gin11  1   Gin1  Gin1  1   Gin

x
x

Equations become …


x n1
x n
Ai  Ain11 
Ai  Ain1  t  Qin11  Qin1  1    Qin1  Qin  0
2
2






x n1
x n
Qi1  Qin1 
Qi1  Qin
2
2
n1
n1
n
n
2
2
   Q2

  Q2







Q
Q
 t   
 gI1   
 gI1    1    
 gI1   
 gI1  

 A

A
A
   A






i

1
i
i

1
i










  I 2  AS f  ASo n1   I 2  AS f  ASo n1 
j
j 1
gtx 

n
n 1
2 

1     I 2  AS f  ASo  j   I 2  AS f  ASo  j 1




  0


More common form

The terms I1 and I2 are often difficult
(expensive / time consuming)


(when originally attempted)
Usual form used in packages
h 1 Q

0
t b x
Q   Q 2 
h
    gA  gAS f  So   0
t x  A 
x
Priessmann

Function represented by,
f 


0.5 < θ ≤ 1
Continuity:


f
2
n 1
i 1
 fi
n 1

1  n

f
n

f
i 1
i
2

h 1 Q

0
t b x






1  hin11  hin1 hin1  hin  2  Qin11  Qin1  1    Qin1  Qin

 

0
n1
n1
n
n
2
t
t  x  bi  bi1  1    bi  bi1


Priessmann

Momentum Equation
Q   Q 2 
h
    gA  gAS f  So   0
t x  A 
x
1  Qin11  Qin1 Qin1  Qin 



2
t
t

    Q 2  n1  Q 2  n1  1      Q 2  n  Q 2  n 


  
  



   

x   A i 1  A i 
 x   A i 1  A i 




1  n
  n1

n 1
 g  Ai 1  Ai 
Ai1  Ain E  0
2
2





Source term

So-Sf
 hin11  hin1
hin1  hin 

E  
 1   
x
x 

n
n
n n
  Qin11 Qin11  Qin1 Qin1
Q
Q

Q
i 1 i 1
i Qi


 
 1   

2
2

   
     
n 2
n
  K n1 2  K n1 2
K

K
i
i
  i1
 1    i1

2
2

2


Unknowns

There are the 4 unknowns
hnj1, hnj11, Qnj1, Qnj11,

Plus k, b, A which are functions of h and Q
n1
j
n1
j 1
n1
j
n1
j 1
h ,h ,Q ,Q ,
Need to linearise

“Linearise” equations
Ahnj1  Bhnj11  CQnj1  DQnj11  RHS

A, B, C, D and RHS are funtion of the
unknowns.

But not strongly
Boundary Conditions

In implicit scheme specify h or Q

Use characteristics to decide appropriately

Or relation between h and Q

2N unknowns (h, Q at each node)


2N – 2 equations from internal points
2 boundary equations
Junction

At a junction each chanel share the same
node


hjuntion 1 = hjuntion 2 = hjuntion 3 …..= hjuntion n = h
Continuity

Sum of inflow and outflow equal to zero
Iterative solution

Need to iterate updating coeficients

Cunge says …


iterates rapidly one or two iterations
Newton-Raphson methods used in packages
Stability

Formally unconditionally stable for all time
steps




0.5 < θ ≤ 1
Further away from Cr = 1 less accurate
Cr = 20 is common
Because of linearisation may fail for extreme
flows or those that are too far from original
assumptions
Explicit Schemes

Not used in simulations of real rivers
Time step limitations.

Important features of explicit schemes




they are simple to implement
allow experiment with weights, time-step and
space-step to understand behaviour of the
solution.
There are MANY schemes
Leap-Frog

Earliest scheme aplied to wave equations

Spatial derivative
f f
f

t
2t
n
i 1

gives
n
i 1
Temporal derivative
G f  Gin1  Gin1

x
2x
n1


A
t
n1
n
i
f i   n1   f i 
Gin1  Gin1
2x
Qi 


Stability


All explict scheme have time step limit
Courant confition (CFL)
dx
 v  c 
dt

Cr < 1
x
t 
v  c 
Lax Explict Scheme

Similar to leap from with weighting, α, in time

n
n


f

f
n 1
n
i 1
i 1

f i   f i 1  1   
2
f



t
t
For 0 ≤ α ≤ 1

Spatial derivative
G f  Gin1  Gin1

x
2x
Lax Explict Scheme

Leads to function solution
n1
n
n


A
f

f
t
n1
n
n
n
i
i 1
i 1
f i   n1   f i  1   

Gi1  Gi1
2
2x
Qi 


Boundary conditions must be applied using
method of characteristics

Some useful texts





Rather old- still basis of many commercial programs.
Cunge, J.A., Holley, F.M. and Verwey, A. (1980):
Practical Aspects of Computational River Hydraulics
Mahmood and Yevjevich (1975): Unsteady Flow in
Open Channels - Fort Collins, Colorado
Liggett & Cunge (1975)
Preissmann, A. (1960): Propogation des
Intumescenes dans les Canaue et Rivieres - 1st
Congress de l'Assoc. Francaise de Calcul, Grenoble