Steady Open Channel Test Problemswith Analytic Solutions: I Macdonald, M J Baines, N K Nichols and P G Samuels 3/95
Steady Open Channel Test Problemswith Analytic Solutions: I Macdonald, M J Baines, N K Nichols and P G Samuels 3/95
Steady Open Channel Test Problemswith Analytic Solutions: I Macdonald, M J Baines, N K Nichols and P G Samuels 3/95
DEPARTMENT OF MATHEMATICS
Steady Open Channel Test Problems with
Analytic Solutions1
Department of Mathematics
University of Reading
PO Box 220
Reading
RG6 2AX
United Kingdom
1 The work reported here forms part of the research programme of the Oxford/Reading
Institute of Computational Fluid Dynamics and was supported by the EPSRC and HR
Wallingford Ltd.
2 HR Wallingford Ltd, Howbery Park, Wallingford, Oxon OX10 8BA, UK.
Abstract: In many areas of computational uid dynamics there are
benchmark test problems which have known analytic solutions. For the case
of steady ow in an open channel there is a need for such test problems, so as
to give a measure of the performance of numerical methods. This is especially
true for problems where the ow changes between subcritical and supercrit-
ical. In this report a method is described that allows the construction of a
wide range of such test problems, where in each case the exact solution to the
full steady Saint Venant equation is known. In particular the report discusses
how to construct problems with solutions having a hydraulic jump. Exam-
ple test problems are given along with their analytic solution. These can be
readily used by modellers to test their own particular software. For some of
the example test problems, results from a commercial software package are
compared with the exact solutions, thus demonstrating how the method is a
valuable tool for validating both steady and unsteady ow solvers.
1 Introduction
The use of numerical methods to compute the water surface prole and dis-
charge, for both unsteady and steady open channel ows, is now very common
in civil engineering hydraulics. The Saint Venant equations, rst formulated
by De St.Venant (1871), are almost always used to model the ow. Deriva-
tions of these equations can be found, for example, in Cunge, Holly and
Verwey (1981) and in Yen (1973). In this report attention is restricted to the
important case of steady ow. Apart from the intrinsic interest of the steady
ow problem, a reason why the computation of steady ows is so signicant
is that they are often required as initial data for unsteady simulations. Un-
der steady conditions the Saint Venant equations reduce to a single Ordinary
Dierential Equation (ODE), which, if the ow is either wholly subcritical or
wholly supercritical, can be eciently integrated using a high accuracy nu-
merical ODE solver. The situation is much more dicult if the ow is mixed
sub-supercritical, having hydraulic jumps or critical sections. Algorithms do
exist which can handle such ows, for example the method of Humpidge and
Moss (1971) which combines the ODE integration with the location of con-
trol points and the tting of hydraulic jumps. Another method is to apply
an unsteady solver and proceed forward in time until all the transients in
the ow have decayed and the ow has reached a steady state. In general
this approach can be very time consuming for trans-critical ows, although
the work of Garcia-Navarro, Priestley and Alcrudo (1994) gives a numerical
method with improved eciency.
Recently in MacDonald, Baines and Nichols (1994) a new approach has
1
been developed which uses shock capturing methods (see Engquist and Osher
(1981)) such as those that have been previously applied to the unsteady
system of equations, but are there applied to the single steady equation.
With the recent emergence of these new steady solvers, some way is re-
quired to compare the performance of these methods with existing schemes.
In many areas of computational uid dynamics there are standard test prob-
lems, for which the exact solution is known. Using some measure of the
dierence between the numerical solution and the exact solution, the perfor-
mance of a particular numerical scheme can be evaluated. Important features
of the solution, for example shocks, can be compared with the exact solu-
tion. This allows the strengths and weaknesses of one scheme over another
to be judged. If a particular scheme gives acceptable performance over a
wide range of such test problems, then, even though there may be no theory
guaranteeing good performance for a practical problem, there will be more
condence in the method. Altogether, the existence of test problems with
known solutions is extremely desirable.
For the steady open channel problem, because of the non-linear nature
of the dierential equation, particularly the friction term, even the simplest
problem of ow in a uniform rectangular channel with zero bed slope cannot
be solved analytically. To obtain a solvable problem it is necessary to assume
zero friction or at least to simplify the friction term signicantly. Features of
the problem will then be lost. For this reason, until now, the performance of
methods for steady computation have mostly been judged only qualitatively.
This report presents a simple method for constructing test problems with
known analytic solutions to the full steady Saint-Venant equation. The
method is an \inverse method" in that some hypothetical depth prole is
chosen and the bed slope that makes this prole an actual solution of the
steady equation is then found. The method can be used to construct test
problems with almost any desired features, including hydraulic jumps. Hence
these test problems can be used to compare the numerical results, for any
algorithm, with an exact solution. The method is also useful for evaluat-
ing unsteady solvers, since, if an unsteady model is given steady boundary
conditions, the limiting steady solution can be compared with the analytic
steady solution. The method presented in this report ts in well with the
validation documentation initiative of the European hydraulics laboratories
(see Dee (1993)), since it enables the creation of benchmark test problems
which can be used as a standard measure for the performance of commercial
software packages.
2
2 Background and Theory
2.1 The Equations
The Saint Venant equations are obtained from the principles of mass and
momentum balance and are given by
@A + @Q = q (1)
@t @x
!
@Q + @ Q2 + gA @y ; gA(S ; S ) = 0: (2)
@t @x A @x 0 f
Here x is the distance along the channel, t is time, y(x t) is the depth, Q(x t)
is the discharge, T (x y) is the free surface width, A(x y) is the wetted area,
S0(x) is the bed slope, Sf (x y Q) is the friction slope, q is the lateral inow
per unit length, is the momentum coecient and g is the acceleration due
to gravity. The bed slope, S0, is given by
dz
S0 = ; dx (3)
where z(x) is the bed level, the elevation of the bed above some horizontal
datum. The friction slope, Sf , is given by
Sf = QKjQ2 j
where K (x y) is the conveyance, which can be taken to be one of the many
commonly used formulae (see Chow (1959)). Here Manning's equation, for
which
K = nPA5=3
2=3
is used, where P (x y) is the wetted perimeter and n is the Manning friction
coecient. However, the method described in this work does not rely on this
choice of conveyance.
In this report only the steady ow problem is considered, so it is assumed
that y = y(x) and Q = Q(x). Under these steady conditions equations (1)
and (2) reduce to
dQ = q (4)
!
dx
d Q2 + gA dy ; gA(S ; S ) = 0: (5)
dx A dx 0 f
3
Although not strictly necessary, the following simplifying assumptions are
made. There is uniform velocity over each cross-section, which implies that
the momentum coecient and the energy coecient are both unity. Also
the lateral inow q is assumed to be zero. With zero lateral inow, equation
(4) becomes trivial with solution Q constant. Since the discharge can have
no jumps it must be constant throughout the entire channel reach. From
now on x will be measured in the direction of this constant discharge and
hence Q > 0.
Dierentiating the momentum term in equation (5) and dividing through
by gA then yields the equation
!
dy ; Q2 @A ; S + S = 0:
1 ; QgAT3 dx
2
gA3 @x 0 f (6)
It is convenient to re-write this equation as
S0(x) = f1(x y(x))y (x) + f2(x y(x))
0
(7)
where
f1 = 1 ; QgAT3 = 1 ; Fr2
2
(8)
and
Q 2 n2 P 4=3 Q 2 @A
f2 = A10=3 ; gA3 @x : (9)
Fr is the Froude number. The assumptions on , and q can be relaxed,
and equation (6) can still be written in the form (7) for dierent functions
f1 and f2.
E = 2Q
2
gA2 + y:
For a rectangular channel the jump condition explicitly yields the required
value of y^R(x?) for any value of y^L(x?) and vice versa. For other cross-
sections it is necessary to solve the jump condition numerically. For many
cross-sections, where there is a unique critical depth, the energy condition
is satised if and only if the jump is from supercritical to subcritical. For
cases where the situation is not clear, it is easy to check the condition for
each individual jump. The bed slope of the channel can now be dened in a
piecewise manner by
8
< S0L (x) 0 x x?
>
S0(x) = >: (12)
S0R(x) x? < x L
where
S0L(x) = f1(x y^L(x))^yL(x) + f2(x y^L(x))
0
For this bed slope, y^ satises the dierential equation (7) everywhere except
at the jump. In general the bed slope is discontinuous at x = x?, i.e.
S0L(x?) 6= S0R(x?):
At rst sight this discontinuity might seem perfectly acceptable, since many
valid test problems have such a feature. However, in general the hydraulic
jump does not occur at the same position as the bed slope discontinuity and
so we would like to construct problems where the jump does not coincide
with a bed slope discontinuity. This can be achieved in the following way.
Having chosen values for y^L(x?) and y^R(x?), choose values for y^L(x?) and
0
(13)
= f1(x y^R(x ))^yR(x ) + f2(x y^R(x )) = S0R(x ):
? ? 0 ? ? ? ?
6
If the functions are smooth enough, equation (12) can be dierentiated to
nd a linear relationship between y^L(x?) and y^R(x?) in order to make the
00 00
chosen with three free parameters. Values for these parameters are chosen so
as to give the required values for y^R, y^R and y^R at x = x?. In general more
0 00
than three parameters are required so that there is some control over the
behaviour of y^R. In particular y^R must always remain positive. The obvious
choice for y^R is a cubic or higher order polynomial in (x ; x?). The examples
in MacDonald (1994) use this form. The diculty with polynomials is that
they can be highly oscillatory and hence dicult to control and keep positive,
a diculty which increases the greater the distance y^R is required to cover.
There are many other possible forms for y^R. The examples in this report
use sums of exponential functions (see Appendix III). The magnitudes of y^R 0
and y^R at the jump are often required to be high. The choice of exponentials
00
allows the restriction of these high magnitudes to the locality of the jump.
Judicious choice of these exponentials allows control over the behaviour of y^R
away from the jump. Even with a reasonable form for y^R, it can take much
trial and error to achieve a specic behaviour.
3 Examples
This section describes specic examples of test problems constructed using
the method given in this report. The full details of the test problems are
given in Appendix III. In Examples 1-4, a 1km long, 10m wide rectangular
channel with a discharge of 20m3/s is chosen. In Example 1 the hypothetical
depth function is the subcritical hump shown in Figure 1. The resulting bed
slope is also shown in the same gure. Figure 2 shows the resulting bed level
and free surface level.
7
In Example 2 the hypothetical depth function is an inverted hump with
the ow being wholly supercritical. Figure 3 shows the resulting bed level
and free surface level.
In Example 3 the hypothetical depth function is chosen to be subcritical
at inow and to change smoothly to supercritical halfway along the chan-
nel, remaining supercritical for the rest of the channel. Figure 4 shows the
resulting bed level and free surface level.
In Example 4 the hypothetical depth function is supercritical at inow,
changes via a hydraulic jump to subcritical halfway along the channel and
remains subcritical. The hypothetical depth function is shown in Figure 5 as
well as the resulting bed slope. Figure 6 shows the resulting bed level and
free surface level.
pIn Example 5 a 5km long 3trapezoidal channel (T = 10 + 4y, P = 10 +
2y 5) with a discharge of 20m /s is chosen. The hypothetical depth function
is the subcritical sine function shown in Figure 8. Figure 7 shows the resulting
bed level and free surface level.
pIn Example 6 a 1km long trapezoidal channel (T = 10 + 2y, P = 10 +
2y 2) with a discharge of 20m3/s is chosen. In this case the hypothetical
depth function, shown in Figure 10, is subcritical at inow, changes smoothly
to supercritical and then returns to subcritical via a hydraulic jump. Figure 9
shows the resulting bed level and free surface level.
9
For two of the examples given, numerical results from a commercial steady
open-channel solver have been compared with the exact solution, demonstrat-
ing how it is possible to judge the performance of the solver.
The method described in this report is a valuable tool for developing,
validating or comparing steady open-channel solvers. The method can also
be used to test the performance of unsteady models as the solution tends to
a steady state.
Acknowledgements
This work was supported by the Engineering and Physical Sciences Research
Council, UK. and HR Wallingford Ltd, UK.
Appendix I: References
V T Chow (1959). \Open Channel Hydraulics". McGraw-Hill, London,
England.
J A Cunge, F M Holly, and A Verwey (1980). \Practical Aspects of Compu-
tational River Hydraulics." Pitman, London, England.
D Dee (1983). \Standard Validation Document:- Denition and Guidelines."
Delft Hydraulics, Netherlands.
B Engquist and S Osher (1981). \One Sided Dierence Approximations for
Nonlinear Conservation Laws." Mathematics of Computation, 36, 321-351.
H B Humpidge and W D Moss (1971). \The Development of a Compre-
hensive Computer Program for the Calculation of Flow Proles in Open
Channels." Proc. of the Inst. of Civ. Engrs., 50, 49-64.
P Garcia-Navarro, A Priestley and F Alcrudo (1994). \An Implicit Method
for Water Flow Modelling in Channels and Pipes." Journal of Hydraulic
Research, 32(5), 721-742.
I MacDonald (1994). \Test Problems with Analytic Solutions for Steady
Open Channel Flow." Numerical Analysis Report 6/94, Department of
Mathematics, University of Reading, UK.
10
I MacDonald, M J Baines and N K Nichols (1994). \Analysis and Computa-
tion of Steady Open Channel Flow using a Singular Perturbation Problem."
Numerical Analysis Report 7/94, Department of Mathematics, University of
Reading, UK.
I MacDonald, M J Baines, N K Nichols and P G Samuels (1995). \Com-
parison of some Steady State Saint-Venant Solvers for some Test Problems
with Analytic Solutions." Numerical Analysis Report (2/95), Department of
Mathematics, University of Reading, UK.
R K Price and Samuels P G (1980). \A Computational Hydraulic Model for
Rivers." Proc. of the Inst. of Civ. Engrs. 69(2), 87-96.
P G Samuels and M P Gray (1982). \The FLUCOMP River Model Package
- An Engineers Guide." Report No. EX 999, HR Wallingford Ltd, UK.
S Wolfram (1988). \Mathematica, A System for Doing Mathematics by
Computer." Addison-Wesley, New York.
De Saint Venant (1871). \Theorie du mouvement non-permanent des eaux
avec application aux crues des rivieres et a l'introduction des marees dans
leur lit." Acad. Sci. Comptes rendus, 73, 148-154, 237-240.
B C Yen (1973). \Open-Channel Flow Equations Revisited." Journal of The
Engineering Mechanics Division, ASCE, 99(5), 979-1009.
Fr Froude Number.
x? Position of hydraulic jump (m).
y^L, y^R Hypothetical depth functions for left and right of jump (m).
S0L, S0R Bed slope functions for left and right of jump.
ymax Maximum depth for hypothetical depth functions (m).
B Width for rectangular channels (m).
a1:::a3 Coecients in depth functions.
k Integer index.
Dummy integration variable.
Example 1
A 1km long rectangular channel of width 10m has a discharge of 20m3/s.
The ow is subcritical at inow and is subcritical at outow with depth
0.748409m. The Manning roughness coecient for the channel is 0.03 and
the bed slope is given by
!
4 y
(2^ ( x) + 10) 4=3
S0(x) = 1 ; gy^(x)3 y^ (x) + 0:36 (10^y (x))10=3
0
where !1=3 x 2 !!
y^(x) = g4 1 + 21 exp ;16 1000 ; 12
12
and
!1=3 !
y^ (x) = ; g4
0
2 x ; 1 exp ;16 x ; 1 2 :
125 1000 2 1000 2
The solution to this problem is given by y y^ and is shown in Figures 1 and
2.
Example 2
A 1km long rectangular channel of width 10m has a discharge of 20m3/s.
The ow is supercritical at inow with depth 0.741599m and is supercritical
at outow. The Manning roughness coecient for the channel is 0.02 and
the bed slope is given by
!
4 y
(2^ ( x) + 10)
S0(x) = 1 ; gy^(x)3 y^ (x) + 0:16 (10^y (x))10=3
0
4=3
where !1=3 x 2 !!
y^(x) = g 4 1
1 ; 5 exp ;36 1000 ; 2 1
and !1=3 !
y^ (x) = 4g
0
9 x ; 1 exp ;36 x ; 1 2 :
625 1000 2 1000 2
The solution to this problem is given by y y^ and is shown in Figure 3.
Example 3
A 1km long rectangular channel of width 10m has a discharge of 20m3/s.
The ow is subcritical at inow and is supercritical at outow. The Manning
roughness coecient for the channel is 0.02 and the bed slope is given by
!
y (x) + 10)4=3
S0(x) = 1 ; gy^(4x)3 y^ (x) + 0:16 (2^(10^
0
y (x))10=3
where
8
>
> 4 1=3
1 ; 13 tanh 3 x
; 12 0 x 500
< g 1000
y^(x) = >
>
: 4 1=3
1 ; 16 tanh 6 x
; 12 500 < x 1000
g 1000
13
and
8
>
<; ; 12 0 x 500
4 1=3 1
> g 1000 sech2 3 x
1000
y^ (x) = >
0
>
:; 4 1=3 1
; 12 500 < x 1000:
sech2 6 x
g 1000 1000
Example 4
A 1km long rectangular channel of width 10m has a discharge of 20m3/s.
The ow is supercritical at inow with depth 0.543853m and is subcritical at
outow with depth 1.334899m. The Manning roughness coecient for the
channel is 0.02 and the bed slope is given by
!
y (x) + 10)4=3
S0(x) = 1 ; gy^(4x)3 y^ (x) + 0:16 (2^(10^
0
y (x))10=3
where
8
>
> 4 1=3 9
; 16 exp ;x
0 x 500
>
<
g 10 250
y^(x) = > P3
k=1 ak exp ;20k ; 12 +
4 1=3
>
> g 1 + x
1000
: exp
1000 ; 1 500 < x 1000
4 x
5
and
8
>
> 4 1=3 1
exp ;x
0 x 500
>
<
g 1500 250
y^ (x) = >
; 501 P3k=1 kak exp ;20k
0
>
>
4 1=3
g
x
1000 ; 12 +
: 1
exp 1000x
;1 500 < x 1000
1250
Example 5
p
A 5km long trapezoidal channel (T = 10+4y, P = 10+2y 5) has a discharge
of 20m3/s. The ow is subcritical at inow and is subcritical at outow with
14
depth 1.125m The Manning roughness coecient for the channel is 0.03 and
the bed slope is given by
! p 4=3
y ( x)) 0
(10 + y
2^ (x )
S0(x) = 1 ; g(10 + 2^y(x))3y^(x)3 y^ (x) + 0:36 (10 + 2^y(x))10=3y^5)(x)10=3
400(10 + 4^
where
y^(x) = 89 + 14 sin 500x
and x
y^ (x) = 2000 cos 500 :
0
Example 6
p
A 1km long trapezoidal channel (T = 10+2y, P = 10+2y 2) has a discharge
of 20m3/s. The ow is subcritical at inow and is subcritical at outow with
depth 1.349963m The Manning roughness coecient for the channel is 0.02
and the bed slope is given by
! p 4=3
400(10 + y
2^ (x )) (10 + y
2^ ( x
S0(x) = 1 ; g(10 + y^(x))3y^(x)3 y^ (x) + 0:16 (10 + y^(x))10=3y^2)
)
(x)10=3
0
where
8
> 0:723449 1 ; tanh 1000 ; 10
> 0 x 300
x 3
>
>
>
>
< 0:723449 1 ; 16 tanh 6 x
; 103 300 < x 600
y^(x) = > 1000
>
> P3 a exp ;20k
>
>
3
4 + k=1 k
x
1000 ; 35 +
:3 expx
;1 600 < x 1000
5 1000
and
8
>
> ;0:723449 10 3 sech2 1000
; x
; 103 0 x 300
>
>
>
< ;0:723449 10 3 sech2 6 1000
y^ (x) = >
0
; x
; 103 300 < x 600
>
> 1 P3
>
> ; 50 k=1 kak exp ;20k 1000 ; 5
x 3
+
: 3 exp 1000 ; 1
x
600 < x 1000
5000
16
1.2
1.0
0.8
0.6
0.4
Depth /m
0.2 Bed Slope x 100
Critical Depth /m
0.0
0.0 200.0 400.0 600.0 800.0 1000.0
x /m
8.0
6.0
Height /m
4.0
17
8.0
6.0
Height /m
4.0
8.0
6.0
Height /m
4.0
18
1.5
Depth /m
Bed Slope x 50
Critical Depth /m
1.0
0.5
0.0
0.0 200.0 400.0 600.0 800.0 1000.0
x /m
8.0
Bed Level
Free Surface Level
6.0
Critical Level
Height /m
4.0
2.0
0.0
0.0 200.0 400.0 600.0 800.0 1000.0
x /m
19
15.0
Height /m 10.0
0.0
0.0 1000.0 2000.0 3000.0 4000.0 5000.0
x /m
1.40
1.20
Depth /m
1.00
FLUCOMP
0.80 Exact Depth
0.0 1000.0 2000.0 3000.0 4000.0 5000.0
x /m
20
5.0
Bed Level
4.0
Free Surface Level
Critical Level
Height /m 3.0
2.0
1.0
0.0
0.0 200.0 400.0 600.0 800.0 1000.0
x /m
1.5
Depth /m
Bed Slope x 50
Critical Depth /m
1.0 FLUCOMP
0.5
0.0
0.0 200.0 400.0 600.0 800.0 1000.0
x /m
21