1

The St Venant Equations

Dr P A Sleigh
Dr I M Goodwill

School of Civil Engineering,
University of Leeds
March 2000

1 THE DERIVATION OF THE CONTINUITY EQUATION ...........................................................................................1
2 THE DERIVATION OF THE DYNAMIC OR MOMENTUM EQUATION. ...............................................................2
3 THE SOLUTION OF THE ST VENANT EQUATIONS .................................................................................................4
3.1 AN MORE CONVENIENT FORM OF THE EQUATIONS...........................................................................................................4
3.2 A DIVERSION . .............................................................................................................................................................5
3.3 THE CHARACTERISTIC FORM OF THE ST VENANT EQUATIONS .........................................................................................5
3.4 THE MEANING OF THE CHARACTERISTIC FORM................................................................................................................7
3.5 AN EXAMPLE COMPUTATION USING THE CHARACTERISTIC EQUATIONS...........................................................................8
3.6 THE DOMAIN OF DEPENDENCE AND ZONE OF INFLUENCE.................................................................................................9
3.7 A FIRST NUMERICAL METHOD OF SOLUTION....................................................................................................................9
3.7.1 The Staggered Grid Method .................................................................................................................................10
3.8 STABILITY AND ACCURACY ..........................................................................................................................................11
3.9 A SECOND NUMERICAL SOLUTION METHOD...................................................................................................................13
3.9.1 The method of characteristics on a rectangular grid ...........................................................................................13
3.9.2 The stability criteria .............................................................................................................................................15
3.10 BOUNDARY CONDITIONS ...............................................................................................................................................16
3.10.1 Depth specified - Upstream boundary..................................................................................................................16
3.10.2 Depth specified - Downstream boundary .............................................................................................................16
3.10.3 Discharge specified - Upstream boundary ...........................................................................................................17
3.10.4 Discharge specified - Downstream boundary.......................................................................................................17
3.11 SUPER-CRITICAL FLOW.................................................................................................................................................17
3.11.1 Boundary conditions in super-critical flow ..........................................................................................................18

1 The derivation of the continuity equation

Consider a short length, x, of channel

Datum
z
Q
1
Q
2
x
bed
Surface
h

Figure 1 - Short length of channel
2

The following symbols are used in this derivation:
A = the cross-sectional area of the section
h = depth of flow at the section
z = elevation of surface above a datum at the section
v = mean velocity at the section
Q = discharge at the section
b = width of the top of the section
x = position of the section measured from the upstream end
t = time
g = acceleration due to gravity
= mass density of the fluid

Others symbols are defined in the text at the point when they are introduced.

Assuming that there is no lateral inflow, then
x
x
Q
Q Q

1 2


This has the partial derivative since Q is changing with both x and time, t.

Now the volume of water between the sections 1 and 2 is increasing as a rate of
x
t
h
b

where b is the top width,

As cross-sectional area A = bh then this is equivalent to
x
t
A

The terms are equal in magnitude but of opposite sign, so
0

x
t
h
b x
x
Q


as
( )
x
Av
x
Q

then


0

t
h
b
x
v
A
x
A
v
(1)
This is the continuity equation

2 The derivation of the dynamic or momentum equation.
By applying Newtons 2
nd
law to our elemental length of channel we have

1
]
1

t
v
x
v
v x A
dt
dv
x A

on accelerati mass Force

3
since v varies with both space (x) and time (t)

Consider the external forces which cause this acceleration. These are, in the simplest case, three


x
H

change in static pressure

F frictional resistance of channel walls and bed
g gravity force (the weight)

If is the bed slope (measured positive as the bed rises from downstream to up see Figure 1) then the
sum of these three forces is

sin cos x gA x F x
x
H
+

For small bed slopes, , then cos = 1 and sin = = i so

xi gA x F x
x
H
+

Now

x
h
gA
x
H

and
gAj F

where j is energy loss / unit length of channel / unit weight of fluid.

Equation these external forces to the change in momentum yields

i x gA j gA x
x
h
gA
t
v
x
v
v x A +

,
_

Rearranging gives


( ) j i g
t
v
x
v
v
x
h
g

(2)
This is the dynamic, or momentum, equation.

Using the Chezy expression, j, can be written

m C
v
j
2
2

where C is the ChezyC and m is the hydraulic mean radius given by

p
A
m
perimeter Wetted
Area


4
3 The solution of the St Venant equations

The St Venant equations are:

The continuity equation

0

t
h
b
x
v
A
x
A
v
(1)
The dynamic, or momentum, equation

( ) j i g
t
v
x
v
v
x
h
g

(2)

The St Venant equations cannot be solved explicitly except by making some very large assumptions
which are unrealistic for most situations. Therefore numerical techniques have to be used.

Equations 1 and 2 are not in a readily usable form to solve, so the first task is to rearrange them.

The problem is that A is a function of h so we need to re-write the equations using h and to arrange them
so we can solve for the two unknowns h and u.

Assume the simplest channel that is, rectangular cross-section of constant slope.
3.1 An more convenient form of the equations

It is easier if we write the equations in terms of celerity, c, rather than h using this equation for celerity of
a small gravity wave in a rectangular channel of depth h:


g
c
h
gh c
2

All the terms involving h of equations (1) and (2) can be replaced by terms containing c making use of


( ) ( )
dc
g
c
dh
dc c c d gh d
2
2
2

For the continuity equation (1), where A = bh, we get


0 2 2
0
2 2
0
2

t
c
x
v
c
x
c
v
t
c
b
g
c
x
v
g
c
b
x
c
bv
g
c
t
h
b
x
v
bh
x
h
bv

(3)
5

For the dynamic equation (2) we get

( ) j i g
t
v
x
v
v
x
c
c

2
(4)

Adding equations (3) and (4) gives

( ) ( ) ( ) j i g
x
c
c v
t
c
x
v
c v
t
v

+ +

+ +

2 2
(5)
Subtracting equation (4) from (3) gives

( ) ( ) ( ) j i g
x
c
c v
t
c
x
v
c v
t
v

2 2
(6)

Equations (5) and (6) can be rearranged to give respectively

( )
( ) ( )
( ) j i g
t
c v
x
c v
c v

+
+

+
+
2 2

(7)

( )
( ) ( )
( ) j i g
t
c v
x
c v
c v

2 2

(8)

3.2 A diversion .

In differential calculus, one of the basic equations of partial differentiation is

t dt
dx
x dt
d

where x and t are independent variables.

In open channel hydraulics we may think of this in a physical sense. If the variable is some property of
the flow e.g. surface level, that varies with both distance (x) and time (t) then if an observer is moving at
velocity
dt
dx
v , the observer will see the surface level change only with time relative to the observers
position.

3.3 The characteristic form of the St Venant equations

As could be any variable, we can write it as = (v + 2c) then


( ) ( ) ( )
t
c v
dt
dx
x
c v
dt
c v d

+
+

+ 2 2 2


Compare this equation with the left-hand side of equation (7)
6

They are equivalent if ( ) c v
dt
dx
+ , so equation (7) can be re-written as the total differential


( )
( ) j i g
dt
c v d

+ 2

(9)
And similarly equation (8) can be re-written in total differential form if ( ) c v
dt
dx
as


( )
( ) j i g
dt
c v d

2

(10)

Physically these can be interpreted as, for equation (9), if an observer is moving at velocity (v+c) then
they see / experience changes given by (9).

And for equation (10), if the observer is moving at velocity (v-c) then they will see / experience changes
given by (10).


These pairs are known as the Characteristic form of the St Venant Equations


( )
( ) j i g
dt
c v d

+ 2
for ( ) c v
dt
dx
+
(11)
and

( )
( ) j i g
dt
c v d

2
for ( ) c v
dt
dx

(12)


You may also see them quoted as one pair using the t symbol i.e.


( )
( ) j i g
dt
c v d

t 2
for ( ) c v
dt
dx
t
The solution of this system of equations is quite straight forward - particularly using a computer.
7
3.4 The meaning of the characteristic form

Consider the graph below which is in the x - t plane.
1
2
3
Zone of quiet
x
t
c5
c4
c2
c1
c0
4
5
6
8
7

Figure 2, Characteristics on an x-t plane.

The paths of these observers that we have talked about can be represented by lines on this graph. That is
to say, lines of gradient
( ) c v dx
dt
+

1
or
( ) c v dx
dt

1
on this graph are the paths of an observer. They may
start at any point and are known as characteristic lines or simply characteristics. Along these lines the
associated characteristic equation applies. It is possible to solve the characteristic equations at the point
where they cross at the same time this is the solution to the St Venant equations. This technique is
known as the Method of Characteristics.

In practice the method of characteristics is becomes cumbersome to apply exactly on a computer,
however there are many different ways that have proven successful, two of which will be shown below.
Before that one very important results can be obtained from looking at the significance of the
characteristics and their slopes that of numerical stability.

Consider again the graph above and a disturbance at the upstream end of a channel (e.g. a flood wave).
The part of the channel some way along from the end will not receive the disturbance for some time. The
time it takes depends on the velocity that the information from the disturbance travels. In the graph of
figure 1 the line C0 represents this velocity. Everything below that line represent the channel waiting for
the disturbance this is known as the zone of quite. C0 is a characteristic line. It is a forward
characteristic as are line C1 and C2, they have gradient of the form v+c, (note though that each has a
different gradient.) Lines C3, C4 and C5 are negative (or backward) characteristics and have gradients
of the form v-c.
The values v and c at point 5 is influenced by events and conditions at points 0, 1, 2 and 4. But not by
values at any other points outside this region the values of v and c do not affect the values of v and c at
point 5.

8
Similarly at point 2, say, the values of v and c will influence the values of v and c at points 5, 6 and 7 etc.
but not at point 4.

3.5 An example computation using the characteristic equations

Assume we know the solution at points 3 and 5 and we want to determine the solution at point 6. We
know that a forward characteristic through 5 intersects with a backward characteristic through 3. If the
conditions at point 3 and 5 are respectively v
3
, c
3
and v
5
, c
5
then from (11) and (12) the forward
characteristic from 5 starts off with gradient
( )
5 5
1
c v +
and the negative characteristic from 3 starts off
with gradient
( )
3 3
1
c v
.

Along the forward characteristic applies

( )
( ) j i g
dt
c v d

+ 2

or
( ) ( ) ( ) j i g t c v c v + +
5 5 6 6
2 2
(13)
and along the backward characteristic

( )
( ) j i g
dt
c v d

2

or
( ) ( ) ( ) j i g t c v c v
3 3 6 6
2 2
(14)

Equations (13) and (14) can be solved simultaneously to give the conditions at point 6.

Adding equations (13) and (14) gives
( ) ( ) j i g t c c
v v
v + +
+

3 5
3 5
6
2

subtracting gives


2 4
3 5 3 5
6
c c v v
c
+
+

i.e. the solution at point 6 has been found.


9
3.6 The domain of dependence and zone of influence

The above concepts of which points influence other lead to two important ideas that is the domain of
dependence and zone of influence. These are shown diagrammatically below.

c1
c2
c1
c2
Q x
t
Domain of dependence of P
P
Zone of influence of Q


Figure 3, Domain of dependence and zone of influence.

Conditions at point P are determined solely by the conditions bounded within the two characteristics C1
and C2. Likewise the conditions at point Q determine the conditions within the area bounded by the
characteristic lines

The figure left shows the same regions when the characteristic lines
are shown extended through the point.





The significance to numerical methods and stability, is that the
numerical method must take only information from within the
domain of dependence of P when determining the value of P. In
practise this limits the size of time step possible hence solutions
usually take longer to compute than desired.




Figure 4

3.7 A first numerical method of solution

Two commonly used numerical methods will be demonstrated in this note. One is based on a staggered
grid of solution point and doesnt actually use the characteristic equations and the second on a rectangular
grid which does make explicit use of these. Both methods are referred to as the Method of Characteristics
although they are both really just approximations. They do use the ideas of where information if
travelling from - which is itself determined by the characteristics.

c1 c2
x
t
Domain of dependence of P
P
Zone of influence of P
10
3.7.1 The Staggered Grid Method

Consider the grid shown in Figure 5, in which the nodes are spaced at regular x intervals along the
direction of the x axis and at t intervals along the t axis, the nodes are staggered a distance x/2 between
t and t+t . (Note that the t spacing may not be equal - more about that below.)

We wish to calculate the solutions at point P which is at time t+t, we know all the solutions values at the
previous time, at time t.

The two points to the left and right of P, but at the previous time level, t, are points L and M respectively.
Drawing two characteristics through these gives the arrangement shown in Figure 5. The partial
derivatives of the dependent variables can be written


x
v v
x
v
L R

x
c c
x
c
L R

and

t
v v
t
v
M P

t
c c
t
c
M P

where point M is mid way between L and R and mean values can be used.

x
t
t
t
t+ t
P
L R
M
C1
C2

Figure 5, A staggered grid

Substituting these into our rearranged St Venant equation (3) gives


( ) ( ) ( )
( ) j i g
t
v v
x
v v
v
x
c c
c
M P L R
M
L R
M

2
(15)
This equation is explicit in v
P
i.e. we can calculate v
P
directly

( ) ( ) ( ) [ ]
M R L M R L M M P
j i x g v v v c c c
x
t
v v + +

+ 2

And similarly for equation (4)

( ) ( ) ( )
0 2 2

t
c c
x
c c
v
x
v v
c
M P L R
M
L R
M

(16)
which is explicit in c
p
giving
11

( ) ( ) [ ]
R L M R L M M P
v v c c c v
x
t
c c +

+ 2
2
1


These expressions for v
p
and c
p
enable the solution to be determined since all variable on the right hand
side are from the known conditions at the starting time level t. they can be used to calculate the solutions
at every advanced time point on the t+t level by moving from point to point using each ones lower left
and right point values.

This statement needs qualifying since the expressions are not true for the upstream and downstream
points (when x=0 and x = L). At these boundaries we must derive special equations - known as boundary
conditions to describe exactly what changes are occurring there. The question of boundary conditions
will be addressed later in this note.

3.8 Stability and Accuracy

The method described above is clearly quite simple to implement - although there are a lot of repetitive
calculations involved if we wish to calculate for every point on a channel. A computer makes it
particularly straightforward to calculate all the point and advance (or march) through time updating all
values at every point and at every time level. However one must be extremely careful that the solution
remains stable and accurate.

Stability and accuracy are controlled by the values t and x. We can determine a stability criterion by
considering how information is being passed in our solution - or as we have seen by looking at the
domain of dependence to see where the information at any point comes from.

Consider Figure 6 which shows the characteristics around the point P from our earlier calculation (see
Figure 5)

t
x
t+ t
P
L
R
M
C1 C2
t
1
=
v+c
dt
dx
1
=
v-c
dt
dx


Figure 6, Characteristics around P

The forward characteristic through L has slope
c v dx
dt
+

1
and the negative characteristic through R has
slope
c v dx
dt

1
. From our earlier arguments point P can only receive information from points L and R if
it lies within the domain of dependence of L and R - i.e. in the region bounded by the two characteristics.

12
Therefore to remain within this domain of dependence this criteria for the forward characteristic must be
satisfied

( ) c v
dt
dx
+ >
(17)

This stability condition is often referred to as the Courant or CFL (Courant-Friedrichs-Lewy) condition

and since dx = x/2 for our staggered grid the stability criteria is

( ) c v
t
x
+ >

2

Considering the backward characteristic we get the criteria for stability
( ) c v
dt
dx
>

but this is unnecessary as it is already taken account of by criteria (17) (as (v+c) > (v-c)).

Clearly to apply (17) values of v and p must be known. But which do we choose? The values at L, R and
P are all different (and those at P are not known beforehand). In practice the values at the known level
must be used and some factor ( 0.9) introduced to take into account that the values at P may be so
different to those at L or R as to change the gradient of the characteristics significantly.

To actually use the stability condition in a calculation it is usual to choose a fixed grid spacing, x, then
determine the maximum value of (v+c) for each point on the know time level then calculating a t. So for
the staggered grid the time step would be


( )
max
2
9 . 0
c v
x
t
+

(18)
This ensures stability of calculation for every point.

For accuracy the stability criteria must be obeyed and ideally the lines joining L and P and R and P should
be as close to the characteristics as possible. This is sometimes not possible if for example there are
greatly differing velocities along a channel when equation (18) would force a very low time step at a
point where is it not required.

13
3.9 A second numerical solution method
While the first solution method shown used a finite difference solution to the actual St Venant equations,
this second method takes the characteristic form of the equations and solves these along the characteristic
paths. The solution to these is the same as the solution obtained from the original St Venant equations.

3.9.1 The method of characteristics on a rectangular grid
For this method we choose a regularly spaced grid arrangement as shown in Figure 7.

x
t
t
t+ t
P
L
R
O E W
x
t

Figure7 - A rectangular grid with characteristics through P.

This figure shows the advanced time solution point P, with the forward and backward characteristics
passing through it. These intersect with the initial time level line at points L (left) and R (right)
respectively.

Recall the characteristic equations:


( )
( ) j i g
dt
c v d

+ 2
for ( ) c v
dt
dx
+
(19)
and

( )
( ) j i g
dt
c v d

2
for ( ) c v
dt
dx

(20)

We will apply these along the forward and backward characteristics. First we must simplify the situation
by assuming that the characteristics are straight lines. This is true only if nothing is changing (steady
uniform flow) - a situation we are not very interested in. However if the time step is small the assumption
usually gives very good solutions. Figure 8 shows a close up of the grid and straight characteristics
associated with point P.

Secondly we must construct the characteristics. We know the slope when the forward characteristic
crosses the known level, t, it is given by ( )
L L
c v
dt
dx
+ and the slope when it passes through P is given by
( )
P P
c v
dt
dx
+ . The problem is that we dont know the values at either L or P! We must make an
approximation. The usual approximation to determine the slope of the characteristics is to use the point at
the same x position as P, but on the known time level i.e. the values at point O. The time step is t so the
slope of the forward characteristic is
14
( )
O O
L
c v
t
x
+

and the backward characteristic is

( )
O O
R
c v
t
x

x
P
W O E R
L
x
x
x
L
R
t
t+ t
t

Figure 8 Grid and associated characteristics for point P


Now we can apply equation (19) along the forward characteristic, L-P, to give:

( ) ( ) ( )
L L L P P
j i g t c v c v + + 2 2
(21)
And equation (20) along the backward characteristic, R-P, to give:

( ) ( ) ( )
R R R P P
j i g t c v c v 2 2
(22)

Solving equations (21) and (22) simultaneously gives the conditions at point P.

Adding equations (21) and (22) gives
( ) ( ) ( ) [ ]
R L R L
R L
P
j i j i
g t
c c
v v
v +

+ +
+

2 2

(23)
subtracting gives

( ) ( ) [ ]
R L
R L R L
P
j i j i
g t c c v v
c

+
+
+

4 2 4

(24)
These are the solution at point P.

To simplify matters we might have chosen to calculate the friction slope at the mid point O (as this is at
the same x position as P) which would give for the forward characteristic

( ) ( ) ( )
O L L P P
j i g t c v c v + + 2 2
(25)
for the backward characteristic
( ) ( ) ( )
O R R P P
j i g t c v c v 2 2
(26)
15
and the solution is
( ) ( )
O R L
R L
P
j i g t c c
v v
v + +
+

2

(27)

2 4
R L R L
P
c c v v
c
+
+

(24)
To calculate this solution we must obtain the values at the point L and R. They are on the initial time level
where we know all the solution, but we only know the solution at the node points (W, O and E).If we
know the position of L and R we can interpolate between these know values to get v
L
, c
L
, v
R
, and c
R
. they
are not at the node

A simple linear interpolation procedure is all that is necessary i.e. for v
L
( )
W O
L
O L
v v
x
x
v v

and for v
R
( )
E O
R
O R
v v
x
x
v v

3.9.2 The stability criteria

The stability criteria discussed in section 3.8 and shown in equation (17)

( ) c v
dt
dx
+ >

is valid for the method of characteristics on a rectangular grid. In this case, t dt and x dx so the
limit on the time step is

( ) c v
x
t
+

< .

This allows a time step twice that of the staggered grid method.

To implement this, as with the staggered grid method the maximum value of (v+c) for each point on the
know time level should be determined then a t calculated. Using a suitable safety factor the time step
would be


( )
max
9 . 0
c v
x
t
+

Choosing a time step like this ensures that the x
L
and x
R
are both less than x.

16
3.10 Boundary conditions
At any boundary there is only one characteristic present because the other is outside the channel. For an
upstream boundary we have the backward characteristic, for the downstream boundary we have the
forward characteristic - as shown in Figure 9.
x
P
W O L
x
L
t
t+ t
t
b
o
u
n
d
a
r
y
P
O E R
x
x
R
t
t+ t
t
b
o
u
n
d
a
r
y

Figure 9 Characteristics at upstream and downstream boundaries

Instead of the second characteristic equation we have a boundary condition equation. This will either
specify the depth (or stage) at P - e.g. from a stage curve possibly from historical records. Or it will
specify the discharge at P - e.g. from a flood hydrograph.

3.10.1 Depth specified - Upstream boundary

This is very simple to solve. The values of depth h
P
is given so c
P
readily calculated (from gh c ) then
the only solution for v
P
is required.

From the backward characteristic equation

( ) ( ) ( )
R R R P P
j i g t c v c v 2 2
so
( ) ( )
R P R R P
j i g t c c v v + 2

3.10.2 Depth specified - Downstream boundary

This is again very simple to solve. The values of depth h
P
is given so c
P
can be calculated and only the
solution for v
P
is required.

From the forward characteristic equation

( ) ( ) ( )
L L L P P
j i g t c v c v + + 2 2
so
( ) ( )
L P L L P
j i g t c c v v + + 2

17
3.10.3 Discharge specified - Upstream boundary

In this case the discharge Q
p
is given, but this must be converted to a value in terms of v
p
or c
p
, so the
shape of the channel section must be known to allow the area of flow to be calculated from depth. For
simplicity we will assume a rectangular channel so

P
P
P P
P P
v
g
c
b
v bh
Av Q
2

or

2
P
P
P
bc
g Q
v
Substitute this in to the backward characteristic gives

( ) ( )
R R R P
P
P
j i g t c v c
bc
g Q

,
_

2 2
2


which can be rearranged to give

( ) ( ) 0 2 2
2 3
+ + +
b
g Q
j i g t c v c c
P
R R R P P


This is usually solved using some iteration technique - e.g. Newton-Raphson.

3.10.4 Discharge specified - Downstream boundary

This is very similar to the upstream boundary except in this case we use the forward characteristic:

( ) ( )
L L L P
P
P
j i g t c v c
bc
g Q
+

,
_

+ 2 2
2


which can be rearranged to give
( ) ( ) 0 2 2
2 3
+ +
b
g Q
j i g t c v c c
P
L L L P P


And again, this can be solved using some iteration technique - e.g. Newton-Raphson.

3.11 Super-Critical flow
When the flow is supercritical the wave speed, c, is greater than the velocity of flow. Using the equation
for the backward characteristic to calculate x
L
we get

( )
L L L
c v t x

but c
L
> v
L
so x
L
is negative. The characteristic around point P then look like those in Figure 10
18
x
P
W O E
R L
x
x
x
L
R
t
t+ t
t

Figure 10: Characteristic in super-critical flow

The backward characteristic has now moved over to the left side and has a positive slope. It will never
move so far as to be to the left of the forward characteristic - but in very fast flow they can be very close.

The solution method is exactly the same as that set out in section 3.9. The only change is that v
R
and c
R

are now interpolated between points W and O (instead of E and O).

3.11.1 Boundary conditions in super-critical flow

At the upstream boundary in super-critical flow both the characteristics are outside the channel. In which
case two boundary condition equations must be given to specify both v and c at point P.

At the downstream boundary there are two characteristics available, as shown in Figure 11, and no
boundary condition is required.

(It is often a common error that a boundary condition is applied to a downstream boundary with super-
critical flow. In this situation you have three equations for two unknowns and the problem over specified
- and usually an error manifesting as instabilities will result.)

b
o
u
n
d
a
r
y
x
P
W O
R L
x
x
L
R
t
t+ t
t

Figure 11: Characteristics at a super-critical downstream boundary

The solution method for the supercritical downstream boundary is the same as the solution for either sub
or super-critical flow in the main stream (see section 3.9).