Cac Phan Mem Quan Ly Tuoi PDF
Cac Phan Mem Quan Ly Tuoi PDF
Cac Phan Mem Quan Ly Tuoi PDF
Institutional Repository
Database-driven hydraulic
simulation of canal
irrigation networks using
object-oriented
high-resolution methods
This item was submitted to Loughborough University's Institutional Repository
by the/an author.
Additional Information:
040317807X
~.
. -. . ·.·. f~glt~.O
Lo
, ·Unrverslty' .•. I'OUgh
.
Sachin Shende
Loughborough University
Doctor of Philosophy
... ,'
April 2006
: "
·i,
,
~ Loughborough
Unh'crsity
P;lki"~lon library
Completing this doctoral work has been a wonderful, challenging and exciting experi-
ence. I have been very privileged to have a Supervisor like Dr Cecil Scott. His wide
mathematical interests, strong object-oriented programming skills, ability to think in a
structured way and organised approach has greatly influenced this work. I have learned
a great deal of Java programming from him and I have been stimulated and excited by
his constant flow of good ideas. He pushed me to delve deeper into all aspects of my
research and helped me to truly understand the object-oriented concepts and numerics
of the software I was developing.
Ian Smout had been undoubtedly the most supportive Supervisor anyone could ask
for. It was a wonderful experience working with him. He influenced the shape and
direction of this research with his helpful comments and gave me a right push in the
forward direction when I needed it.
I would like to express my deepest gratitude to Dr Sunil Gorantiwar who gave me a
prudent advice of pursuing the doctoral research programme at Lougborough University
and actually helped me settle here during flrst couple of months. His guidance as an
External Supervisor for this research has been invaluable.
Throughout my three years, I was supported by the funding from Water Engineering
Development Centre (WEDC), Loughborough University. This research would have
never got off the ground without WEDC's funding.
I am grateful to my Director of Research, Dr M. Sohail for his support, guidance
and inspiration. A special thanks to him for arranging a JBuilder Enterprise license
for this research.
The other Faculty member of the Department who had a strong influence on me is
Prof. Koji Shiono. His short, witty lectures on hydraulics were enlightening, especially
ACKNOWLEDGMENTS
interesting were his talks on physics behind various flow phenomenon. I am grateful to
him for his terse advice and for keeping my research on the track.
I must thank our Postgraduate Research Administrator, Helen, for all her support
she provided for the last three years. She is surely a kindest and most helpful person
in our Department.
I am grateful to people from outside this University who have helped and provided
advice when required, namely, Prof. Forrest Holly, University of Iowa, Dr Albert
Clemmens, USWCL and Dr Richard Crowder, Halcrow Group Ltd. I am indebted
to the team from DEFRA/EA which worked on the benchmarking of hydraulic river
modelling software project. The datasets and documentation they provided as a part
of this project has been extremely useful. I am especially grateful to Prof. Randall
Leveque, University of Washington, who developed and distributed the CLAWPACK
software for the benefit of a wider research community.
Typesetting of this thesis has been done in Jb.'IE;X 20 using MiKTeX, TeXnicCenter
and a modified version of CUED Thesis class file. This has been a very challenging,
fascinating, at times frustrating, and often sweaty affair. In doing so, CTAN website,
Jb.'IE;X primer by the Indian TUG, TeX FAQ on http://www.tex.ac.uk/ and discus-
sion on comp.text.tex had been extremely useful. This research has benefitted from
quite a few open source softwares. I would like to recognise and support all those
people who developed applications and systems like CLAWPACK, CVS, DBDesigner,
GIMP, Jb.'IE;X20, MiKTeX, TeXnicCenter, MySQL database, Poseidon and made them
available for free.
Thanks are also due to my research friends Atmanand, Keyur and Deepak for
their valuable discussion and continued support. A special thank you to Harshal for
all his help and support during last three years. Mukesh deserves a big thank you
for being a willing resource for all matters and for his guidance and encouragement
throughout. This acknowledgment will not be complete without mention of Mukund
and Mohammed. The research room where I worked has fostered an open, friendly and
collaborative research group in the Department. This research has certainly benefited
from the cultural and academic diversity of the place.
ii
ACKNOWLEDGMENTS
Finally, and of course far more importantly, the family. I do not know how to
express my feelings of gratitude towards my Dad for providing me all the opportunities
and been there for me when I needed him. He always wanted to see this work but at the
end he could not. Thanks to Mom and bro for their support. A huge debt of gratitude
to Sudeshna for her unending support and understanding during last three difficult
years. She has seen my best and my worst, and provided support, hugs, and taken me
places I never imagined. Also a dollop of thanks is due to my in-laws especially my
sister in-law, Sulogna, who has supported me in countless ways. Without their support
this research would not have been possible.
iii
Abstract
Canal hydraulic models can be used to understand the hydraulic behaviour of large
and complex irrigation networks at low cost. A number of computational hydraulic
models were developed and tested in the early 1970s and late 80s. Most were developed
using finite difference schemes and procedural programming languages. In spite of
the importance of these models, little progress was made on improving the numerical
algorithms behind them. Software development efforts were focused more on developing
the user interface rather than the core algorithm.
This research develops a database-driven, object-oriented hydraulic simulation model
for canal irrigation networks using modern high-resolution shock capturing techniques
that are capable of handling variety of flow situations which includes trans-critical flow,
shock propagation, flows through gated structures and channel networks. The technol-
ogy platforms were carefully selected by taking into account a multi-user support and
possible migration of the new software to a web-based one which integrates a Java-based
object-oriented model with a relational database management system that is used to
store network configuration and simulation parameters.
The developed software is tested using a benchmark test suite formulated jointly by
the Department for Environment, Food and Rural Affairs (DEFRA) and the Environ-
ment Agency (EA). A total of eight tests (seven of them adapted from the DEFRAjEA
benchmark suite) were run and results compiled. The developed software has outper-
formed ISIS, REC-RAS and MIKE 11 in three of the benchmark tests and equally well
for the other four. The outcome of this research is therefore a new category in hydraulic
simulation software that uses modern shock-capturing methods fully integrated with a
configurational relational database that has been fully evaluated and tested.
Contents
Contents xi
Nomenclature xxviii
1 Introduction 1
1.1 Genesis . 2
1.1.1 Rationale 4
1.1.2 Aims and Objectives . 5
1.2 Thesis Structure ... . . . . 6
v
CONTENTS
2.6.10 Hivel-2D . . . . 37
2.6.11 ISIS .. :37
2.6.12 Mike 11 38
2.6.13 MODIS 38
2.6.14 NOAH . . . . . 39
2.6.15 SIC .. ,10
2.6.16 Sobek . 41
2.6.17 TELEMAC . . . . . 41
2.6.18 UNET . . . . . 42
2.6.19 USM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
2.6.20 Summary of Hydraulic Simulation Softwares . . . . . . . . . . . 43
2.7 Testing and Validation of the Software . . . . • . . . . . . . .. 44
vi
CONTENTS
4 CanalFlow - Theory 66
4.1 Introduction.... 66
4.2 One-Dimensional Unsteady Open Channel Flow. 67
4.2.1 Integral Form of the Equations .. 67
4.2.2 Differential Form of the Equations 73
4.2.3 Matrix Form of the Equations ... 76
4.2.4 Eigenstructure of the Governing Equations 77
4.2.5 Characteristics Theory . 79
4.3 Rlemann Problem . . . . . . . 80
4.4 Finite Volume Method - Formulation. 84
4.4.1 Godunov Type Scheme . . . . 86
4.4.1.1 The Wave-Propagation Form of Godunov's Method 88
4.4.2 Exact Riemann Solver 91
4.5 High Resolution Methods 95
4.5.1 Limiters . . . . . . 96
4.6 Initial and Boundary Conditions 97
4.6.1 Initial Conditions ... 97
4.6.2 Boundary Conditions 98
4.6.2.1 Method of Characteristics Approach. 98
4.6.3 External Boundaries . . . . . . . 10,5
4.6.3.1 Transmissive Boundary 105
4.6.3.2 Wall Boundary . . . . . 105
4.6.3.3 Discharge Hydrograph . 106
vii
CONTENTS
viii
CONTENTS
ix
CONTENTS
x
CONTENTS
Glossary 224
References 288
xi
List of Tables
5.1 Different software tools used during various phases of the research. 119
6.1 The benchmark test suite used to test CanalFlow. Listed are the eight
test cases classified under three broad categories. . . . . . . . . . . . . . 142
6.2 Definition and type of boundary condition used for the sub-parts N.1.1
to N.L6 in the sub-critical, super-critical and transitional Hows test case. 144
6.3 Boundary condition details for sub-parts N.L1 to N.L5 in the sub-
critical, super-critical and transitional Hows test case. . . . . . . . . . . 145
6.4 Boundary condition data for parts N.2.1 and N.2.2 of the triangular
channel test case. . . . . . . . . . . . . . . . . . . . . . . . . . . 147
6.5 Channel geometry details for the looped system (C.1) test case. 1f>4
6.6 Broad crested weir parameters used by CanalFlow. . . . . . . . 156
6.7 Details of the rectangular gates used in the irrigation network (C.3) test
case. Cd is the coefficient of discharge used for the free and submerged
How condition. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159
xii
LIST OF TABLES
6.8 Gate operation schedule for the rectangular gates used in the irrigation
network test case (Refer Figure 6.8). The gate number 1, 2, 5 and 8
are opened in first hour, kept at the same opening during second hour
and then closed in third hour. Between 3 and 6 hours they are kept at
the same gate opening. Gates 3, 4, 6 and 7 are kept at constant gate
opening and are not operated at all. .. . . . . . . . . . . . . . . . . . . 159
7.1 Network parameters set for the sub-critical, super-critical and transi-
tional flows (N.l) test case. . . . . . . . . . . . . . . . . . . . . . . . . . 170
7.2 Triangular channel test (N.2): Comparison of water levels generated by
CanalFlow, ISIS, HEC-RAS and MIKE 11 run in a resistance radius
(RR) and hydraulic radius (HR) mode for the sub-critical flow test case
(N.2.1). The results for softwares other than CanalFlow presented in
this table and Table 7.3 are adapted from Crowder et al. (2004d). . . . . 180
7.3 Triangular channel test (N.2): Comparison of water levels generated by
CanalFlow and other softwares used in the EA benchmark suite for the
super-critical flow test case (N.2.2) . . . . . . . . . . . . . . . . . . . . . . 180
7.4 Network parameters set for the Monoclinal wave (N.4) test case. The
limiters used are Minmod (MM), Monotonised centered (MC), Superbee
(SB) and van Leer (VL) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 189
7.5 Minimum and maximum errors between CanalFlow and analytical so-
lution using various limiters. The limiters used are Minmod (MM),
Monotonised centered (MC), Superbee (SB) and van Leer (VL). The
data clearly shows that the van Leer limiter has a minimum range be-
tween minimum and maximum error while Superbee has got maximum
error on the negative side and Minmod has got maximum error on the
positive side with an exception at 10 Hr where Monotonised centered
has got maximum positive error. The last row in the table shows MIKE
11 error for the same test case. ... . . . . . . . . . . . . . . . . . . . . 198
xiii
LIST OF TABLES
7.6 Looped system (C.1) test: Flow and stage results for the QS1 boundary
condition. Flow is expressed in m 3 / sec while chainage and stage is in
m. HEC-RAS, ISIS, MIKE 11 results are adapted from Crowder et al.
(2004c). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202
7.7 Looped system (C.1) test: Flow and stage results for the QS2 boundary
condition. Flow is expressed in m 3 / sec while chainage and stage is in
m. HEC-RAS, ISIS, MIKE 11 results are adapted from Crowder et al.
(2004c) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203
7.8 Weir flow (C.2) test case: Steady state results for free and submerged
flow condition. The output shows that the head loss in SSl is equal to
the ISIS and for SS2 it is closer to ISIS than MIKE 11 and HEC RAS. . 205
xiv
LIST OF TABLES
B.1 Initial data for the sub-critical (N.l.1) sub-part in the sub-critical, super-
critical and transitional flows test case. Chainage, bed, depth are given
in m while discharge is in cumec. . . . . . . . . . . . . . . . . . . . . . . 238
B.2 Initial data for the super-critical flow (N.l.2) sub-part in the sub-critical,
super-critical and transitional flows test case. Chainage, bed, depth are
given in m while discharge is in cumec. . . . . . . . . . . . . . . . . . . . 240
B.3 Initial data for the super-critical to sub-critical flow (N.l.3) sub-part in
the sub-critical, super-critical and transitional flows test case. Chainage,
bed, depth are given in m while discharge is in cumec. . . . . . . . . . . 242
BA Initial data for the sub-critical to super-critical to sub-critical flow (N.l.4)
sub-part in the sub-critical, super-critical and transitional flows test case.
Chainage, bed, depth are given in m while discharge is in cumec. . . . . 244
B.5 Initial data for the super-critical to sub-critical to super-critical (N.l.5)
flow sub-part in the sub-critical, super-critical and transitional flows test
case. Chainage, bed, depth are given in m while discharge is in cumec. . 246
B.6 Initial data for the sub-critical and super-critical flow test in triangular
channels (N .2.1 and N .2.2). Chainage, bed, depth are given in m while
discharge is in cumec. . . . . . . . . . . . . 248
B.7 Initial data for the Ippen wave test (N.3). Chainage, bed, depth are
given in m while discharge is in cumec. . . . 249
B.8 Initial data for Monoclinal rising wave test (NA). Chainage, bed, depth
are given in m while discharge is in cumec. . 252
B.9 Initial data for the Looped System test(C.1). dx (distance step), bed,
depth, bottom are in m, discharge is in cumec and side slope is expressed
as H:V. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 255
xv
LIST OF TABLES
B.lD Channel geometry details used in the irrigation network (C.3) test case.
The length, bottom and depth are in m, slope is expressed as 1 in column
value, side slope is specified as horizontal to vertical while the discharge
is in cumec. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 256
B.n Initial data for test R.I - Dambreak in converging and diverging chan-
nel section. The chainage, bed, depth and bottom are given in m and
discharge is in cumec. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 258
xvi
List of Figures
2.1 Plan view of a typical irrigation network showing a reservoir, main canal,
control structures, distribution network and a field application system. 10
xvii
LIST OF FIGURES
xviii
LIST OF FIGURES
4.11 Illustration of the computational domain [a, bJ of the length L and cell
average values U. F l / 2 is the flux entering the domain at a, upstream
side, and FN+1/2 is the flux leaving the domain at b, downstream side.
Also shown is the right and left going fluctuations A + and A - , respectively. 99
4.12 Illustration of the external boundaries at channel ends, (A) C_ charac-
teristic is used to solve the boundary condition at the upstream end and
(B) C+ is used at the downstream end . . . . . . . . . . . . . . . . . . . . 100
4.13 Illustration of an irrigation network showing nodes, channels, control
structures, junctions and external boundary points. Internal boundary
conditions are solved at the junctions. Shown are the junctions; J-1 with
three channels and a structure, J-2 as a three channel open junction, J-3
with two channels and a structure and J-5 as a two channel open junction.108
4.14 Junctions J-l, J-2, J-3 and J-5 as seen in the Figure 4.13. . . . . . .. 109
4.15 Channels connected through a structure and h2 < ho < hl (free flow). III
4.16 Channels connected through a structure and ha < h2 < hl (submerged
flow) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
5.1 Components of canalfiow package. Shown are various classes and their
relations with other classes. For example in a network package, rue-
mann1DSolver is a sub-class of a Channel class additionally it imple-
ments OneDSolver interface. Also, Network consists of many Channel,
Node and Reservoir objects while a Node might/might not have a Junction.118
5.2 Illustration of a ChannelStructure interface. All classes in canalflow.channelstructures
implement this interface and as seen in the figure this interface is also
used by Channel and ruemann1DSolver class. It defines necessary method
signatures for geometric calculations of a channel section. . . . . . . . . 121
5.3 T'rapezoidalStructure class diagram implements ChannelStructure inter-
face and implements all methods related to geometric calculations of
trapezoidal, rectangular and triangular channel sections. . . . . . . . . . 122
xix
LIST OF FIGURES
5.4 DBConnection class diagram. This class is used by other classes for
database related functionalities. . . . . . . . . . . . . . . . . . . . . . . . 123
5.5 Illustration of a Canal class diagram. Canal has one or more Channel
objects and is a part of a network package. . . . . . . . . . . . . . . . . 123
5.6 OneDSolver interface. The new solvers developed in future could imple-
ment this interface. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
5.7 RiemannlDSolver, a sub-class of Channel class implements OneDSolver
interface. This class implements exact Riemann solver and high-resolution
methods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
5.8 Illustration of a Network class diagram. This is a central class which
initiate all processes pertaining to building network, simulation process
and time keeping operations. The concept of multithreading is used
to run the simulation and visualisation processes simultaneously as two
separate threads. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
5.9 Shown is the Node class diagram. Node objects are building blocks of a
network. Important role of a Node class is to initiate right kind of bound-
ary condition based on the available data but the actual application of
boundary condition is done by the Channel class. 128
5.10 Channel class diagram. Continued... . . . . . . . 130
5.ll Junction class diagram. It takes care of solving boundary condition
problem at the junctions. 132
5.12 Reservoir class diagram . 132
5.13 NodeStructure interface. It defines methods required to initialise and
handle flow situation at structures. The structure could be rectangular,
radial, circular gate or broad weir. . . . . . . . . . . . . . . . . . . . . . 13~~
xx
LIST OF FIGURES
5.16 CanalFlow data model. Shown are all tables, their attributes, data types,
relations and primary keys. Table names prefixed with 'm' are master
and 't' are transaction tables. Master tables are in which data trans-
actions (inserts and deletes) are limited while in transaction table they
happen frequently. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
xxi
LIST OF FIGURES
6.6 Schematic illustration of the looped system (C. 1) test. There are four
channels: Channel A, Channel B (B1 and B2 combined), Channel C
and Channel D. At the downstream end of Channel A the system di-
verges (splits) into Channels B and C, which then converge to Channel
D at their downstream sections. In CanalFlow test setup, Channel B is
divided into 2 sub-channels of equal length. . . . . . . . . . . . . . . . . 153
6.7 Schematic illustration of the weir (C.2) test case configuration. A weir is
placed between Channel A and Channel B, both of which have two cross
sections at a spacing of 100 m, a constant Manning's n of 0.014 and a
bed slope of 0.005 and 0.001, respectively. The fiow is from left to right. 155
6.8 Schematic illustration of the Irrigation network test (C.3) case. The
channels are numbered with prefix C, control structures with S and junc-
tions with J. Rectangular gates are used as control structures and they
are operated as per the given schedule. Constant depth of water is maln-
tained in the reservoir and water flows through main channels Cl to C4
and then through other channels connected to one of these channels. . . 158
6.9 Schematic illustration of the Contractions and expansions test (R.I).
The laboratory model consists of a 19.30 m long rectangular channel
with a removable sluice gate at 6.10 m and a constriction at 13.8 m.
The constriction is 0.1 m wide and 1.0 m long. Also shown are the
gauging stations SI to S4. . . . . . . . . . . . . . . . . . . . . . . . . . . 162
7.1 CanalFlow main screen. This is a basic GUI developed to help build
the network, run simulation and visualise the flow in a channel. It also
shows run-time simulation parameters such as current time, time step
(L\.t), Courant number, initial and final mass of the system and overall
progress. After the simulation is over, the screen on the right hand side
bottom corner displays the total time taken to execute the task. 166
xxii
LIST OF FIGURES
7.2 Front·end of the MySQL storage engine (refer Table 5.1, page 119). The
screen is divided in two parts; the left hand part shows list of available
tables and the right hand part shows columns of the selected table in
which data can be entered. Shown is the example of network parameter
table and corresponding data. The users aware of SQL can also use SQL
editor to insert/modify/delete the data. SQL scripts are generated for
building network and retrieving output data (which is then exported to
Excel sheet using MySQL Front) for each test case. . . . . . . . . . . . . lG7
7.3 The figure shows bed elevation and comparison of longitudinal water
surface profile generated by CanalFlow and as given by analytical solu·
tion plotted on the primary y-axis. It also shows the percentage error
between CanalFlow and analytical solution plotted on the secondary
y-axis for test (A) sub-critical flow (N.l.l) (B) super-critical flow (N.l.2)172
7.4 The figure shows bed elevation and comparison of longitudinal water
surface profile generated by CanalFlow and as given by analytical solu·
tion plotted on the primary y-axis. It also shows the percentage error
between CanalFlow and analytical solution plotted on the secondary
y-axis for test (A) super-critical to sub-critical flow (N.1.3) and (B)
sub-critical to super·critical to sub-critical flow (N.l.4) . . . . . . . . . . 173
7.5 (A) Bed elevation and comparison of longitudinal water surface profile
generated by CanalFlow and as given by analytical solution plotted on
the primary y-axis and percentage error between CanalFlow and ana-
lytical solution plotted on the secondary y-axis for super-critical to sub-
critical to super.critical flow test (N.l.5) (B) Longitudinal water surface
profile at the end of 6, 12, 18 and 24 hour. The flow boundary at the
upstream side (on the right) is kept constant while at the downstream
side (on the left) the depth boundary is kept constant for first 6 hours
and then increased linearly from 0.602 m to 2.0 m till 18:00 hour and
then kept at that level for another 6 hours. . . . . . . . . . . . . . . . . 175
xxiii
LIST OF FIGURES
7.6 Comparison of RMS errors for five sub-parts of the sub-critical, super-
critical and transitional flow (N.1) test case. It is evident from the
bar graph that amongst ISIS, HEC-RAS, MIKE 11, and CanalFlow,
CanalFlow has the lowest RMS error for all parts of the test. The av-
erage RMS error is approximately 4 times less than that of ISIS and
HEC-RAS and by approximately 6 times less than MIKE 11. . . . . . . 176
7.7 Comparison of normal depth and CanalFlow depth for (A) the sub-
critical and (B) the super-critical part of the test (C) shows the percent-
age error in depths generated by CanalFlow, ISIS, HEC-RAS, MIKE11
and the normal depth for the sub-critical flow, (D) shows the similar
percentage error for the super-critical flow test case. . . . . . . . . . . . 178
7.8 Ippen wave test: Water level and velocity time series at (A) and (B)
o km, (C) and (D) 25 km, (E) and (F) 50 km, and (G) and (H) 75
km. The Figure (C) and (E) shows that the calculated water levels
follow analytical solution closely between 2.0 and 2.5 hours and after
that it break from the analytical solution. At 75 km (G) calculated
solution follows the analytical solution between 3.0 and 3.5 hours and
after that variations are seen. The Figure (B,D,F and H) shows that the
velocities produced by CanalFlow follow the analytical solutions fairly
closely. In (D) and (F) it is seen that the velocities are identical to
analytical solution till 2.4 hours and after that they start breaking up.
While in (H) they are identical between 2.5 hours and 3.5 hours. .. . . 183
7.9 Ippen wave test: Water level and velocity profiles after (A) and (B) 2
hours, (C) and (D) 2.5 hours, (E) and (F) 3 hours. . . . . . . . . . . . . 185
7.10 Ippen wave test: Water level and velocity profiles after (A) and (B) 3.5
hours, (C) and (D) 4 hours. . . . . . . . . . . . . . . . . . . . . . . . . . 186
xxiv
LIST OF FIGURES
xxv
LIST OF FIGURES
xxvi
LIST OF FIGURES
7.23 Irrigation network (C.3) test case: Discharge hydrograph and mass bal-
ance at junctions at (A) J-1 where channel C-1 is split into C-2 and
0-5 with a rectangular gate on C-5 (B) J-2 where channel C-2 is split
into C-3 and C-20 with a rectangular gate on 0-20 and (0) J-3 where
channel C-3 is split into 0-4 and 0-33 with a rectangular gate on 0-33.
The percentage error at each junction is a error between discharge of the
upstream channel and sum of the discharges of split channels. . . . . . . 208
7.24 Dam-Break in a Rectangular Channel: The figure shows depth along the
length of the rectangular channel with unit width at 6 seconds. . . . . . 210
7.25 Dam break in Contractions and expansions test case (R.1): CanalFlow
predicted water depth time series using Minmod limiter at (A) 81, (B)
82, (C) 83 and (D) 84 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 212
7.26 Dam break in Contractions and expansions test case (R.1): CanalFlow
predicted depth time series using Monotonised Centered limiter at (A)
81, (B) 82, (C) 83 and (D) 84. . . . . . . . . . . . . . . . . . . . . . . . 213
7.27 Dam break in Contractions and expansions test case (R.1): CanalFlow
predicted water depth time series using 8uperbee limiter at (A) 81, (B)
82, (C) 83 and (D) 84. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 214
7.28 Dam break in Contractions and expansions test case (R.1): CanalFlow
predicted water depth time series using van Leer limiter at (A) 81, (B)
82, (C) 83 and (D) 84. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 215
7.29 Dam break in Contractions and expansions test case (R.1): CanalFlow
predicted water depth time series using no limiter at (A) 81, (B) 82, (C)
83 and (D) 84. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 216
xxvii
Nomenclature
Greek Symbols
Superscripts
n Time co-ordinate
For left going fluctuation from the right edge of the cell
+ For right going fluctuation from the left edge of the cell
Subscripts
i Space co-ordinate
xxviii
NOMENCLATURE
L Left state
n Time co-ordinate
R Right state
Units
GB Giga Byte
Ha Hectare
hr Hour
km Kilo Meter
min Minutes
m Meter
Other Symbols
A Jacobian matrix
xxix
NOMENCLATURE
Cd Coefficient of discharge
Fp" F;" Fp" Fp2 Pressure forces acting on a control volume (F)
F Vector of fluxes
L Left eigenvector
m1 Left side-slope
m2 Right side-slope
xxx
NOMENCLATURE
Q Discharge (L 3T- 1 )
R Right eigenvector
u Velocity (LT- 1 )
W Wave
I-DOne-Dimensional
3-D Three-Dimensional
AD Above Datum
xxxi
NOMENCLATURE
AI Artificial Intelligence
BC Boundary Condition
xxxii
NOMENCLATURE
CV Control Volume
DB Database
EA Environment Agency
FK Foreign Key
xxxiii
NOMENCLATURE
GC Grid of Characteristics
IC Initial Condition
xxxiv
NOMENCLATURE
JIT J ust-In-Time
MC Monotonised Centered
00 Object-Oriented
xxxv
NOMENCLATURE
PK Primary Key
QS Quasi-steady
xxxvi
NOMENCLATURE
WS Water Surface
XP eXtreme Programming
xxxvii
CHAPTER 1
Introduction
Today, in the most countries suffering water shortages, at the heart is the question of
whether a water crisis can be averted or whether water can be made more productive
by using it efficiently. Increasing the productivity of water is central to producing food,
fighting poverty, reducing competition for water and in ensuring adequate water for
nature. The more food produced with less water and/or with the same amount of water
will result in the less need for infrastructure development, the fewer conflicts among the
sectoral water uses, the greater local food security and the more water for agricultural,
household and industrial uses. However, to achieve such goals, major improvements are
still required in water resource use and irrigation technology management. The efficient
management of existing water resources at all levels therefore assumes a crucial role
in the intensification and sustainability of irrigation. One way of managing water
resources in a better way is to modernise the irrigation systems by implementing better
technologies. Typical project irrigation efficiencies1 are in the range of 25-35% (Burt
and Styles, 1999). Thus, there is a vast scope for increasing water use efficiency in the
irrigation sector. Raising efficiency of water use is a 'do or die' task and cost-effective
means for achieving this is through adoption of modern techniques in canal irrigation
management. Use of computer modelling in hydraulic simulation is one of the tools
IThe percentage of water delivered to the field that is used beneficially
1
1.1 Genesis
1.1 Genesis
Computer solution of the unsteady-flow equations began with single in-line canals and
solution with the method of characteristics. During the 1960s and '70s, significant
advances were made in simulating unsteady open-channel flow with computers (Clem-
mens, 1993). Finite Difference Methods (FDM) were developed in the 1970s (Liggett
and Cunge, 1975) and applied to branching and looped networks (Cunge et al., 1980).
Most of the applications of these techniques were for river- and sewer-system hydraulics.
Irrigation canal networks are somewhat different from other open-channel networks
because of the type of structures in use and the types of operations of these struc-
tures (A8CE, 1993). Few computer programs for canal- and river-system modelling
were developed by consulting firms and marketed as a part of their services to clients,
e.g., I8IS from HR Wallingford and MIKEll from DHI.
Unsteady flow of water in a one-dimensional (I-D) approximation is governed by
the shallow water or St Venant equations which are expressions for conservation of
mass and momentum (Garcia-Navarro et al., 1999). These are coupled first order,
hyperbolic, non-linear, partial differential equations. Most unsteady, irrigation canal
hydraulic simulation programs use FDM (see Table 2.1, page 45). While finite differ-
ence schemes such as the Preissmann (Cunge et al., 1980) and Abbott-Ionescu (Ab-
bott and Basco, 1989) schemes are widely used in hydraulic engineering, they are not
2
1.1 Genesis
suitable for certain types of problems, particularly flows which exhibit strong discon-
tinuities (Erduran et al., 2002). These limitations can be dealt with by introducing
shock capturing schemes (Glaister, 1995). Shock capturing capabilities are much easier
to implement using Finite volume methods (FVM) because they use the conservative
form of governing equations (Zhao et al., 1996). The FVM is based on discretisation
of the computational domain into smaller sub-domains over which the equations are
integrated. Finite volume methods are closely related to finite difference methods, and
a finite volume method can often be interpreted directly as a finite difference approxi-
mation to the differential equation. However, finite volume methods are derived on the
basis of the integral form of the conservation law, a starting point that turns out to
have many advantages (Leveque, 2002). FVM with shock capturing schemes applied
to the shallow water equations is thus suitable for modelling a variety of types of flow;
subcritical and supercritical; steady and unsteady; continuous and discontinuous and it
is thus emerging as a powerful tool in computational hydraulics (Erduran et al., 2002).
Majority of the numerical algorithms of the existing hydraulic models discussed
in the Section 2.6 were developed in the 1970s using sequential programming lan-
guages. The recent efforts on these hydraulic models were focused on the development
of a Graphical User Interface (GUI) with visualisation tools and Geographic Infor-
mation System (GIS) components (Kutija and Murray, 2002; Murray, 2003), without
any change in a programming paradigm of the core numerical algorithms. Thus, these
numerical algorithms have failed to take advantage from the advances in the software
programming. Further, Kutija (1998) has pointed out the lack of efforts in the develop-
ment of new numerical algorithms. In general, it was seen that the technology adoption
in the field of hydraulic simulation of irrigation canals has been quite slow and thus
the need was felt to bring in some industry-wide software development standards like
object-oriented programming, relational databases and configuration management to
build a highly-scalable software for hydraulic simulation of a canal irrigation network
using the latest numerical algorithm.
Originality of this research is the use of FVM in association with a high-resolution
shock- capturing technique to develop a database-driven, object-oriented hydraulic sim-
3
1.1 Genesis
ulation model for canal irrigation networks. This is the first known attempt where a
relational database structure has been developed and integrated into the core hydraulic
simulation model. The issues related to canal irrigation networks and some of the com-
mon flow scenarios occurring in a network are discussed in the Section 2.2 and 2.3.
It is envisaged that the newly developed model will be able to simulate various flow
scenarios in a canal network and could also be used to determine how best to operate
the canal to deliver water in an optimal manner to the users.
1.1.1 Rationale
Feldman (1992) discussed the initiative called Next Generation (NexGen) simulation
models and proposed the idea of providing the user with a better means to visualise
and understand the process being simulated, and to build more engineering expertise
into the models themselves. With recent advances in programming languages, data,.
bases, and hardware it is possible to develop models with tools to easily work with
input data, simulation processes, and output results. These models could also benefit
from web/intemet technology to run in a networked environment. Additionally, the
database-integrated system offers flexibility in storing and retrieving the network con-
figuration, operational and output data. Due to object-oriented design of the model,
it would be relatively easy to build an integrated modelling system and perform en-
hancements to the code with minimum effort. For example, a comprehensive irrigation
management system could be developed around the canal hydraulic model which takes
inputs from reservoir operation and/or canal scheduling modules such as Area and
Water Allocation Model (AWAM) (Smout and Gorantiwar, 2005). The 00 nature of
the developed canal hydraulic model code provides reusable components for the fu-
ture development and allows it to readily benefit from the new advances in numerical
schemes.
After the first web page written in 1990, the world has seen an exponential growth
in the use of intemet (http://www.zakon.org/robert/internet/timelinel).This
proliferation of intemet and the World Wide Web (WWW) calls for the development
of web-based hydraulic models which could be run using standard web browsers and
4
1.1 Genesis
a robust back-end to manage the data. Velickov et al. (1998) have described briefly
the use of Common Gateway Interface (CGI), ActiveX and Java technologies in remote
modelling, distributed computing and collaborative decision-making. Jolma (1998) has
reported development of the Rosa system for a reservoir operation using a Java-based
client. These kind of technologies will help end-users to set up their own models, execute
them and analyse the model's results from their offices and homes using the internet.
Thus, the philosophy of a development of the new model stems from the possible
migration of such a stand-alone model to a web-based one. It will be relatively easy
task to develop a web interface around the existing database-driven Java-based model
enabling multiple end-users to run the model simultaneously over the internet. Keeping
this in mind, it was decided to use Java! and MySQL2 as a preferred development
platform with a Concurrent Versions System3 (CVS) back-bone for the version control
of the code.
Java has attracted the scientific computing community because of its simple object-
oriented model and features like portability, safety, automatic garbage collection, and
multi-threading (Biount and Chatterjee, 1999). An effective database management
system has always been a high priority for modelling systems. MySQL is the world's
most popular open-source relational database management system (RDBMS). Its speed, .
scalability and reliability makes it suitable for small to medium sized databases (http:
/ /www.mysql.com). Discussion on object-oriented programming languages, relational
databases, CVS and a rationale behind selecting Java and MySQL is presented in the
Chapter 3.
5
1.2 Thesis Structure
capturing techniques. The model has been developed to specifically address the needs
of rotational irrigation systems.
Thus, the study aims at development of a 00 model for hydraulic simulation of
canal irrigation networks, referred to as CanalFlow, with the following objectives:
• Design and develop a FVM based model for hydraulic simulation of irrigation
canals using modern shock capturing techniques and object-oriented technology;
• Simulate flow conditions in the canal network for both steady and unsteady state,
based on a water delivery schedule;
• Design and develop a database strncture for the efficient storage and retrieval of
the necessary data and integrate the developed model with the database; and
• Test and validate the developed model for various flow scenaxios.
• Chapter 1 provides an introduction and sets a stage for the subsequent chapters.
6
1.2 Thesis Structure
• Chapter 6 describes the test strategy and description of the test cases used for
this research with detailed results and discussion in Chapter 7.
7
CHAPTER 2
2.1 Introduction
The objective of hydraulic simulation software is to create a flow situation in the model
which is dynamically similar to a flow in the prototype. The prototype behavior is then
inferred from appropriate quantities estimated by the software. The starting point of
an hydraulic model is a mathematical model. A wide variety of flows are described
by the mathematical models of a shallow-water type. These mathematical models
consist of partial differential equations (PDE) along with the initial and boundary
conditions. The mathematical description of the prototype is then transformed into
concrete numerical algorithms and procedures in the model and then solved numerically
by a computer.
This chapter discusses the necessity of a canal hydraulic modelling, various problems
faced in managing the canal networks, existing numerical methods, available hydraulic
simulation softwares, testing and validation of these softwares.
8
2.2 Why Canal Hydraulic Modelling Software?
The canal irrigation system delivers water from its source(s) to one or more points
of outlets on the downstream side. A typical irrigation system consists of a main
intake structure, a conveyance system, a distribution network, and a field application
system as shown in Figure 2.1. Operation of sluice gates in a canal network controls the
movement and behaviour of water in a canal system. The primary function of operation
is to manage the changes in flow and depth throughout the canal system. Several
methods for canal operation exist and they are described in detail in Clemmens and
Replogle (1989). The method of operation determines the water level to be maintained
in a canal pool in order to satisfy the user demands on the downstream side. A canal's
recovery characteristics - the speed and manner in which the canal recovers to a steady
state flow after a flow change - depends on the method of canal operation.
One of the root causes of the poor performance of irrigation systems is the lack
of scientiflc approach to their management. On most command areas served by a
canal, water is poorly distributed over area and time. A common shortcoming is that
tail-end users do not get sufficient water. If they are supplied, it may be insufficient
and unreliable. Conversely, head-end users often get too much water, either because
they have no choice or deliberately take water when they can and often more than
needed (Mandavia, 1999). This phenomenon is common in rotational irrigation systems
in South Asia and was reported by Mishra et al. (2001).
The operation of an irrigation network becomes complicated due to the following:
• Main canals range from tens to hundreds of kilometers and monitoring of dis-
charge and water level along the length of the canal is difficult. For example
Sardar Sarovar Project (SSP) in India has a main canal of length 478 km and a
distribution network of over 2750 km (http://www . sardarsarovardam. org);
• The command area spreads across a large area and is located far away from the
reservoir. In case of SSP, the distance of the farthest service area from the head
works is about 1000 km;
9
2.2 Why Canal Hydraulic Modelling Software?
Figure 2.1: Plan view of a typical irrigation network showing a reservoir, main canal,
control structures, distribution network and a field application system.
10
2.2 Why Canal Hydraulic Modelling Software?
• Changes are made in the releases at the head of the canal in anticipation of the
change of demand (closing ofrunning distributaries and opening of new ones, e.g.,
as for rotational water supply). The first aspect of this dynamic behaviour is the
varying time delay between an upstream operation and the corresponding change
in different locations of the system. Another aspect of this dynamic behaviour
is the fact that transient flow takes place most of the time, leaving only short
periods for steady flow. As an illustration, measurements made by International
Irrigation Management Institute (HMI) on Kirindi Oya Right Bank Main Canal,
Sri Lanka during 1991/92 on a 19 weeks period showed that steady flow periods
accounted for less than 20% of the time (Kosuth, 1996);
• Gates along the irrigation canal network are operated more or less simultaneously
and generally independently by various operators, thus inducing upstream and
downstream fluctuations in water levels and flows; and
• Pumps may be operated along the canals to irrigate the command area at a higher
elevation.
11
2.2 Why Canal Hydraulic Modelling Software?
• Minimise adverse hydraulic transients during flow changes. The transients hap-
pen due to opening and closing of the gates. It was found that the wave reflection
effects from fixed gates are significant and, in some situations, may cause embank-
ment erosion and bank overtopping (Misra et al., 1992). Gichuki et al. (1990)
developed a hydraulic model to simulate transient flow conditions in branching
canal networks;
12
2.2 Why Canal Hydraulic Modelling Software?
various structures manually. This leads to inequity in water delivery, and a mismatch
between supply and demand. Canal hydraulic models can help user to anticipate these
inequities in a canal network. Some of the other performance indicators related to
adequacy, equity, reliability and efficiency of the irrigation networks are explained in
Ha1crow (1998).
The main applications envisioned for unsteady-flow models are (ASCE, 1993; Hal-
crow, 1998) :
2. Gate scheduling: Unsteady flow models determine flow rates and water levels
from known boundary conditions and gate settings. Gate openings are calculated
to achieve desired water levels and flows in the system.
13
2.3 Problems Unique to Modelling Irrigation Networks
5. Flood routing: To analyse the impact of floods which may enter irrigation systems
and test the effectiveness of the available alternatives to route the flood waves
through the system in order to prevent or minimise the damage.
8. Training: To train design engineers and system operators on the basic principles
of unsteady flow in open channels and the consequences of changes made in sys-
tem design and operation on the flow and water levels in the system. The better
understanding of such issues by design and operation engineers should help them
make better designs and plan more effective and achievable operational proce-
dures.
Most of the unsteady simulation engines developed previously (listed in Table 2.1, page
45) were meant for river and sewer studies. Irrigation canal networks are somewhat
different from river or other open-channel networks for the following reasons:
• Canals generally branch out in the downstream direction while rivers converge;
• Type of structures in use and the type of operations of these structures. These
control structures are extremely important model elements and nearly all branch
points are controlled by structures; and
14
2.3 Problems Unique to Modelling Irrigation Networks
Thus, the computational environment for canal simulation is much more demanding due
to extreme variability and unsteadiness of the flow. Further, complications arise because
of the following flow scenarios which could be encountered (adapted from ASCE (1993);
Holly and Merkley (1993) to explain various flow scenarios):
1. Dry-bed flow: Filling and emptying of the canals occur regularly in many
irrigation systems. Dry bed flow involves both the advance of water in an empty
canal, and the subsequent dewatering of the canal after inflow is discontinued.
This flow situation is common in a rotational irrigation system where canals
are operated for a certain fixed duration and then closed until the next cycle of
irrigation (own experience). Thus, this situation is an important consideration in
unsteady hydraulic modelling.
15
2.3 Problems Unique to Modelling Irrigation Networks
such gate, and the conditions upstream and downstream, the flow regime may be
free or submerged. A discharge discontinuity usually results when the flow regime
through a gate changes from free to submerged, or vice-versa. This discharge dis-
continuity can introduce hydraulic instabilities in a real system. It can also cause
numerical instability in an unsteady-flow model. Another problem occurs when
a gated structure is opened to a level above the water surface. In this case the
structure behaves as a channel constriction (non-orifice flow condition) instead of
an orifice. Submerged weirs also behave as channel constrictions. The modelling
of these constrictions requires the application of a different algebraic equation to
relate flow depth with discharge at the structure, and the form of the equation is
such that convergence of the overall numerical solution is complicated. When the
regime changes back to that of orifice flow, there will be a sharp discontinuity in
simulated flow levels and discharge through the structure.
5. Inverted siphons and culverts: Inverted siphons and culverts usually appear
in a canal model as short computational reaches that obey specialised hydraulic
laws for the particular structure. These laws are the rather straightforward pipe-
flow equations with appropriate entrance, exit, and friction losses when siphons
or culverts are flowing full with submerged inlets and outlets. The flow may be
partially surcharged and partially free surface when the inlets/outlets are par-
tially submerged or free flowing. The computational problem arises from the
near impossibility of incorporating such detailed and dynamically changing hy-
draulics in the required two equations (continuity and momentum) between the
computational points at the ends of the computational reach. It is especially im-
portant to be able to reproduce the energy loss associated with a hydraulic jump;
and this cannot properly be done through ad hoc adjustment of loss or discharge
coefficients, especially in a dynamic environment. The only solution may be to
implement as complete a description of the conduit hydraulics as possible, as a
subcomponent of the simulation code.
16
2.3 Problems Unique to Modelling Irrigation Networks
17
2.4 Numerical Methods
ideal for training purposes. But in such operations it is essential to maintain the
computational efficiency and robustness of the underlying computational scheme.
The majority of numerical methods to solve the shallow water equations (SWE) fall
in one of the following categories: Method of characteristics (MoC), Finite difference
method (FDM), Finite element method (FEM) and Finite volume method (FVM).
Numerical solution of the equation provides answers at discrete points in the domain
called as grid points. The governing equations are approximated in order to obtain the
solution. Anderson, .Jr. (1995) defines discretisation as:
In short, it is the process of expressing general flow laws, written for a continuous
medium, in terms of discrete values at finite number of points in the flow field. Dis-
cretised flow laws are then numerically solved to furnish solutions, i. e., flow depth and
discharge.
The objective of this section is to review the numerical methods, in brief, before
discussing the hydraulic simulation softwares. The softwares discussed in the Section
2.6 use one of these numerical methods. It also sets a background for much of the
model theory discussed in Chapter 4.
18
2.4 Numerical Methods
In FDM, the partial derivatives in the governing equations are replaced by approximate
algebraic difference quotients, where the algebraic difference quotients are expressed
strictly in terrns of the flow-field vaxiables at two or more of the discrete grid points.
Thus the laws describing the evolution of the continuum are replaced by algebraic finite
19
2.4 Numerical Methods
ii. mass and a reduced Bernoulli equation with a one point boundary condition at
each end;
iii. mass and fully-reduced momentum equation with a one point boundary condition
at each end; and
20
2.4 Numerical Methods
iv. mass and partly-reduced momentum equation with a one point boundary condi-
tion at each end.
It was concluded that results were worst (as compared with other three methods)
when the momentum equation was used with a fully-reduced convective term, while
the method with a part-reduction of convective term (corresponding to the convective
term of the Bernoulli equation) yielded good results.
Several finite-difference schemes satisfying conservative and shock-capturing prop-
erties have appeared in literature, originally for Euler equations of gas dynamics (Roe,
1981). These schemes when extended to a higher order of accuracy for enhanced shock
resolution, employ sophisticated techniques to obtain smooth solutions and one class
of such schemes is referred to as Total Variation Diminishing (TVD) schemes (Fiedler
and Ramirez, 2000). Delis et al. (2000) showed that the modified flux linearised conser-
vative implicit (LCI) and monotonic upstream-centered scheme for conservation laws
(MUSCL) TVD were able to cope with all steady cases of sub-critical and mixed regime
flows with strong discontinuities.
MacCormack and TVD finite-difference schemes have been used by several authors:
Fiedler and Ramirez (2000) reported a MacCormack finite difference scheme to simulate
two-dimensional overland flow with spatially variable infiltration and micro-topography
using the hydrodynamic flow equations; Burguete and Garcia-Navarro (2001) pre-
sented high-resolution TVD schemes to solve shallow water flow with a source term;
Delis and Skeels (1998) used variety of finite difference TVD schemes, viz., second-
order symmetric, second-order upwind, predictor-corrector and the MUSCL scheme
to model flow in open channels; Liska and Wendroff (1999) presented finite-difference
schemes on uniform rectangular mesh and trapezoidal mesh for two-dimensional shal-
low water equations. The authors have applied Lax-Wendroff, Lax-Friedrichs and
weighted essentially non-oscillatory (WENO) composite schemes to the SWEs (Liska
and Wendroff, 1999). Macchione and Morelli (2003) compared the performance of
first-order and TVD second-order upwind flux difference splitting schemes, first-order
and second-order space-centered schemes with the TVD artificial viscosity term. The
21
2.4 Numerical Methods
schemes were applied to the dam-break cases in a dry frictionless horizontal chan-
nel, in a dry, rough and sloping channel and in a non-prismatic channel. Mohapatra
and Bhallamudi (1996) presented a MacCormack finite-difference scheme to simulate
a dam-break flow in channel transitions. Tseng (1999) presented Wgh-resolution non-
oscillatory shock-capturing finite-difference schemes to solve steady and unsteady one-
dimensional flow with steep waves in channel. Further, Tseng and Chu (2000b) reported
a finite-difference predictor-corrector TVD scheme for the computation of unsteady one-
dimensional dam-break flows. Numerical modelling of open channel flows using explicit
finite difference schemes is constrained by the Courant, Ftiedrichs, and Lewy (CFL)
stability criteria which limits the time step. To overcome this limitation Yost and Rao
(2000) introduced multiple grid algorithm coupled with a second-order accurate Mac-
Cormack scheme. In a multiple grid algorithm, time marching is accompanied by a
spatial marching. At any time level, the solution is marched over progressively from a
fine grid to a coarse grid level. TWs was the first attempt to use this algorithm to solve
a flow problem with moving shock fronts.
In the FEM, solutions are approximated by either eliminating the differential equation
completely (steady state problems), or rendering the PDE into an equivalent ordinary
differential equation, which is then solved using standard techniques such as finite dif-
ferences. The FEM relies upon a variational formulation (Zienkiewich and Taylor, 2000)
of the governing equations. Its main advantage stems from its rigorous mathematical
foundation that allows error to be estimated from the results (explained in detail in
Zienkiewich and Taylor (2000)). However, it is conceptually more difficult and far fewer
authors choose this option (Alcrudo, 2002). The program TELEMAC (Hervouet and
Petitjean, 1999; Hervouet et al., 2000) was one of the first to be applied to dam-break
modelling using the finite element method. Hanke (1998) has successfully used a semi-
discrete finite element method for the computation of transient shallow water flows
and applied the method to dam break flood propagation in experimental flumes and to
solve irrigation problems (Hauke, 2002). A discontinuous Galerkin FEM is described
22
2.4 Numerical Methods
by Aizinger and Dawson (2002) for the supercritical, river inflows and tidal flows as
well as for contaminant transport. Aronica et al. (1998) presented a hyperbolic model
for the simulation of flood wave propagation on initially dry land and Heniche et al.
(2000) developed a model to simulate two-dimensional free surface flow in rivers and
estuaries with drying and wetting using the Eulerian method to predict the position of
the moving boundary. Quecedo and Pastor (2002) showed how classical finite element
formulation can be used to solve the wetting drying problem in an efficient manner and
Ribeiro et al. (2001) implemented the generalised residual minimum (GMRES) itera-
tive solver using an edge-based data structure. When used to perform matrix-vector
product operations, the edge-by-edge formulation results in a remarkable improvement
in the performance over a conventional element-by-element technique, in terms of the
memory and CPU requirements. Sheu and Fang (2001) developed a high-resolution
Taylor-Galerkin finite element model to simulate SWEs for bore wave propagation in
2-D. Shock-capturing has been achieved using flux-corrected transport filtering tech-
nique which is similar to the one discussed in the Section 2.4.5.
Finite-volume schemes solve the integral form of the shallow water equations in com-
putational cells. Instead of point-wise approximations at grid points, the FVM divides
the domain into grid cells and approximates the cell average of a function, which is
the total integral divided by the volume of the cell (Leveque, 2002). These values are
modified in each time step by the flux through the edges of the grid cells. It is essential
to determine good numerical flux functions that approximate the correct fluxes reason-
ably well, based on the approximate cell averages. One of the methods of evaluating
fluxes at cell faces is by solving a Riemann problem, which is discussed in the Section
4.3.
FVMs have several advantages over FDM and FEM approaches. FVMs combine
the simplicity of FDMs with the geometric flexibility of FEMs (Mingham and Causon,
1998). Thus, the methods can be applied using an unstructured grid system as FEM;
generally FVM needs less computational effort than FEM. FVM is based on the integral
23
2.4 Numerical Methods
form of the conservation equations, thus a scheme in conservative form can easily be
constructed to capture shock waves (shock-capturing property). By discretisation of the
integral form of the conservation equations, the mass and momentum remain conserved,
exactly. As discussed earlier, the key problem in FVM is to estimate the normal
flux through each cell interface and methods which estimate this flux should handle
the intrinsic directional property of signal propagation in the SWEs. For example,
in I-D unsteady, sub-critical flow, information propagates from both upstream and
downstream, whilst, in super-critical flow, information propagates from upstream side
only. Methods capable of doing so are those which employ Riemann solvers (which
solve the Riemann problem and discussed in the Section 4.3). When used in FVM, the
Riemann solver is especially suitable for capturing flow discontinuities which arise in
a dam break, bore wave propagation in a channel, or rapid opening and closing of a
sluice gate (Mingham and Causon, 2000). As already discussed in the Section 2.3 of
this chapter, these flow situations are common in canal irrigation networks.
In last ten years, FVM has been widely used to solve problems in the field of
hydraulic engineering. Few examples in the literature which refer to solving I-D, 2-D
problems using FVM are:
.• Dam~break (Brufau and Garcia-Navarro, 2000; Caleffi et al., 2003; Lin et al., 2003,
2005; Nujic, 1995; Tseng et al., 2001; Valiani et al., 2002);
• Hydraulic jump and bore waves (Glaister, 1995; Hu et al., 1998; Yoon and Kang,
2004; Zhou and Stansby, 1999);
• Actual river flow problem (Caleffi et al., 2003; Capart et al., 2003; Goutal and
Maurel, 2002; SchuIz and Steinebach, 2002); and
24
2.4 Numerical Methods
• Simulation of flow over wet/dry beds (Yoon and Kang, 2004; Zoppou and Roberts,
1999) and in prismatic/non-prismatic channels (Sanders, 2001).
As seen from the list, the FVM has been applied to solve a variety of hydraulic
engineering problems but no known reference in the literature exist where a FVM
based high-resolution model has been developed for the hydraulic simulation of a canal
irrigation network. However, Zoppou and Roberts (1999) and Erduran et al. (2002)
have applied their model for a single irrigation channel to test the ability of the model
to simulate a wetting and drying phenomena.
The shallow water mathematical model is a set of hyperbolic PDEs that typically
represents propagation phenomena. In order to correctly reproduce such physical
phenomena, where flow variables may be discontinuous, upwind schemes were devel-
oped (Hirsch, 2002). A family of methods that belong to this class is the Godunov-type
ones. In 1959, Godunov gave an exact solution of the Riemann problem. The Riemann
problem is defined as a standard shock-wave problem describing the interaction of two
different constant-property states (for example depth and velocity) and the methods in
which a numerical solution is built up from a series of elementary Riemann problems
at cell interfaces are called Riemann solvers. Riemann solvers are discussed in detail
in Hirsch (2002), Leveque (2002), and Toro (2001). Riemann solvers can be classi-
fied based on their solution as exact and approximate solver. Exact Riemann solvers
are more accurate but computationally more expensive. Roe's approximate Riemann
solver (Roe, 1981) has led the way to many others (Brufau and Garcia-Navarro, 2000;
Delis et al., 2000; Hu et aI., 1998; Nujic, 1995; Sanders, 2002; Valiani et al., 2002) while
the exact Riemann solver is used by Alcrudo and Benkhaldoun (2001), Karelsky et al.
(2000), and Zanuttigh and Lamberti (2002a). The exact Riemann solver used in this
research is described in the Section 4.4.2, page 91. The Riemann solvers like Osher,
HLL, HLLC exist and are discussed elsewhere (Erduran et al., 2002; Toro, 2001).
25
2.4 Numerical Methods
26
2.5 Attributes of Numerical Solution Techniques
and Chu, 2000a,b). Thus, the concept of limiters is introduced. High-resolution limiters
used in this research are discussed in the Section 4.5.1, page n6.
In this section, various attributes of the numerical solution techniques, viz., applicabil-
ity, accuracy, convenience and robustness are described:
• Applicability refers to the range of system geometries and flow regimes that the
technique is designed to handle. Examples of the system geometry are control
gates, cross-section, pumps, channel branches, and looped networks. Examples of
flow regimes are backwater effects, entry of water into a dry channel, dewatering a
sloping channel with upstream inflow cut off, flow overtopping gates, a hydraulic
jump downstream of a free flowing gate, supercritical flow, and bores.
• Accuracy in a model means that the values of depth and discharge it predicts
are sufficiently close to the physical reality to be useful. The model, if allowed to
run indefinitely with constant and equal inflow and outflow must converge to a
steady-state solution, without adding or losing water. Ultimately, the best test
of accuracy for a theoretically sound model run within its range of applicability
is agreement with field and/or laboratory tests.
27
2.6 Hydraulic Simulation Softwares
Thus, the numerical solution technique used to develop the model should be applicable
to a variety of flow conditions, should be accurate in estimating the flow parameters
and convenient to use. Finally, the technique should be easy to implement using some
programming language, the developed code should execute quickly, should have small
memory requirements, use modular construction, and provide convenient method to
input the data and present the output which can be interpreted in relatively short
time (Lai, 1986).
In 1987, the Irrigation and Drainage Division of ASCE presented a symposium on Plan-
ning, Operation and Rehabilitation of Irrigation Water Delivery Systems (Zimbelman,
1987). One topic of focus was irrigation water distribution system management, and
basic technical details on water-control hardware and strategies were provided. As a
result of this symposium, an ASCE task committee was formed in 1989 to evaluate and
examine unsteady flow simulation models. The committee was interested primarily in
the models' ability to simulate water level and flow variations in canal systems with a
variety of gates and weirs. A special issue of the Journal of Irrigation and Drainage
Engineering (119(4), July/August, 1993} presented the committee's findings, including
a description and evaluation of the various unsteady flow programs. This special issue
carried a review of five hydraulic models in use at that time. Since then new versions
of some of these models were released and many more new softwares were available
for the hydraulic simulation in an open-channel network. Thus, it was thought to be
appropriate to review the existing hydraulic simulation softwares, look at their features
and known limitations before embarking on the development of a new model. In short,
28
2.6 Hydraulic Simulation Softwares
this section is a snapshot of hydraulic models developed and used in the last three
decades.
In this section, various hydraulic simulation softwares are discussed with a brief de-
scription of the software and underlying numerical solution technique provided. These
models are classified based on their numerical method, system platform and availability.
2.6.1 BRANCH
29
2.6 Hydraulic Simulation Softwares
Canal CAD (Holly and Parrish Ill, n.d.) was the first unsteady flow simulation soft-
ware developed primarily to test automatic canal control algorithms (Clemmens et al.,
2005). It is essentially a menu-driven front end to the unsteady simulation model
CARlMA (Holly and Parrish III, 1993), which is discussed in the Section 2.6.4. Canal-
CAD simulates both steady and unsteady flow in canal systems with manual or au-
tomatic gates. Hydrodynamic channel computations are based on the full, dynamic,
one-dimensional St Venant equations. The hydrodynamic simulator is a mature appli-
cation designed originally to model extensive riverine systems. CanalCAD comprises
the synthesis of a mature, reliable dynamic-equation solver; a menu-driven interface
for canal definition and results processing; and user-customized access to gate-control
algorithms. CanalCAD is designed for the design, analysis, and operation of irrigation
canals comprising sub critical flow in a single in-line system of pools and appurtenant
structures including turnouts, in-line weirs, check structures, culverts and off-line stor-
age reservoirs. The system supports a high degree of user guidance in canal description
and simulation post-mortem diagnosis.
Canal CAD utilises simulation technology developed based on the algorithms de-
veloped by Cunge, Preissmann, Chevereau, Holly, and others at Societe Grenobloise
d'Etudes et d'Applications Hydrauliques (SOGREAH), France. Computations are
based on the full, dynamic, one-dimensional St Venant equations and utilise the Preiss-
man four-point implicit finite difference scheme.
CanalCAD allows seamless inclusion of a user-supplied control algorithm for the
control of gates by means of compiling and linking the user's FORTRAN subroutine
into the simulator's binary code. This is done entirely from within CanalCAD by means
of a few simple commands, so that the user need not attend to compiler set-up and
management issues (other than installation).
30
2.6 Hydraulic Simulation Softwares
2.6.3 CanaIMan
CanaIMan (Canal Management Software) was developed for performing hydraulic sim-
ulations of unsteady flow in branching canal networks (Merkley, 1997). The model was
designed for use in design, analysis, operational and training activities, thus the name
CanalMan. The model can be used to simulate canal operations in a manual mode, and
it can generate proposed operating schedules through a centralised automatic mode.
Several common local gate automation schemes are also included in the model, and
these can be easily selected and calibrated through the model interface.
CanalMan implicitly solves an integrated form of the St Venant equations of conti-
nuity and momentum for one-dimensional unsteady open-channel flow. Computational
nodes are used internally by the model, and they are automatically inserted along the
length of a canal reach, between the system layout nodes.
The model is highly interactive and includes integrated data editing capabilities,
with numerous options for canal system configuration, hydraulic simulations, and out-
put of results. Canal networks are bnilt interactively by inserting and arranging nodes
graphically in a system layout window on-screen, where nodes represent locations of
flow control structures and channel bifurcations. Simulations can be started by filling
an empty canal system, continuing a previous simulation, or from a specified steady
or unsteady flow condition. The model will directly simulate the layout of most canal
systems, including branching canals.
The model cannot handle channel de-watering, rapid flow changes (such as bores
and surges), negative flow at the structures, hydraulic jump and a supercritical flow.
Looping canal systems are also not handled. The computational time step can be
from one to ten minutes, in whole minutes. Channel cross-sections can be either circu-
lar, trapezoidal or specified coordinates. Trapezoidal sections can be non-symmetrical
(that is, different side slopes on either side), and by specifying up to 16 cross-section
coordinates one can work with cross-section of an arbitrary shape. However, the cross-
sectional shape and size can change from reach to reach, but not within a reach. Thus,
the model does not handle non-prismatic or natural channel sections. However, in
31
2.6 Hydraulic Simulation Softwares
many cases natural channels can be adequately approximated (for hydraulic modelling
purposes) by a trapezoidal section, and "section change" structure types can be used
to model some kinds of non-prismatic channels (Merkley, 1997).
2.6.4 CARIMA
The Calcul des Rivieres MailJees (CARlMA) code simulates unsteady free-surface flow
in simple or multiply-connected systems ofrivers or canals (Holly and Parrish IH, 1993).
CARlMA was originally developed as a batch code, and has most often been used in
this form. However, an unsteady flow simulator based on CARlMA's one-dimensional
algorithms, restricted to a single suite of prismatic canal pools, has been developed
and embedded in the CanalCAD (Section 2.6.2) system for the design and analysis of
canals with automatic control.
The simulation uses the Preissmann implicit finite difference method for the so-
lution of the complete St Venant equations and appropriate equations for hydraulic
structures. The topological generality of CARlMA (branched and looped network)
naturally accommodates features like downstream influence across structures with no
modification to governing hydraulic equations. The nonlinear algebraic equations re-
sulting from application of the Preissmann method are solved in a Newton-Raphson
context, although a single iteration is generally sufficient. Special features of the code
include automatic topology recognition, a user interface for fully implicit incorporation
of specialised supervisory and global control schemes, without limitation on network
complexity or size of data tables.
2.6.5 DUFLOW
DUFLOW (from "Dutch flow") is a software package for simulating one-dimensional un-
steady flow in open-channel systems (Clemmens et al., 1993). The main computational
procedures in DUFLOW are based on IMPLIC, a main-frame-based FORTRAN pro-
gram. The program is designed for simple networks of channels with simple structures.
The program is intended for use in all types of open channels, e.g., rivers, navigation
32
2.6 Hydraulic Simulation Softwares
channels, drainage canals, and not just irrigation canals. It has a user interface for con-
veniently entering a description of the network conditions. Although the program can-
not handle some of the more-sophisticated modelling needs, viz., dry beds, automatic
gates, it is useful for the first-time users of canal network software. DUFLOW is now
part of the DuHow Modelling Studio (DMS), which is an integrated water management
modelling system (http://www.mx-groep.nl/producten/Duflow /Duflow-web/).
A four-point implicit Preissmann scheme is used to solve the complete St Venant
equations of continuity and momentum. The user has an option of selecting the so-
lution either by linearised or fully non-linear versions of the equations, the non-linear
being solved with a Newton-Raphson type scheme that starts from the linearised re-
sults. External and internai boundary conditions are expressed and solved within the
Preissmann scheme. The linear set of equations is solved by Gaussian elimination.
DUFLOW allows the user three options related to the flow inertia. The user can
choose to retain the Froude term (inertia fully accounted for), omit the Froude term
(zero-inertia), or restrict the Froude-term effect to not exceed frictional-resistance ef-
fects (damped). This term has a small impact for large, slow-moving irrigation canals,
but it can be signiflcant for higher velocity canal flows. Removing the inertial term
allows waves to travel only in the main flow direction, and it removes reflected waves
from downstream structures as the flow is dominated by gravitational forces.
2.6.6 FEQ
The Full Equations (FEQ) model for the simulation of one-dimensional unsteady flow
in open channels and through control structures was first developed in 1976 to model
the Sanitary and Ship Canal, Chicago. FEQ simulates flow in a open-channel network
by solving the full, St Venant equations for one-dimensional unsteady flow in open
channels (Franz and Melching, 1997). The structure of the program is designed to
follow the structure of a open-channel system while providing maximum generality and
flexibility of description. A channel system that is simulated with FEQ is subdivided
into three broad classes of flow paths: (1) channel reaches (branches), (2) parts of the
channel system for which complete information on flow and depth are not required
33
2.6 Hydraulic Simulation Softwares
(dummy branches), and (3) level-pool reservoirs. These components are connected by
special features or hydraulic control structures, such as junctions, bridges, culverts,
dams, waterfalls, spiJIways, weirs, side weirs and pumps.
The principles of conservation of mass and momentum are used to calculate the flow
and depth throughout the channel network resulting from known initial and boundary
conditions with an Preissmann four-point implicit finite-difference scheme.
2.6.7 FESWMS-2DH
34
2.6 Hydraulic Simulation Softwares
Galerkin's finite element method weights the governing equations over the entire
solution domain. The weighting process uses integration, which is carried out numeri-
cally using Gaussian quadrature on a single element. Repetition of the integration for
all elements that comprise a solution region produces a system of nonlinear algebraic
equations when the time derivatives are discretised. As the system of equations is non-
linear, an iterative solution procedure is needed and is computationally expensive since
it is 2-D.
2.6.8 FourPt
2.6.9 HEC-RAS
35
2.6 Hydraulic Simulation Softwares
36
2.6 Hydraulic Simulation Softwares
2.6.10 Hivel-2D
2.6.11 ISIS
ISIS consists of range of software tools for modelling rivers and their catchments. Its
complete open channel modelling solution models open channel systems containing
loops, branches, floodplain conveyance and storage. It incorporates standard equa,-
tions and modelling techniques for structures including weirs, sluices, bridges, culverts,
pumps, siphons, orifices and outfalls. Logical control of moving structures systems
can be simulated either in full unsteady mode, or as an advanced steady state simula-
tion (http://www.wallingfordsoftware.com/produets/isis!).ISIS systems (qual-
37
2.6 Hydraulic Simulation Softwares
ity, sediment and probability distribution module) are jointly owned, developed and
supported by Wallingford Software Ltd and Halcrow Group Ltd.
ISIS is a one-dimensional, full hydrodynamic simulator for modelling flows and lev-
els in open channels and estuaries. It uses a classical four-point implicit Preissmann
scheme to solve the St Venant equations. ISIS can model complex looped and branched
networks, and is designed to provide a comprehensive range of methods for simulating
flood plain flows. It incorporates both unsteady and steady flow solvers, with options
that include simple backwaters, flow routing and full unsteady simulation. The simu-
lation engine provides a direct steady-state solver and adaptive time-stepping methods
to optimise run-time and enhance model stability.
2.6.12 Mike 11
2.6.13 MODIS
MODIS (from Modelling Drainage and Irrigation Systems), developed at Delft Uni-
versity, was developed to investigate the hydraulic performance of dynamic controlled
38
2.6 Hydraulic Simulation Softwares
2.6.14 NOAH
39
2.6 Hydraulic Simulation Softwares
to support the Preissmann scheme on a non-staggered grid and the NewC scheme (K u-
tija and Hewett, 2002) on a staggered grid are also available (Kutija and Murray, 2002).
The GEA implements concepts derived from graph theory, and gives a general, uncon-
ditionally stable, accurate and fast method of solution for solving a flow problem in
a network. This general network algorithm provides a fully-automated instantiation
and control procedure for combined dendritic and looped networks providing solution
speeds independent of network complexity. This novel approach consists of two main
parts: the definition of the control mechanism and the actual calculation of the flow
conditions. The control mechanism required to steer the calculation procedure is ob-
tained on the basis of the graph-theoretical approach. The calculation of flow conditions
consists of two parts, the one for flow conditions in the dendritic part of the graph and
the other for flow conditions in the cyclic part of the graph.
2.6.15 SIC
The SIC software (Simulation of Irrigation Canals) is one of the hydraulic models
developed by Centre National du Machinisme Agricole du Genie Rural des Eaux et des
Forets (CEMAGREF). The first developments on hydraulic numerical modeling started
at CEMAGREF in the early 1970s (Malaterre and Baume, 1997). The very first version
of the SIC model was developed for the International Irrigation Management Institute
(HMI) for a canal located in the south coast of Sri Lanka (Kirindi Oya Right Bank
Main Canal). One purpose was to create a model easily usable by canal managers as
a decision support tool, in order to help them in the daily operation and maintenance
of their network. Since this first application was promising, CEMAGREF, with other
partners, decided to develop a new standard version of this software, that could be used
on most of the irrigation canals world-wide.
The SIC hydraulic model solves the complete St Venant equations using the classical
implicit Preissmann scheme. It consists of a topographical module used to enter the
topological and geometrical data along with a steady and an unsteady flow module.
The unsteady flow calculations can be performed from a steady state results or previous
40
2.6 Hydraulic Simulation Softwares
unsteady flow state. SIC does not handle advances on a dry bed, channel dewatering
and hydraulic jumps but automatic control of canal structures is possible.
2.6.16 Sobek
In 1993, Delft Hydraulics introduced a new unsteady flow simulation software package,
Sobek - an integrated software package which links river and sewer systems. Sobek
can handle a canal/sewer network of any size and has a robust numerical core (http:
//""'" sobek.nl!).
Sobek uses the Delft Hydraulics Laboratory implicit, staggered grid, finite differ-
ence scheme to solve the 1-D St Venant equations (Cunge et al., 1980). The scheme is
based on the concept of nodes where the water surface elevations are computed. These
computational nodes are connected to adjacent nodes on the left and right through
discharge equations. Thus, a computational grid is set up that consists of water sur-
face elevations (h-points) and flow rates (Q-points). The h-points are located at the
computational nodes and the Q-points are the connections between the nodes (Cunge
et al., 1980).
2.6.17 TELEMAC
41
2.6 Hydraulic Simulation Softwares
software is supplied in modular form allowing users to create the package most suited
to their needs.
TELEMAC-2D solves, through a finite element method over non-structured grids
consisting of triangular elements, the St Venant shallow water equations in two horizon-
tal dimensions, primary solution being two horizontal components of depth-averaged
velocities and water depth. Since its early inception, TELEMAC-2D has been applied
and validated extensively. The research community, as well as practicing engineers,
have used TELEMAC in variety of applications, for example, floods (Asaro and Paris,
2000), dam breaks (Hervouet, 2000), tides (Malcherek, 2000) and meandering chan-
nels (Rameshwaran and Shiono, 2002). TELEMAC-3D solves the three-dimensional
Reynolds Averaged Navier-Stokes equation (RANS) with free surface and both hy-
drostatic and non-hydrostatic modes. Mostly suited for the flows which are vertically
heterogeneous or purely three-dimensional, TELEMAC-3D allows the study of 3-D flow
structures generated in compound meandering channels.
2.6.18 UNET
UNET simulates one-dimensional unsteady flow through a full network of open chan-
nels. One basic element of a full network problem is the split of flow into two or more
channels. For sub-critical flow, the division of flow depends on the stages in each of the
receiving channels. These stages are a function of channel geometry and downstream
backwater effects (USACE, 2001).
UNET uses the four-point implicit finite difference scheme for solving I-D unsteady
flow equations. The resulting non-linear algebraic equations are solved using a Newton-
Raphson iteration technique. For a reach of river there are N computational nodes
which bound N - 1 finite difference cells. From these cells, 2N - 2 finite difference
equations are developed. Because there are 2N unknowns (AQ and Az for each node),
two additional equations are needed. These equations are provided by the boundary
conditions for each reach, which for sub-critical flow, are required at both the upstream
and downstream ends. For supercritical flow, boundary conditions are only required at
the upstream end. However, UNET is limited to sub-critical flow conditions.
42
2.6 Hydraulic Simulation Softwares
2.6.19 USM
USM (from "Unsteady Model") developed by the United States Bureau of Reclamation
(USBR) is an hydraulic simulation computer program that models gradually-varied un-
steady flow in canal systems. The primary purpose and application of USM has been
hydraulic analysis during the design of new canals and canal-control systems (Burt,
1990; Rogers and Merkley, 1993). USM analyses canal systems with a variety of struc-
tural boundary conditions, including gates, weirs, siphons and pumps. The gates can
have either manual or automatic gate contro!' Branching canals, supercritical flow, and
channel dewatering cannot be analysed.
USM solves the St Venant equations using the MoC, thus, providing high degree
of computational accuracy. Two variations of the MoC solution are available in the
program: complete grid of characteristics and specified time interval. In the complete
grid of characteristics solution, the calculation time step varies as water depth and
wave speed change. The specified time interval solution uses a fixed time step, requir-
ing interpolation to yield depth and flow values at regular time intervals. Numerical
accuracy is particularly high when the grid of characteristics option is selected.
Although, the MoC yields an accurate numerical solution, it requires a large number
of computations largely because the calculations must also satisfy the Courant condi-
tion, typically resulting in a maximum time step of about 60 seconds. Therefore, the
solution using the method of characteristics can be slower than other solution meth-
ods. USM is more efficient for solving problems of short duration involving rapid flow
changes than for problems of long duration with gradual changes.
It is apparent from the Table 2.1 that majority of the models (15 out of 20) fall un-
der the broad category of implicit finite difference scheme, which includes commercially
available software ISIS, Mike 11 and the widely used REO-RAS. NOAR and CanalFlow
are the only softwares implementing object-oriented numerics while only CanaIFlow im-
plements high-resolution shock capturing methods. The CanalFlow model is presented
43
2.7 Testing and Validation of the Software
in the Chapter 4 and 5. As a direct result of this research, a new category of canal
hydraulic simulation software has been created (Category 6 in the Table 2.1).
44
2.1 Testing and Validation of the Software
Table 2.1: Summary of hydraulic models (Notation: SC - Source code; Exe - Executa-
bles; PR - Proprietary; PU - Public). Included as a new category (Category 6) is a
FVM based high-resolution shock capturing methods.
Model Type Platform Affiliation I Domain
I se I Exe
Category 11 Impllclt Finite Difference (Prel8smann four-point scheme)
BRANCH 1-0, Unsteady Flow and UNIX or DOS USGS,USA PU PU
transport FORTRAN
CanalCAD 1-0, Unsteady DOS IIHR, University of Iowa, PR PR
FORTRAN USA
CanalMan 1-0, Unsteady MS Windows Utah State University, USA PR PU
CARIMA 1-0, Unsteady DOS SOGREAH, France PR PR
DUFLOW 1-0, Unsteady DOS Delft University, PR PR
IHE, Rijkswaterstat, Th.
Netherlands
FEQ 1-0, Unsteady DOS or UNIX USGS,USA PU PU
FORTRAN
FourPt 1-0, Unsteady DOS or UNIX USGS,USA PU PU
FORTRAN
HEC-RAS 1-0, Unl:lteady MS Windows HEe, USACE, USA PR PU
Hydraulic Multi-User
Modelling System
ISIS 1-0, Unsteady, Hy· MS Windows Walllngford Software and PR PR
draulic/Hydrodynamic Halcrow, UK
Modelling System
MODIS 1-0, Unsteady"Hydrody_ DOS Delft University, PR PR
namic Model The Netherlands
SIC 1-0, Unsteady MS Windows CEMAGREF, France PR PR
UNET I-D Unsteady MS Windows HEC, USACE, USA PU PU
Category 21 Implicit Finite Difference (Abbott-Ionescu scheme)
Mike 11 1-0 Steady, MS Windows DHI, Denmark PR PR
2-D Unsteady
Hydraulic/Hydrodynamic
Modelling System
NOAHID 1,2-0 Unsteady MS Windows Newcastle University, PR PR
Hydraulic model Delphi UK
Category SI Implicit Finite Difference (Delft Hydraulics Laboratory SCheme)
Sobek 1-0, Unsteady MS WlndowlJ WL Delft Hydraulics, PR PR
Hydraulic modelling The Netherlands
system
Category 4= Method of Characteristics
USM I 1-0, Unsteady I DOS I USBR, USA I PR I PR
Category SI Galerkln Finite Element Method
FESWMS 2-D, Unsteady UNIX USGS,USA PU PU
Hydrodynamic FORTRAN
Hivel-2D 2-D, Unsteady MS WindowlJ CHL, USWES, USA PR PR
TELEMAC 2-D, 3-D Unsteady MS WindowlJ EDF-LNH, France PR PR
Hydrodynamic UNIX
Category 61 FVM based High-Resolution Shock Capturing Methods
CanalFlow 1-0, Unsteady Independent Loughborough University, Not Not
Hydraulic Java and MySQL UK decided decided
Multi-User
45
2.1 Testing and Validation of the Software
Borthwick et al. (2001), and Sheu and Fang (2001). In few cases only the numerical
behaviour has been tested, and therefore only truncation and round-off errors have
been evaluated. In other cases the validation has been made against experimental data
obtained in laboratory models of a flood event (Aureli et al., 2000; Jha et al., 2001;
Tseng and Chu, 2000a).
Validation work against real world data have been performed in Hervouet and Pe-
titjean (1999), and Hervouet (2000) regarding the Malpasset dam break. Few authors
have tested the canal hydraulic models using real canal data (Clemmens et al., 1993;
Holly and Parrish IH, 1993; Malaterre and Baume, 1997; Merkley and Rogers, 1993;
Rogers and Merkley, 1993; Schuurmans, 1993). Primarily, data from the Central Ari-
zona Project (CAP) (Rogers et al., 1993) and CalPoly Model (Parrish III and Burt,
1993) was used for the testing. Rogers and Merkley (1993) used data from CAP to
test the USM model. Direct comparison of water levels calculated by USM to field
measurements was carried out. USM was found to be fairly accurate in estimation, pri-
marily because of its use of MoC scheme. But USM was more suitable for simulating
fiow conditions with a long-term gradual changes in a canal network. Holly and Parrish
III (1993) applied CARIMA model to the CalPoly model canal. Use of CARIMA was
tested with the control algorithm CARDD for canal automation purposes. Merkleyand
Rogers (1993) tested CANAL model using a hypothetical data for which known theo-
retical solutions were available and also using field data from CAP. CANAL was found
to be a good tool for modelling daily canal operations. Clemmens et al. (1993) tested
DUFLOW using test specified by Contractor and Schuurmans (1993) and CalPoly data.
DUFLOW was primarily developed for simulating flow in drainage networks hence it
did not handle the hydraulic structures as conveniently as other models, viz., CANAL,
CARIMA, and USM. Malaterre and Baume (1997) tested SIC model for the Kirindi
Oya Right Bank Main Canal in Sri Lanka and for the test caseS specified by Contractor
and Schuurmans (1993). Kumar et al. (2002) applied CanalMan (Merkley, 1997) to
the Right Bank Main Canal of the Kangsabati irrigation project in India. CanalMan
was calibrated for the first irrigation period and then validated for the fourth irriga-
tion period in Kharif season. The observed and simulated values of inflow at system
46
2.1 Testing and Validation of the Software
source were found to be in a close agreement for both calibration and validation pe-
riod. Further, MacDonald et al. (1997) developed analytic benchmark solutions for the
open-channel flows.
However, there are few shortcomings with the above-mentioned test cases for irri-
gation canals:
• Do not cover broad range of modelling scenarios, viz., flow transitions and looped
systems;
Thus, there was a serious need for exhaustive benchmark test cases in order to test
newly developed hydraulic simulation software packages. This gap was recently filled
by new benchmark tests developed under the Department for Environment Food and
Rural Malrs (DEFRA)/Environment Agency (EA) Flood and Coastal Defence R & D
Programme (Crowder et al., 2004a,i). The EA used three software packages to carry
out these tests: ISIS, MIKE 11 and HEC-RAS. A total of twelve test specifications were
published specificaliy to provide methods for assessing numerical accuracy, capability,
and reproducibility of the models. The tests are well documented and input/output
datasets along with analytical solutions were made available to the research community.
After careful consideration, it was decided to apply CanalFlow, the canal irrigation
network simulation model outlined in this thesis, to this sequence of test cases to allow
a ready-made comparison with the ISIS, MIKE 11 and HEC-RAS models. As has been
indicated, CanalFlow delivers the first of a new category of irrigation canal network
simulation models using modern high-resolution shock-capturing methods. As may be
expected, CanalFlow outperforms ISIS, MIKE 11 and HEC-RAS in almost all cases.
This will be discussed in Chapter 7. However, before doing so, a discussion outlining
the reasons behind the choice of technology platforms (Java and MySQL) is discussed
in Chapter 3, followed by a theoretical background behind the CanalFlow in Chapter
4.
47
CHAPTER 3
3.1 Introduction
48
3.2 Programming Languages
LABoratory), the language invented by Cleve Moler in the late 1970s (MATLAB, n.d.);
GNU Octave (a clone of MATLAB), a free computer program for performing numerical
computations (Octave, n.d.); R Programming Language (sometimes called GNU S), a
mathematical language and environment used for statistical analysis and display (R,
n.d.); Scilab, a scientific software package for numerical computations providing a pow-
erful open computing environment for engineering and scientific applications (Scilab,
n.d.); and DADiSP (Data Analysis and Display, pronounced day-disp), a computer
program which allows to display and manipulate data series, matrices and images in
an environment similar to a spreadsheet (DADiSP, n.d.).
This chapter takes a short tour of history of programming languages, then describing
object-oriented programming and Java. Further, the chapter has a section on Java's
performance in scientific computing and a rationale behind selecting Java. It also
discusses relational databases, their role and the reasons behind using MySQL.
Computers are very powerful tools that can store, organise, and process a tremendous
amount of data. However, detailed set of instructions are required to perform the given
task. These instructions are in a series of binary coded instructions, called machine
language. In early days, software engineers used to write their programs using a simple
language, called assembly language. In assembly language, programmers had to manu-
ally translate, or assemble, each line into a machine code. This type of translation was a
difficult and tedious task until one software engineer wrote a program, called an assem-
bler. Later, programmers used assembly language to write programs and then used an
assembler to translate the program into machine language. In this section, the history
of programming languages is discussed, in short, before dealing with object-oriented
languages and Java.
With advent of technology, programming languages became more and more con-
venient for programmers to use. Over the years a series of high-level programming
languages have been devised. These languages are attempts to let programmers write
49
3.2 Programming Languages
in something that is easy for them to understand and that is also precise and simple
enough for computers to process it. Early high-level languages were designed to han-
dle specific types of applications. For example FORTRAN was designed for number
crunching; COBOL, for writing business reports; and PASCAL, for student use. Table
3.1 shows classification of some of the more popular high-order programming languages
in generations arranged according to the language features they first introduced (Booch,
1994).
Language Features
First-generation languages (1954-1958)
FORTRAN I Mathematical expressions
ALGOL 58 Mathematical expressions
Flowmatic Mathematical expressions
IPLV Mathematical expressions
Second-generation languages (1959-1961)
FORTRAN II Subroutines, separate compilation
ALGOL 60 Block structure, data types
COBOL Data description, file handling
Lisp List processing, pointers, garbage collection
Third-generation languages (1962-1970)
PL/I FORTRAN + ALGOL + COBOL
ALGOL 68 lligorous successor to ALGOL 60
Pascal Simple successor to ALGOL 60
Simula Classes, data abstraction
The generation gap (1970-1980)
Many different languages were invented
50
3.2 Programming Languages
In 1957, the first of the major languages appeared in the form of FORTRAN (FOR-
mula TRANslating system). The language was designed at IBM for scientific comput-
ing. The components were very simple, and provided the programmer with a low-level
access to the computers innards. The commands like IF, DO, and GOTO statements,
at the time, were a big step forward. The basic types of data in use today had their
start in FORTRAN, these included logical variables (TRUE or FALSE), and integer,
real, and double-precision numbers.
Although FORTAN was good at handling numbers, it was not so good at handling
input and output, which mattered most to business computing. Business computing
started to take ofi'in 1959, and because of this, COBOL (COmmon Business Oriented
Language) was developed. Its only data types were numbers and strings of text. It also
allowed for these to be grouped into arrays and records, so that data could be tracked
and organized better.
In 1958, John McCarthy of Massachusetts Institute of Technology (MIT) created
the LISt Processing (or LISP) language. It was designed for Artificial Intelligence
(AI) research. During the same time, Algol language was created by a committee for
scientific use. Its major contribution is as a precursor to such languages as Pascal, C,
C++, and Java. It was also the first language with a formal grammar, known as Backus-
Naur Form or BNF. Pascal was begun in 1968 by Niklaus Wirth. Its development was
mainly out of necessity for a good teaching tool. Pascal introduced dynamic memory
allocation, which meant that memory could be allocated while a program was being
run, through the NEW and DISPOSE commands. Later, C was developed in 1972
by Dennis Rltchie while working at Bell Laboratory in New Jersey. The transition in
usage from the first major languages to the major languages of today occurred with
the transition between Pascal and C. Its direct ancestors are Band BCPL, but its
similarities to Pascal are quite obvious.
Simula was developed in the 1960s at the Norwegian Computing Centre in Oslo,
primarily by Ole-Johan Dahl and Kristen Nygaard. Syntactically, it is a superset of
Algol 60, adding features that are close to the modern idea of classes and objects,
plus coroutines. Simula introduced the object-oriented programming paradigm and
51
3.3 Object Oriented Programming
thus can be considered the first object-oriented programming language and a prede-
cessor to Smalltalk, C++, Java, and all modern class-based object-oriented languages.
In 1970s, Smalltalk was developed at Xerox Palo Alto Research Centre (PARC) by
by Alan Kay, Dan Ingalls, Ted Kaehler, Adele Goldberg, and others. Smalltalk is a
dynamically-typed object-oriented programming language that has been a great influ-
ence on the development of many other computer languages, including Objective-C,
Actor, Java and Ruby. Many software development ideas of the 1980s and 1990s
like Model-View-Controller, Class-Responsibility-Collaboration card, design patterns
(as applied to software), eXtreme Programming (XP) and refactoring came from the
Smalltalk community.
Objects are pieces of data and instructions that can be packaged and manipulated
by the programmer. Bjarne Stroustroup liked this method and developed extensions to
C known as "C with Classes". This set of extensions developed into the full-featured
language C++, which was released in 1983. C++ was designed to organize the raw
power of C using OOP, but maintain the speed of C and be able to run on many
different types of computers. Sun Microsystems developed Java in 1995 which became
very popular. Microsoft presented the Visual Basic, VC++, and C# programming
languages.
A brief history of programming languages is available at (http://www .levenez.
com/lang/) and O'Reilly has published a timeline ofprograrnming languages (O'Reilly,
n.d.).
Currently, there are over 2000 different high-order programming languages (Booch,
1994). These languages were shaped by the particular requirements of its perceived
problem domain. Furthermore, the existence of each new language enabled developers
to move on to more and more complex problems. With each previously unexplored
application, language desiguers learned new lessons that changed their basic assump-
tions about what was important in a language and what was not. This evolution of
52
3.3 Object Oriented Programming
languages was also heavily influenced by progress in the theory of computing, which
has led to a formal understanding of the semantics of statements, modules, abstract
data types, and processes.
Most programming languages allow their user to quickly create programs that give
satisfying results at the expense of rigorous structure. Languages such as C, Perl and
FORTRAN have been used to create some of the most complex systems. The programs
written in these "structured programming languages" is sometimes referred as spaghetti
code because of their size and complexity (Naughton, 1996). With the increase in the
capacity of computers, programs grew in size and the software to run them became
geometrically more complex and harder to manage. To handle this complexity a new
style of programming was required and thus, concept of an object-oriented programming
was born. This paradigm shift from functional structured programming to object-
oriented programming began almost 30 years ago with a language called Simula 67.
As shown in Table 3.1, programming languages may be grouped into four gener-
ations, according to whether they support mathematic, algorithmic, data, or object-
oriented abstractions. The most recent advances in programming languages have been
due to the influence of the object-model. A language is considered object-based if
it directly supports data abstraction and classes. An object-oriented language is one
that is object-based, but also provides support for inheritance and polymorphism. The
common ancestor of almost every contemporary object-based and object-oriented pro-
gramming language is Simula, developed in the 1960s. Thus, to be truly considered
object-oriented, a programming language should support at the minimum four charac-
teristics:
• Inheritance: can define new classes and behavior based on existing classes to
obtain code re-use and code organization; and
53
3.3 Object Oriented Programming
• Dynamic binding: objects could come from anywhere, possibly across the net·
work. The message is passed on to objects without having to know their specific
type at the time of writing the code. Dynamic binding provides maximum flexi-
bility while a program is executing.
These concepts are not discussed in detail here as the description of these terms can
be found in Booch (1994); Naughton (1996); Schach (2004).
OOP helps writers and programmers use real-world paradigms. 00 analysis is a
methodology used to study real·world systems by identifying systems components, com-
ponentbehaviour, and the relationship between components. This approach provides
tools for programmers to represent elements in the problem space. The representation
is general enough that the programmer is not constrained to any particular type of
problem. Elements in problem space and their representations in solution space are
referred to as 'objects' (Eckel, 2002). The 00 concepts support the development of
robust and scaleable programs because of mechanisms like encapsulation, inheritance,
and polymorphism. Well-designed hierarchy of classes allows reusability of the code
and polymorphism helps to create clean, sensible, readable and future-proof code.
Most numerical codes are written in procedural type languages (such as FORTRAN
and C). This involves writing a large array of procedures which perform operations on
matrices and vectors. Each procedure has to be provided with the correct data, which
usually requires a knowledge of not just what the procedure does but how it is done. As
codes increase in complexity and size, such procedures tend to become deeply nested
and a knowledge of their working is more difficult to acquire. This hampers not only
productivity but heightens the probability of erroneous coding.
The object-oriented programming philosophy does not suffer from the above de-
ficiencies and thus has a significant advantage as compared to procedural systems in
terms of maintaining programming productivity and maintainability. The benefits of
OOP include better data abstraction, improved modularity, greater maintainability, a
54
3.3 Object Oriented Programming
higher degree of code reuse, and enhanced robustness in the face of changing require-
ments (Veldhuizen, 2000). Significant advantage of OOP is it allows rapid software
development through component sharing and code reuse (Kolditz and Bauer, 2004).
Johnson and Foote (1988) have discussed in detail the 00 concepts and reusability of
the software.
In the past decade, the object-oriented technology has been widely used in scien-
tific computations. Early use of OOP was in developing graphical applications and data
visuaiisation tools. Use of OOP in numerical analysis was started in early 1990s focus-
ing mainly on the algorithmic side (Kong and Chen, 1995; Mackie, 1992; Ross et al.,
1992a,b). Its importance to software is analogous to that of the integrated circuit to
hardware. Liu et al. (1996) have used the hierarchical formulation of OOP to present
a generic formula that can be reformulated into the respective formulas suitable for
FEM, FVM and LSFEM (Least-squares FEM). This generic formula serves as a base
formula which in turn could be inherited by those methods in their own reformula-
tion. Malan and LewIs (2004) used C++ to develop algorithms for efficient explicit
edge-based finite-volume code. Kolditz and Bauer (2004) presented the design and im-
plementation of an 00 approach to the computation of multi-field systems. This is a
good example of a 'process-oriented approach' in which the physical processes govern-
ing reality are also the objects of the program structure. Processes are implemented as
objects organised in a class containing all information and data required to treat the
connected physical problem.
Computational hydraulics, bydrologic, environmental and hydraulic modelling as a
whole, are mature disciplines, with hundreds of algorithms in use. It is usual practice
to use the same algorithm to solve different set of problems but still the high degree
of code reuse is normally not achieved. Solomatine (1996) discussed advantages of rep-
resenting a hydraulic process in an object (having data and methods). This paper is
one of the few attempts to present the use of 00 in an architecture of the modelling
systems. Kutija (1998) in her paper challenged the common perception of inability to
code computational hydraulics algorithms in an 00 way. The principal argument was
to reassess the computational hydraulics algorithms with 00 ideas in mind. This was
55
3.3 Object Oriented Programming
the first instance where an 00 numerical algorithm was developed for an open-channel
network. Adriana et al. (2004) provide conceptual micro architectures for flexible hy-
drologic modelling on the basis of a general conceptual architecture that facilitates a
simulation model's integration with environmental information systems. Alfredsen and
Saether (2000) implemented an 00 flood routing, production and impact assessment
model using C++. The design is based on reusable components and a framework to
provide a foundation for future developments. Beinhorn and Kolditz (2005) used OOP
for groundwater modelling and Hatakeyama et al. (1998) proposed and applied a new
00 modelling and programming paradigm to the fluid flow simulation system. Lal
et al. (2005) used 00 with Extensible Markup Language (XML) and GIS to model
complex regional hydrologic system, and organise and present large amounts of com-
plex data. Further, McKinney and Cai (2002) used 00 method to link GIS and water
resources model. Murray (2003) has discussed how 00 techniques can be used to en-
hance and bring real benefits to hydroinformatics models' cores and algorithms; Kutija
and Murray (2002) presented a fast hydrodynamic network modelling system using
object-oriented numerics (discussed in the Section 2.5.14 on page 39); Tisdale (1996)
has presented 00 analysis conducted for South Florida hydrologic system; Wang et al.
(2005) used 00 approach to develop a new watershed based hydrological model; and
Yang et al. (2002) integrated 00 methodology with a 1-D river model.
3.3.2 Java
56
3.3 Object Oriented Programming
57
3.3 Object Oriented Programming
pointers, pointer arithmetic, mallac, and free do not exist. Automatic garbage collec-
tion is an integral part of Java and its run-time system. While Java technology has a
new operator to allocate memory for objects, there is no explicit free function. Once
you have allocated an object, the run-time system keeps track of the object's status
and automatically reclaims memory when objects are no longer in use, thus freeing
memory for future use. This feature makes programming tasks easier, eliminates en-
tire classes of bugs, and in general provide better performance than obtained through
explicit memory management (Gosling and McGilton, 1996).
Java is portable at both the source and object format levels. The most important fe&-
ture Java has to offer is pervasiveness (runs on virtually every platform). But this very
feature caused early interpreted versions of Java Virtual Machine (JVM) to perform
poorly, especially when compared with compiled languages like FORTRAN and C. The
idea that results produced on every JVM should be bitwise identical (Gosling et al.,
2005) on all platforms threatens the usability of Java for high performance scientific
computing. Java does not have a complex data type, although this is not a fatal flaw
since new types are easy to define (Boisvert et al., 1998). However, since Java does not
support operator overloading, these types cannot be made to behave like the primitive
types (float or double). Therefore, performance and convenience is a concern in use of
Java for numerical computations.
Inspite of all these reasons, the combination of Java programming features, per-
vasiveness, and performance could make Java the language of choice for numerical
computing (Boisvert et al., 2001). Indeed, a community of scientists and engineers
developing new applications in Java has been increasing slowly. Things have changed
in the new versions of Java. Today, JVM uses just-in-time (JIT) compiler technology.
JITs operate as a part of the JVM, compiling Java bytecode into native machine code
at runtime. Once the JVM generates the machine code, it executes it at raw machine
speed. Modern JITs perform sophisticated optimisations, such as array bounds check
elimination, method devirtualisation, and stack allocation of objects (Boisvert et al.,
58
3.3 Object Oriented Programming
59
3.4 Relational Database and their Application
Earlier references of use of Java in the development of numerical libraries are found
in Smith et al. (1997) and Russell et al. (1997). Few attempts were made to convert
existing FORTRAN-based numerical libraries to Java (Casanova et al., 1997; Fox et al.,
1997). Blount and Chatterjee (1999) developed and evaluated JLAPACK - the Java
version of the high-performance FORTRAN 77 library, Linear Algebra PACKage (LA-
PACK). Boisvert et al. (1998) describes the Jaya Numerical Toolkit (JNT), a set of
basic numerical functions and kernels developed using Java. This paper reviews ad-
vantages and disadvantages of using Java in scientific computation, and various ways
to overcome these disadvantages. Moreira et al. (2000) developed a high-performance
numerical library, written entirely in Java. The paper has discussed, in depth, program-
ming techniques that can lead to Java numerical code with performance comparable to
FORTRAN or C. The paper demonstrated Java performance in the range of 55-90% of
the best FORTRAN performance for a variety of benchmarks for a high-performance
numerical computing. Hu et al. (2000) concluded that in comparison with Fortran and
C, the I/O and computing performance of Java varies from 30% slower to about 3 times
slower, depending on the platforms and the compilers (and the Java Virtual Machine).
The best relative performance was achieved on a Pentium II, where it was found that
the IBM Java (Version 1.3) yields code that is about 30-40% slower in computing and
I/O performance. But, given the special advantages of Java, and the continuing in-
crease in hardware speeds, this slightly inferior performance of Java, on balance, would
be acceptable for many applications.
60
3.4 Relational Database and their Application
benefits of using a RDMBS come from its ability to store data in a normalised format,
a concept originally presented by Codd (Codd, 1970), who mathematically developed
the relational model to provide a better structure for databases. Data normalization is
simply a way of organising data so that it allows for increased efficiency of data storage
and retrieval. While spreadsheets can store data in a normalised format, it is very
difficult to retrieve in a simple and timely manner.
Codd's twelve rules call for a language that can be used to define, manipulate,
and query the data in the database, expressed as a string of characters. The language,
Structured Query Language (SQL), originally developed in the research division ofIBM,
has been adopted by all major relational database vendors. SQL has been adopted as
an ANSI and ISO standard and it is the most popular computer language used to
create, modify and retrieve data from relational database management systems.
The use of RDBMS has many advantages over using spreadsheets or flat/text files.
RDBMS provide data-program independence, where the method of storing the data,
the order of the stored information, and how the data is managed on disk are inde-
pendent of the software that accesses it. Managing complex relationships is difficult
in a spreadsheet or text file. For example relation between a canal and its reaches
could be maintained by unique reach keys. Spreadsheets or text files do not work well
when there are associations or relationships between stored data items. In contrast,
RDBMSs are designed to manage complex relational data. RDBMS also permit multi-
user transactions. Medium and large scale databases include features that control the
writing of data by multiple users in a methodical way. Additional benefits of databases
include fast and efficient searching, data security, and data recovery support. With a
multi-user support it would be quite easy to convert the existing stand alone model
into a web-based application. Further, the database schema can be used to define an
XML data store and data can be exported directly to the XML file.
There is no known previous reference for a hydraulic model which uses relational
database. Few references have been found for application of relational databases in
environmental and watershed modelling. These instances have primarily used data-
bases in hydrological modelling and real-time data collection systems. Carleton et al.
61
3.4 Relational Database and their Application
3.4.1 MySQL
MySQL is the most popular open source SQL database, developed and provided by
MySQL AB . MySQL is an open source relational database system and is relatively
fast, reliable, and easy to use (http://wml.mysql.com). It is a client/server system
that consists of a multi-threaded SQL server that supports different backends, sev-
eral different client programs and libraries, administrative tools, and a wide range of
application programming interfaces (APIs) (Widenius et al., 2002).
MySQL offers users a wide variety of choices and options in terms of table types, op-
erating system and hardware. Figure 3.2 shows a logical view of MySQL's architecture.
It is easy to install and the configuration file can be changed easily to run on a wide
range of hardware. At the application development level, a variety of data types can
be chosen while creating tables to store records. Different types of tables can be mixed
and matched in the same database. MySQL offers an option of using MyISAM (de-
fault), InnoDB, HEAP, and NDB (Widenius et al., 2002) tables. CanaIFlow database
uses InnoDB and MyISAM table types. Additionally, MySQL is ACID (Atomicity,
Consistency, Isolation, and Durability) compliant. These satisfy four tightly related
criteria that are required in a well-behaved transaction processing system (Zawodny
and Balling, 2004).
A disadvantage of a database approach is that it lacks buffering and hence can slow
down the system especially while inserting records. This problem can be overcome
62
3.5 Concurrent Versions System
.
,,:"'::<;',,:~onnection management. securi~:
I' ". .,/.... . ..
~,-----/---
, . .
-- ----~- ~-- -- . ,
,,
--,
--------------------------------------------~
Figure 3.2: A logical view of MySQL's architecture. The topmost layer shown in the
figure deals with the network-based client/server services like connection handling, au-
thentication and security. The functionalities like query parsing, analysis, optimisation,
caching, stored procedures and all the built-in functions reside in the second layer. Fi-
nally, the third layer is made up of storage engines, responsible for the storage and
retrieval of all data stored in MySQL.
Concurrent Versions System (CVS) is a version control system and by using it one can
record the history of source files. For example, bugs sometimes creep in when software
is modified, and it might not get detected until a long time after the modifications are
63
3.5 Concurrent Versions System
done. CVS can easily retrieve old versions to see exactly which change caused the bug
and if necessary revert back to the previous version. This is of big help especially in a
research project of this kind. CVS stores all the versions of a file in a single file in a
clever way that only stores the differences between the versions and not the different
versions itself. It also helps if a group of people are working on the same project. In
a group, it is all too. easy to overwrite each others changes unless one is extremely
careful. CVS solves this problem by insulating the different developers from each other
as every developer works in his/her own directory, and CVS merges the work when the
developer is done. This is just a short introduction to the CVS and more details on
its use and features could be found on the CVS home http://www . nongnu. org/ evs/,
in the CVS manual http://ximbiot . eom/evs/manual/ and in Fogel and Bar (2001).
This section is an attempt to share, with other research community, some of the software
development practices using CVS.
To start with, a CVS repository was set-up on a departmental web server. At the
end of the day (or as and when required) the code was committed to this repository.
This ensured that the latest version of the code was always placed on the CVS server
and there was absolutely no need to take backups on the local diskl. This saved a
lot of development time and was more organised approach toward the development.
CVS was extremely helpful during code reviews with the Supervisor. It ensured that
the latest copy of the code was reviewed and then where possible minor changes were
done immediately to the code and changes committed to the repository. In many in-
stances, Supervisor reviewed the code by checking out the code from the repository,
making changes if necessary and immediately updating the repository for further de-
velopment. Best part of this was all these actions were performed remotely without
meeting the researcher. This ascertained timely reviews, corrections and thus saving
total development time.
1It is always advisable to back-up the CVS repository itself at regular intervals
64
3.6 Concluding Remarks
As described in this chapter technology platforms selected for the software development
were Java and MySQL with a CVS backbone. Java was an obvious choice because of
its simple 00 structure, platform independence and its flexibility in converting the
existing model in a web-based application in future. The section 3.3.2.1 has discussed
advantages and disadvantages of use of Java for the numerical computing. The major
disadvantage is its speed and we decided to compromise its lesser speed for all other
advantages it offered over other programming languages. The consideration was done
because of the fact that Sun is totally committed to bringing new and better JDKs in
future and also because the advances in hardware is not going to limit the computing
power.
Right from the beginning our aim was to design an appropriate data structure
to store and retrieve the data required for hydraulic models. MySQL was primarily
considered because of its speed, low memory footprint and its suitability for small to
medium size databases. It is also the most popular open source database.
The design and development procedure of the software is discussed in detail in the
Chapter 5. This chapter also discusses the design of relational database model used in
CanalFlow while actual structure of tables is presented in the Appendix A. Before going
into details of the software development, Chapter 4 lays out the theoretical background
of the CanalFlow software.
65
CHAPTER 4
CanalFlow - Theory
4.1 Introduction
66
4.2 One-Dimensional Unsteady Open Channel Flow
The flow in open channels and in other bodies of water is classifled as steady or un-
steady flow. The flow of water in river channels, canals, reservoirs, galleries, tunnels
and culverts, for which velocities change with time, is deflned as unsteady flow (non-
permanent, non-stationary, or time-variable free-surface water flow). Water flow in a
natural channel is almost always unsteady (Yevjevich, 1975).
The equations which solve the unsteady open channel flow are based on the St
Venant hypotheses (Cunge et al., 1980, p.8). The one-dimensional unsteady flow in
channels, assuming that the density is constant, can be described by two dependent
variables; the cross-sectional area A and the discharge Q. These dependent variables
define the state of the fluid motion along the channel with time, i.e., as a function of
two independent variables (x for space and t for time). We will now look at the I-D,
unsteady flow equations in an integral, differential and matrix form in the next few
sections.
Consider the control volume (CV) in the (x, t) plane between cross sections x = Xl and
x = X2, and between times t = tl and t = t2 shown in Figure 4.1. The flow is assumed
to be nearly horizontal, i.e., the angle", between the channel bed and the x-axis is
sufficiently small so that cos Cl< = 1. Based on these hypotheses, the basic equations can
be formulated by using the principles of conservation of mass and momentum within a
control volume. The following derivation has been adapted from (Cunge et al., 1980)
and is presented here for completeness.
The net inflow of mass into the volume is defined by the time integral of the differ-
ence between the mass flow rates entering (puA)x, and leaving (puA)." the control
volume:
(4.1)
67
4.2 One-Dimensional Unsteady Open Channel Flow
Figure 4.1: Longitudinal section view of a control volume showing external forces acting
on a control volume in a x-direction - pressure forces F;. and F;. at boundaries Xl and
X2 respectively, gravitational force Fg which is weight component along the channel,
and frictional resistance force Ff manifested by means of shear along the bottom and
sides of a channel.
This net inflow must be equal to the change of storage in the reach during the time
interval:
(4.2)
(4.3)
where Q = uA
68
4.2 One-Dimensional Unsteady Open Channel Flow
The conservation of momentum in the x-direction requires that the change of mo-
mentum in the control volume between times tl and t2 be equal to the sum of the net
inflow of momentum into the control volume and the time integral of the sum of the
external forces acting on it over the same interval.
• Momentum is the product of mass and velocity, and momentum flux through the
flow section is the product of the mass flow rate and velocity, or
The net momentum flux into the control volume (momentum entering through
section x = Xl minus the momentum leaving through the section x = X2) at tune
is
(4.5)
• Let us assume that the external forces acting upon the control volume in the
x-direction are pressure, gravity, and frictional resistance.
- The pressure force Fp" is the difference of pressure forces F;. and F;.,
applied at boundaries XI, X2 of the reach. At any cross section x with free
surface elevation y(x) the pressure force is expressed under the hydrostatic
distribution hypothesis by
rh(x)
F;. = g io p [h(x) - '7J a(x, '7) d'7 (4.7)
69
4.2 One-Dimensional Unsteady Open Channel Flow
A)
I,
rb
..--c.-+-------------f-- --
yb x
z
B) ,.
F'
I
F' ,. F"
"
0 • ~
x
, , ,, F"
"
70
4.2 One-Dimensional Unsteady Open Channel Flow
Thus the time integral of the net pressure force Fp" when F~. is expressed
as in Equation (4.7), is
(4.8)
For an infinitesimal channel length dx, the increase of the pressure force due
to the width variation is represented by the increase of the wetted area da d1]
for constant h = ho times the distance of its centroid from the free surface
or
(4.9)
where I2 = Joh(x) (h - rJ) [~lh=ho drJ, a pressure force term due to variation
in a longitudinal width. Equation (4.9) is not valid if a sudden width change
occurs between sections Xl and X2.
71
4.2 One-Dimensional Unsteady Open Channel Flow
- In that case supplementary forces will act upon the control volume, and they
must be taken into account (see Figure 4.2). In any case, in such a situation
the curvature of streamlines will not be negligible, violating one of the basic
hypotheses in stating that the streamline curvature is small and vertical
accelerations are negligible (hence pressure distribution is hydrostatic).
- The force Fg due to gravity, that is to say, the weight component along
the channel axis, is evaluated by assuming that the channel bottom slope
So = -~ = tana is small (Yb being bottom elevation above datum), so
that tan a "" sin a:
It,
h
Fgdt = It' lX'
tl Xl
pgA Sodxdt (4.10)
(4.11)
(4.12)
where f>M is the net increase in momentum contained in control volume and M f
is the net momentum contained in the control volume. For constant density p
72
4.2 One-Dimensional Unsteady Open Channel Flow
Xl
' plzdxdt + 9 t'l
Jtl Xl
x
'A (So - Sf) dxdt (4.13)
Equation (4.13) is the integral form of the momentum conservation equation for un-
steady one-dimensional flow in natural channels of an arbitrary shape. Equations (4.3)
and (4.13) together are the integral form of the unsteady flow relations based on the
St Venant hypothesis.
Equations (4.3) and (4.13) are integral relations established without the requirement
that the dependent flow variables, viz., A, Q, u, be continuous or differentiable. Also, in
the above section, the distance X2 - Xl is not considered to be infinitesimally small. The
differential equations of gradually varied unsteady flow may be obtained from integral
equations if one assumes that the dependent variables are continuous, differentiable
functions (Cunge et al., 1980). Then, by Taylor series expansions:
&A &2 A At2
(A)t, = (A)tl + at At + &t2 "'2 + ...
(4.14)
By retaining only first derivatives and assuming that Ax and At approach zero,
(4.15)
73
4.2 One-Dimensional Unsteady Open Channel Flow
Substitution of the first terms of expansions (4.17) into Equation (4.13) and then pas-
sage to the limit (D.x -> 0, D. t -> 0) leads to
(4.18)
If the relations (4.16) and (4.18) are to hold everywhere in the (x, t) plane, then they
must hold for an infinitely small volume, and we can write two differential equations:
Continuity equation:
8A 8Q
.-+-.
at ' 8-=D
x:
(4.19)
Momentum equation:
(4.20)
74
4.2 One-Dimensional Unsteady Open Channel Flow
Equations (4.19) and (4.21) are written in a special form, often called the 'divergent'
form of partial differential equations (also called as conservation form). If the right
hand sides of Equations (4.19) and (4.21) are equal to zero, these equations express
nil divergence of the mass and momentum vector functions in any closed contour in
the (x, t) plane; and mass and momentum are conserved. When the right hand side
of Equations (4.21) is different from zero, momentum is no longer conserved, the free
terms acting as momentum sources or sinks.
h and h given in the Equation 4.8 and 4.9, respectively, can be derived for the
trapezoidal and rectangular channel cross-sections. Referring to the Figure 4.3,
hhl~(-- B
"""t
T '( ---~)>lI/ h -IJ
h ....L
IJ
L T
Figure 4.3: Illustration of the Trapezoidal channel cross-section to calculate the hy-
drostatic pressure force term (h) and the pressure force term due to variation in a
longitudinal width (12). B is the channel bottom width, ml and m2 are left and right
side slopes, respectively and (J is the cross-sectional width.
TJ
= and TJ
=
h (4.22)
h mlh h m2h
75
4.2 One-Dimensional Unsteady Open Channel Flow
(4.24)
(4.25)
(4.27)
The flow is expressed by two dependent variables; cross-sectional area A and the dis-
charge Q, thus, the St Venant equations for one-dimensional unsteady flow in a non-
prismatic channel of arbitrary cross-section, can be written as:
76
4.2 One-Dimensional Unsteady Open Channel Flow
where U is the vector of conserved variables, F(U) is the vector of fluxes and each of
its components is a function of the components of U such as,
U= [~]
F(U) = [g:
A
Q
+gh
]
(4.29)
and
S(U) = [g A(8 0
- ~/) + 9IJ
where t is time, x is the horizontal distance along the channel, A is the wetted cross-
sectional area, Q is the flow discharge, 9 is the gravitational acceleration, and 8 0 is
the bed slope. The frictional slope 8, is defined by either Manning's or Chezy's rela-
tion (Cunge et al., 1980),
• Manning's relation:
(4.30)
• Chezy's relation:
(4.31)
OU A(U)oU =S A= of (4.32)
at + ox ' OU
where A(U) is the matrix
A(U) = ;~ = [~ ~] 8q eq
(4.33)
77
4.2 One-Dimensional Unsteady Open Channel Flow
called the Jacobian matrix corresponding to the flux F(U). The entries of A(U) are
partial derivatives of the components of the vector F with respect to the components
of the vectors of conserved variables U. U can be expressed in terms of components
ql and q2
A(U) _ [ 0 1] (4.36)
- _u 2 + !gh 2u
where u = Q/A
(4.37)
where I is the identity matrix. Expanding Equation (4.38) we obtain the polynomial
which are real and distinct. Physically, these eigenvalues represent speeds of propaga,-
tion of information. The speeds with positive sign are in the direction of increasing x
and negative otherwise.
A right eigenvector of a matrix A corresponding to an eigenvalue Ai of A is a vector
Ri = [rj, r~l satisfying ARi = AiRi. In a similar way, a left eigenvector of a matrix A
78
4.2 One-Dimensional Unsteady Open Channel Flow
RI = [ 1 ] R2 = [ 1 ] (4.40)
u-c' u+c
The system of equations 4.~l2 is hyperbolic at a point (x, t) as A has two real
eigenvalues Al and A2 and a corresponding set of linearly independent right eigenvectors
RI and R 2. Further, it is strictly hyperbolic as the eigenvalues Al and A2 are distinct.
Region inflrnCCd by Q
, /
P
C+ C-
C
V
C+
XL A!
n
I~XR
Q x " x
Initial data upon which P depends
(A) (B)
Figure 4.4: Region of influence and domain of dependence for the solution of hyperbolic
SWE. (A) Illustrates the range of influence of a point Q with characteristic speed C+
and C_, and (B) Illustrates the domain of determinacy of an interval [XL, XRJ with
characteristic speed C+ and C_
channel flow propagates along the channel in time and in two directions: downstream
79
4.3 Riemann Problem
and upstream. For a 1-D, unsteady flow, consider a point Q in (x - t) plane in Figure
4A(A). The perturbation which takes place at point Q influences the shaded region
limited by curves C+, C_ which represent the paths of disturbance propagation. This
region above point Q is called as region of influence. Considering a point P, referring
Figure 4.4(B), one can define backward in time the domain in which disturbances can
influence the conditions at point P. These two characteristics extended backward to
the x-axis are intercepted by the portion on x-axis, XL and XR. Assume that the
dependent variables A and Q are known on the x-axis at t = 0 then the solution can
be obtained by "marching forward" in the time space starting from the given boundary.
The solution for A and Q at point P will depend only on that part of the boundary
between XL and x R and included between the two retreating characteristics. For this
reason, the region below point P is called as the domain of determinacy and the domain
of dependence is an interval [XL, XRJ on the x-axis that is subtended by characteristics
passing through the point P. The projection of the position of these small discontinuities
into (x - t) plane produces lines cailed characteristics (Cunge et al., 1980). The C+
curve is called as a right-going and C_ as a left-going characteristic.
Riemann problem is the initial value problem (IVP) for the hyperbolic conservation
laws with the simple initial condition consisting of two constant states. The Riemann
problem for the shallow water equations is a generailsation of the dam-break prob-
lem (Toro, 1997,2001).
PDE. BU BF(U) = 0
. Bt + Bx
(4.41)
lC: U(x, 0) = {UL if x < 0,
UR if x > o.
Considering Equation 4.28 and 4.29, the initial states U Land U R with,
(4.42)
80
4.3 Riemann Problem
are constant vectors and represent conditions at time t = 0, to the left of x = 0 and
to the right of x = 0, respectively. The possible wave patterns emerging from the
solution of Riemann problem given in Equation 4.41 is illustrated in Figure 4.6 and the
structure of the Riemann problem is as shown in Figure 4.5. If Ui-I and Ui are the
cell averages in two neighbouring grid cells on a finite volume grid, then by solving the
Riemann problem with UL = Ui-I and UR = Ui, one can obtain a numerical fiux and
update the cell averages over a time step. For hyperbolic problems the solution to the
Riemann problem is typically a similarity solution, a function of xlt in an expansion or
rarefaction wave region, and consists of a finite set of waves that propagate away from
the origin with constant wave speeds (Leveque, 2002).
Figure 4.5 shows that, from origin (0,0) in the (x, t) plane there are two waves
travelling with speeds equal to the characteristic speeds AI and A2 where (AI < A2).
The solution to the left of dxldt = AI is simply the data state UL = alRI + a2R2 and
to the right of dxldt = A2 the solution is the constant data state U R = ,8IRI + ,82R2
where a and ,8 are constant coefficients. The wedge between AI and A2 is called as the
star region and its value is due to the passage of two waves emerging from the origin of
the initial discontinuity (Toro, 1997). The characteristics with speeds AI and A2 can be
traced back from the point P*(x, t), which basically are parallel to those passing from
the origin and they pass through the initial points xa = x - A2 t and xA = x - AI t.
,
Referring to Figure 4.7 and
• By selecting a time t* and a point XL to the left of the slowest wave such that
U(XL, t') = UL, the solution at the starting point (XL, to) is
2
UL(X,t) = 2:aiRi (4.43)
i=l
All coefficients are a's and the point (XL, to) lies to the left of every wave .
• By moving to the right of (XL, t') on the horizontal line t = to, the point cross the
wave dxldt = At, hence x- Alt changes from negative to positive, thus, changing
the coefficient from al to ,81. The solution of the star region is given by
(4.44)
81
4.3 Riemann Problem
~ =u-c
Figure 4.5: Structure of the general solution of the Riemann problem for one-
dimensional shallow water equations. The two wave families separating three constant
states denoted, from left to right, by UL, U. and UR. The region between the left and
right waves is called the star region.
__________ ~~---------+x
__________~~---------+x
o o
(A) (B)
__________~~----------_+x
o
(D)
Figure 4.6: Possible wave patterns in the solution of the Riemann problem for the
one-dimensional shallow water equations; (A) The left wave is a rarefaction wave and
the right wave is a shock wave; (B) The left wave is a shock and the right wave is a
rarefaction; (C) Both left and right waves are rarefaction waves, and (D) Both left and
right waves are shock waves.
82
4.3 Riemann Problem
~ =u-c
U.
0·. . · · . · · · · · . · ·. ····. · . · ·. · · · ·
;
i, U
;
L
~;--------~t- __________-+x
o
Figure 4.7: The Riemann problem solution found by travelling along dashed horizontal
line t = t'. U L is the left state and U R is the right state separated by U, .
• By moving further right and crossing the ),2, the value of x - ),2t changes from
negative to positive and hence the coefficient changes to f32, thus, solution to the
right of the fastest wave is
(4.45)
From the above equations, the jump in U across the whole wave structure in the solution
of the Riemann problem is
(4.46)
This is an eigenvector expansion with coefficients that are the strengths of the waves
present in the Riemann problem. The wave strength of wave i is f3i - "i and the jump
in U across wave i, denoted by (<lU)i is the Rankine-Hugoniot jump condition and is
given by,
<IF = Si<lU (4.47)
83
4.4 Finite Volume Method - Formulation
This section deals with the explicit FVM of the Godunov type for solving the time-
dependent, non-lineax, St Venant equations in a one space dimension for a canal irri-
gation network. A typical canal irrigation network is shown in the Figure 2.1, page
10. This network is basically a network of chaxmels, connected to each other with or
without a control structure. From a computational perspective the channel looks as
shown in the Figure 4.8. Each chaxmel runs between a fromNode (at from side) and
a toNode (at to side) with or without a control structure at from side and to side.
The solution for a network is an iterative process where at time, t = 0, unsteady flow
equations axe solved for each channel in a network, before incrementing the time to
t = 1. In each time step values axe updated using approximation to the flux through
~ Channel ~
fromNode ••-------------------------------------------.. roNode
i ,,i
1 ,,,
,
1=1
i
I.
L
••
the endpoints of the intervals. The ilh grid cell is denoted by (see Figure 4.9),
84
4.4 Finite Volume Method - Formulation
Now, the basic form of the numerical method is obtained referring to Leveque (2002)
t
U ,Il + 1
F;~112 F;:112
i-I --i ;+1
t"
x
U", U ,':,
Figure 4.9: Illustration of a FVM, shown in the x-t space, for updating the cell average
ut by fluxes, F?_1/2' F?+1/2' at the cell edges, Xi-l/2, Xi+1/2'
and Toro (2001). Recalling equations 4.3 and 4.13, the integral form of the Equation
4.28 without a source term can be written as
F i - 1/ 2 = A
1 i'n+!
F(U(Xi_l/2, t)) dt, Fi+l/2 = A1 l'n+1 F(U(Xi+l/2' t)) dt (4.51)
ut tn L.J.t tn
85
4.4 Finite Volume Method - Formulation
where At == t n+1 - tn. With these definitions the integral form of Equation 4.28 can
be written as,
This is a basic form of the numerical method where Fi+I/2 term is called as the
numerical fiux corresponding to the intercell boundary at x = Xi+!/2 between cells i
andi+l.
Now we look at the steps followed in order to solve the flow equations for a charmel:
ii. Initialise the state variables (A and Q), at t = 0, throughout the domain and
apply boundary conditions at each node as explained in the Section 4.6;
• The state variables are updated using an exact Rlemann solver described in
the Section 4.4.2
• Godunov's update is applied in the next step which is given in the Section
4.4.1.1
• High resolution methods and limiters are used as described in the Section
4.5
• Source term is applied using Equation 4.30 or 4.31
In its simplest form the Godunov type scheme consists of three steps called as reconstruct-
evolve-average (REA):
86
4.4 Finite Volume Method - Formulation
1. Reconstruct a piecewise polynomial function Un(x, t n) defined for all x, from the
cell averages Ur. In the simple case this is a piecewise constant function that
takes the value Ur in the ith grid cell,
2. Evolve the hyperbolic equation exactly (or approximately) with this initial data
to obtain Un(x, t n+1 ) at a time t.t later. In this case, the exact Riemann solver
is used.
3. Average this function over each grid to obtain new cell averages
u n+1 =
_l_JUn(x tn+1)dx (4.53)
t.b..x '
e.
In step 1 a function Un(x, t n ) is reconstructed from the discrete cell average. In
Godunov's original approach this reconstruction is a simple piecewise constant function.
This leads to Riemann problem, but gives only a first order accurate method. To obtain
better accuracy, a better reconstruction can be used, for example a piecewise linear
function that is allowed to have a nonzero slope in the ith grid cell. This idea forms
the basic of the high-resolution methods as seen in the Section 4.5.
So the exact solution at time tn+1 can be constructed by piecing together the Rie-
mann solution, provided that the time step t.t is short enough that the waves from
two adjacent diagram of this process have not yet started to interact. But this type of
scheme requires that
Smaxt.t 1
- t.x <
--""?'-- - -2 (4.54)
where Smax is the maximum wave speed. The quantity smax ~; is Courant, Friedrichs,
and Lewy (CFL) condition and the number is Courant number. It is apparent that the
condition is limited to Courant number less than!.
The first version of Godunov's method defines new average values Ui+ 1 at time
t n+1 = t n + t.t via the integrals given in Equation 4.53. The integrand U(x, t) is an
87
4.4 Finite Volume Method - Formulation
exact solution of the conservation laws. The second version of the Godunov's method
is written as:
U in +1 = uni + Ax
At [F i-I/2 - Fi+1/2 1, (4.55)
with intercell flux given by Fi+I/2 = F(Ui +1/2(0)), if the time step At satisfies
At::; Ax (4.56)
Smax
This is a CFL condition and second version allows a larger time step resulting in a more
efficient time-marching scheme. Here a wave is allowed to travel, at most, a complete
cell length Ax in a time At. Thus,
Smax
At < 1
(4.57)
Ax -
The proof of this form is given in Toro (1997) and hence not discussed here.
The Figure (4.10) shows a linear system of two equations solved assuming Al < 0 < A2.
The function un (x, tn+1) will typically have two discontinuities in the grid cell Ci, at
the points Xi-I/2 + A2At and xi+I/2 + AIAt. As seen in the Section 4.3, for a linear
J1f-112
Figure 4.10: lilustration of the REA algorithm. The Riemann problem is solved at
each cell interface, and the wave structure is used to determine the exact solution time
At later. The wave WLI/2' for example, has moved a distance A2At into the cell.
system the solution to the Riemann problem can be expressed as a set of waves,
2 2
Ui - U i _ 1 = 2: af_l/2 RP == 2: W;_1/2 (4.58)
p=l p=l
88
4.4 Finite Volume Method - Formulation
Consider the wave denoted by WLI/2 for example. It consists of a jump given by
2
Wi-I/2 2 R2
= (Y;i_I/2
propagating at speed A2t.t. This wave modifies the value of V over a fraction of the
grid cell given by A2t.t/t.x. It follows that the effect of this wave on the cell average
of V is to change the average value by the amount
The minus sign arises because the value W~_1/2 measures the jump from right to left.
Each of the waves entering the grid cell has an analogous effect on the cell average, and
the new cell average can be found by simply adding up these independent effects.
The 2-wave originating from Xi-I/2 and the I-wave originating from xi+I/2, based on
the presumed wave speeds, as shown in Figure 4.10, are used. This can be written in
a form that generalises to arbitrary hyperbolic systems of m equations. Let
and suppose the solution of the Rlemann problem consists of 2 waves WP travelling at
speeds Ap, either positive or negative. The cell average can thereby updated by
(4.61)
The cell average is affected only by all right-going waves from Xi-I/2 and only by all
left-going waves from Xi+1/2' In a shorthand notation,
m
A-t.Ui _1/ 2 = I: (AP)-W;_1/2'
p=1
m (4.62)
A+ t.Ui _I / 2 = I: (AP)+W;_1/2'
p=l
89
4.4 Finite Volume Method - Formulation
n+1= un.·
Ui b.t [Jl+
~
i - "",x AV:i-l/2 + Jl"",UH1/2
.. ""'. -A 1 (4.63)
The symbol Jl +b.Ui - 1/ 2 should be interpreted as a single entity that measures the
net effect on the cell average of all right going waves from Xi-l/2, while Jl- b.UH1 / 2
measures the net effect on the cell average of all left-going waves from xHl/2' These net
effects are usually referred as fluctuations. Note that within cell ei , it is the right-going
fluctuation from the left edge, Jl+ b.Ui _ 1 / 2 , and the left-going fluctuation from the right
edge, Jl- b.UH1 / 2, that affect the cell average.
The wave-propagation form (Equation 4.6:3) of Godunov's method can be related
directly to the numerical flux function,
F~_1/2 =
1
b.t L
tn + 1
!(VI(W_l, Vf))dt (4.64)
= !(VI(V:'-l, W))dt
Where V I is the similarity solution to a Riemann problem along x It = 0 when left and
right state is specified. The value of V in the Riemann solution along x = Xi-l/2 is
90
4.4 Finite Volume Method - Formulation
Similarly, there are two ways to express F[tl/2. Choosing the form
m
F[tl/2 = AUi + 2:: (,\P)- W;+1/2
p=1
and combining this with 4.66 in the flux-differencing formula 4.52 gives
since the AUi terms cancel out. This is exactly the same expression obtained in 4.61.
m
Fj:l/2 = !(Ui - I ) + 2:: (,\P)- W;_1/2 == !(Ui-ll + A-flUi_1/2 (4.68)
p=l
or
m
F;'':..1/2 = !(Ui) -2:: (,\P)+ W:'-1/2 == !(U
p=l
i) - A+ flUi _ 1/ 2,
.
(4.69)
corresponding to 4.65 and 4.66 respectively, where the speeds '\p and waves WP are
obtained from the solution to the lliemann problem.
91
4.4 Finite Volume Method - Formulation
As discussed in the Section 4.3 and in Figure 4.5, two waves separate three constant
states in terms of the vector of primitive variables,
these are UL (left data), U. and UR (right data). The constant values of the water
depth and velocity in the star region are denoted by h. and u.. The exact lliemann
solver used in CanalFlow is adapted from the conservation-laws package (CLAWPACK)
software developed at the University of Washington (http://,,,,,, . amath. washington.
edu/-claw/). The theoretical background for this solver is as given in Leveque (2002)
and could be described using following steps:
1. Calculate left state velocity, Ut, right state velocity, Ur, and locate the occurrence
of dry state:
UL = {*"o if AL > Acritica/,
if A L < Acritica/
(4.70)
Acritica/ in this case is taken as IE-20. The dry state appears to the left if UL =0
and I-wave is zero strength shock while 2-wave is a rarefaction.
UR=
~R if AR > Acritica/, (4.71)
{
. - 0 · - - if AR < Acritica/
Dry state appears to the right if U r = 0 and I-wave is rarefaction while 2-wave is
a zero strength shock. If both states are dry then both waves are zero-strength
shock. Wave speed in left state is given by
where hLand hR is depth in left and right state calculated from A L and A R
respectively.
2. Compute intermediate state Urn (urn and hrn): Depending on the specific data,
any combination of shocks and rarefactions are possible, for general values of UL
92
4.4 Finite Volume Method - Formulation
and UR. In general, two functions <PL and <PR are defined, to find the intermediate
state Um (Leveque, 2002, p. 281-282),
if h > hL (shock),
(4.72)
if h < hL(rarefaction)
and
3. Compute wave speeds: The left and right speed of the I-wave and 2-wave is
calculated as:
(a) I-wave
AL =UL - J9hL
AR = (hL UL - h m u m) (4.74)
hL-hm
• Speed AL and AR when I-wave is rarefaction
AL = UL - JghL
(4.75)
AR = Um - Jgh m
(b) 2-wave
93
4.4 Finite Volume Method - Formulation
AL = Urn + Jgh m
(4.77)
AR = UR + JghR
4. Set h., u. to Godunov values along x It = 0: Based on the left and right speed
of I-wave and 2-wave, five configurations are possible for the Riemann problem
between left and a right state. These configurations are as follows:
(a) I-wave transonic rarefaction (Condition: For I-wave, AL < 0 and AR > 0)
1
h. = 99(UL+2 J9 hL)2
(4.78)
u. = UL - 2 (.fih. - J 9 hL)
(b) 2-wave transonic rarefaction (Condition: For 2-wave, AL < 0 and AR > 0)
(4.79)
(c) I-wave fully supersonic and right going (Condition: For I-wave, AL > 0)
(4.80)
U", = UL
(d) 2-wave fully supersonic and left going (Condition: For 2-wave, AR < 0)
(4.81)
(e) I-wave is left going and 2-wave is right going (Condition: For all other cases)
h. = h m
(4.82)
U. =U m
5. Compute the fluxes F and left-going and right-going fluctuations, A±: Fluxes are
given by the following relations and are based on the equations described in the
94
4.5 High Resolution Methods
Section 4.2.3:
FI = ALuL,
FA =ARuR,
F; =A.u.,
FL2 =ALUL+gh
2 2 (B hL) ' (4.83)
L 2" +(ml +m2 )6
2 2
FR=ARUR+gh 2 (B
hR)
R 2"+(m1 +m2)6 '
F; = A. U; + 9 h; (~ + (ml + m2) ~. )
and the fluctuations are given by following equation
A+ = FR - F"
(4.84)
6. Compute waves and their speed: The wave strength 'W" can be deflned as the
total jump across the wave (Leveque, 2002), thus,
(4.85)
where V~ and Vi, are the states just to left and right of the wave. The charac-
teristic speed )..P varies continuously through the rarefaction wave, thus, average
speed is used and is given by
(4.86)
The second-order accurate methods give better accuracy on smooth solutions but fail
near discontinuities, where oscillations are generated (e.g. Leveque, 2002, p.lOl). Even
in the case of smooth solution, oscillations appear due to dispersive nature of these
methods. The upwind methods have the advantage of keeping the solution monoton-
ically varying in regions where the solution should be monotone but the accuracy is
95
4.5 High Resolution Methods
not very good. The idea with high-resolution methods is to obtain second-order ac-
curacy where possible but not to insist on it in the regions where the solution is not
smooth (such results achieved by the high-resolution methods are shown in Leveque,
2002, p.104).
The Godunov's method (discussed in the Section 4.4.1) is first-order accurate and
introduces numerical diffusion, yielding poor accuracy and smeared results (Leveque,
2002). The correction terrns are introduced in the Equation 4.63 to obtain a method
of the form,
(4.87)
where, the fluxesF are based on the waves resulting from the exact Riemann solution
which are calculated in the process of determining the fluctuations A ±. F is given by,
2
F -l/2 =
i ~ L ISf-l/21 (1- ~~ISL/21) WL/2
p=l
(4.88)
W is the pth wave arising in the solution to the Riemann problem and W is a limited
version of this wave.
4.5.1 Limiters
As discussed earlier the second order accurate solutions produce unwanted local maxima
(overshoot) or minima (undershoot). In order to avoid this behaviour, limiters are
used (Erduran et al., 2002). There are two types of limiters; flux limiters and slope
limiters. Their name comes from the fact that they limit either the fluxes themselves or
slope of the Riemann data. In this model flux limiters given below are used (Leveque,
2002):
minmod: </>(8) = minmod(l, 8),
superbee: </>(8) = max(O, mine!, 28), min(2, 8)),
(4.89)
MC: </>(8) = max(O, min((1 + 8)/2,2,28)),
van Leer: </>(8) = (J + 181
1+ 181'
96
4.6 Initial and Boundary Conditions
where the function <1>(0) is the flux limiter junction, whose value depends on 0 which is
a measure of smoothness of the data and MC is monotonised centered. The 0 is defined
in such a way that it takes into account the degree of alignment of the wave vector as
well as their magnitudes. This is accomplished by projecting the vector WL1/2 onto
the vector W;_1/2 to obtain a vector Of-1/2 W;_1/2 that is aligned with W;_1/2 (Leveque,
2002). Thus, 0f-1/2 is given by,
WL1/2· Wl'-t/2
Of-1/2 = PP' (4.90)
Wi_1/2,Wi_1/2
where. is the dot product. The index I is used to represent the interface on the upwind
side of Xi-1/2:
I = {i -
i+1
1 if it > 0
ifit<O
(4.91)
where it is the constant velocity. In the region where the data is smooth 0;'..1/2 "" 1
(except at extrema) and near a discontinuity 07-1/2 is far from 1. Thus, to obtain a
high-resolution method, the wave W;_1/2 is replaced by a limited version in Equation
4.88,
(4.92)
The Equations 4.19 and 4.20 given in the Section 4.2.2 govern the flow in an open-
channel. These are the operation of the laws of physics (law of conservation of mass
and momentum). In addition to these laws, the state of the flow is governed by the
initial conditions (lCs) and boundary conditions (BCs) at the domain boundaries (at
a and b in the Figure 4.11). This section describes the way CanalFlow handles initial
and boundary conditions.
In this process, state variables A and Q are initialised (at t = 0) throughout the domain.
These variables can be initialised in two ways,
97
4.6 Initial and Boundary Conditions
• by specifying h (calculate A based on h) and Q at the start and end of the channel
and then values at each computational nodes are linearly interpolated; or
• by specifying h and Q at each computational node. This is true for the elevation
z as well.
As seen in Figure 4.8, the domain is divided into finite volume cells. The initial con-
ditions are specified for each of this finite volume cells and the solution is progressed
with these initial values.
The computations are carried out on a finite set of grid cells covering a bounded domain.
As the required boundary information is not available in the first and the last cell,
physical boundary conditions are provided that update these boundary cell values based
on the fluxes entering and leaving the domain. Figure 4.8 on page 84, shows such
computational domain of the length L and the boundary interfaces.
where feU) is a cell-centered flux. In the computational domain [a, b] the lliemann
solver provides A+ and A- at all interfaces except Ai/2 and A N+ 1/ 2. These are calcu-
98
4.6 Initial and Boundary Conditions
t"+ 1
1';"
1I A,';, - - -- - - -- AN+~
-~ -
/"
a U, U2 Ut
I(
Figure 4.11: Illustration of the computational domain [a, b] of the length L and cell
average values U. Fl/2 is the flux entering the domain at a, upstream side, and FN+J/2
is the flux leaving the domain at b, downstream side. Also shown is the right and left
going fluctuations A + and A -, respectively.
lated based on the boundary fluxes F{i2 and FN+J/2 (Equation 4.93). The methodology
to estimate these boundary fluxes is now explained.
Boundary fluxes are estimated using the Method of Characteristics (MoC). The St
Venant equation in characteristic form is written as (Strelkoff, 1970),
where Q is discharge per unit time, A is cross-sectional area, h is hydraulic depth de-
fined as ~, T is top-width, So is longitudinal slope and S, is frictional slope while x
and t represent the along-channel distance and time respectively. The forward char-
acteristic (0+) is used at the end of the reach and the backward characteristic (0_)
is used at the beginning as shown in the Figure 4.12. For any boundary, the number
of pieces of information which must be specified relate to the number of character-
istics entering/leaving that region. Consider the upstream boundary. If the flow is
sub-critical, then the 0+ characteristic points into the computational domain, whereas
the 0_ characteristic leaves the region. The result is that one of the flow variables
(typically discharge) must be specified, from which a value of A can be generated using
characteristic theory. If the flow is super-critical, then both characteristics enter the
99
4.6 Initial and Boundary Conditions
tu
p .1
t"-#-1 " L t,,+l
~o
flt
C+
t t"
Xo v,,, X, UI+ I12 X, X X X _
N 1 UN_Ill XN U N + l12 XN+l
(A) (B)
Figure 4.12: Illustration of the external boundaries at channel ends, (A) C_ character-
istic is used to solve the boundary condition at the upstream end and (B) C+ is used
at the downstream end.
region at the upstream boundary, and so values of A and Q must be specified. For the
downstream boundary, sub-critical fiow requires the specification of one variable (usu-
ally A) as in the upstream case. If the flow is super-critical then both characteristics
should propagate information from the upstream direction downstream.
We refer to the Figure 4.12, an intermediate point Q is added on the upstream
side and point M is added on the downstream side to make the solution progress
to point P and L, respectively. The values of A and Q are specified at the points
XQ, Xlj2, Xl+lj2, XN-lj2, XN+l/2 and XN+l at time t n and those at points P, Q, L, and
M are to be determined in order to progress the solution. Following are the steps
involved:
1. First Approximation:
(a) A first value of xQ and XM is found by using the slope of the characteristics
at points XQ, t n and XN+b t n
(4.95)
similarly
XM=XL- (~ +v'Bh)n !:"t (4.96)
N+l
100
4.6 Initial and Boundary Conditions
U
Q
= {u
U
1+1/2 -
(u 1+1/2 -
(u 1/2 - U 0 )
U)
1/2
X,/,-XQ
X1+'!s.,-xQ
x
'f
1 XI/2
'f
< xQ < XI+1/2,
1/2 - O.5l>x 1 Xo < XQ < XI/2
(4.97)
similarly state variable QM at the point M is given by
QM
= {u
U N+1
N+1/2 -
(u
(uN+1
N+I/2 - U N-I/2 ) XN+'i.,-XM
U N+1/2 ) XN+l-xM
x
'f
1 XN-I/2
'f
< XM < XN+I/2,
- - o.5Lix 1 XN+1/2 < XM < XN+1
(4.98)
(c) Equation (4.94), in the difference form, is used to obtain the state variables
at point P
. (4.101)
where
(4.102)
and
(4.103)
(d) In a similar way Equation (4.D4), in the difference form, is used to obtain
the state variables at point L
101
4.6 Initial and Boundary Conditions
(4.106)
where
(4.107)
and
(4.108)
This provides the required relationship that completes the information given by
the boundary condition. Provided h = h(t) or Q = Q(t), other state variables
can be estimated based on the Equations 4.101 and 4.106.
(4.109)
and
(b) Interpolation of state variables at point Q and M is carried out as per the
step (b) given above
102
4.6 Initial and Boundary Conditions
(c) Values of the state variables at point P are recalculated by averaging values
at points Q and the last approximation at P. This is given by:
(4.111)
(4.112)
Ap
Qp = TAl +BI (4.113)
where
and
103
4.6 Initial and Boundary Conditions
(4.116)
(4.117)
(4.118)
where
(4.119)
and
(4.120)
Steps (a) to (d) are repeated until a suitable accuracy is achieved (the accuracy is
defined here as an absolute difference in previous and current value less than equal to
lE-3). Thus, at the boundaries when one of the state variables is known then other
can be obtained using the relations 4.101, 4.106, 4.113 and 4.118.
104
4.6 Initial and Boundary Conditions
This is the free outfall boundary condition on the downstream side. It is an outlet
boundary which imposes no influence on the fluid in the domain where depths are
maintained and waves travel through without reflection. This type of boundary is
specified by fixing the boundary flux using the internal flow conditions. Thus boundary
flux is set equal to the internal cell-centered flux. Referring Equation 4.93:
F{f2 = J(Ui),
(4.121)
F'JV+1/2 = J(U'j,)
Based on Equations 4.24, 4.29 and Figure 4.11 the cell-centered flux is given by,
In the case of a solid wall boundary the velocity is set to be zero. Thus, for the solid
wall at x = a in the Figure 4.11, u(a, t) = 0 and the flux F at the boundary a and b
could be written as:
105
4.6 Initial and Boundary Conditions
As discussed above, the theory of characteristics is used to apply the stage hydrograph
boundary condition. One of the states, h in this case is known and system of equations
given in the Section 4.6.2 is solved to obtain Q.
In case of a discharge-stage hydrograph, both the state variables (A and Q) are specified.
Referring to Equations 4.24 and 4.29 the fluxes at the boundaries are estimated as,
n -
F1/2 - [ Q;peo
A(h spec )
Q.pec
+ 9 (h spec )2 (B"2 + (mt + m2 ) h,peo
6
) ] (4.124)
FJV+1/2 = [ Q;""
A(hspec )
Qspec
+ 9 (h spec )2 (B"'2 + (mt + m2 ) h,peo
6 ) ]
where Qspec and hspec are specified Q and h, respectively in the hydrograph. Typically,
this boundary condition is applied in case of a super-critical flow on an upstream
side (Garcia-Navarro and Saviron, 1992a, see).
Internal boundaries are defined at the channel junctions and structures. Interesting
literature has been found in solving these junction and structure boundary conditions.
Kutija (1995) presented an algorithm to solve a problem of computing flows in combined
dendritic and looped networks of channels. Underlying hydraulic properties were taken
into account to express the relations in simplest of graph-theoretical terms. This kind
of general solution for solving the network problem appeared for the first time and
106
4.6 Initial and Boundary Conditions
the algorithm is unconditionally stable, accurate and fast method of solution. This
algorithm has been used in NOAH (see Section 2.6.14) in association with implicit finite-
difference scheme. Garcia-Navarro and Saviron (1992b) presented the work related
to unsteady fiow simulation at a junction of open channels. The proposed solution
is valid for supercritical flows, hydraulic jumps and shock propagation through the
junction. Other solutions by different authors are: (Aral et al., 1998; Capart et al.,
1999; Hsu et al., 1998; Ji, 1998; Misra, 1998; Misra et a/., 1992; Naidu et al., 1997;
Noto and Tucciarelli, 2001; Reddy and Bhallamudi, 2004; Sen and Garg, 2002; Sobey,
2001; Swain and Chin, 1990). In CanaIFlow, the theory of characteristics discussed in
Section 4.6.2 is extended to solve the junction problem. This is also discussed in the
paper by Garcia-Navarro and Saviron (1992a).
An irrigation scheme consists of a network of canals organised in a hierarchical
fashion (Section 2.2). For example, a main canal takes off from the reservoir and acts
as a main distribution line. Distributaries or majors take off from the main canal and
in turn distributaries/majors have minors on them. The off-taking point (point at
which a child is connected to its parent) is through a structure or without a structure.
Figure (2.1, page 10) shows a typical schematic of an irrigation scheme. In order to
simulate flow conditions in the canal network, canal is divided in smaller sections called
channels. These channels are then assembled to simulate the flow in a network and
are connected to each other through nodes. Thus a node will have a channel and/or
multiple channels connected to it.
Boundary conditions are to be applied at the node junction where different channels
meet in order to solve the flow simulation problem for a network. The BC to be applied
at the end of channel connected to node depends on:
• Type of control structure present for example rectangular gate, circular or radial
gate;
107
4.6 Initial and Boundary Conditions
• Node
EB
rn Control Structure
C·6
C· Channel
1-4
llll
r[j C·7 J - Junction
EB • External Boundary
C·S
EB
Figure 4.13: Illustration of an irrigation network showing nodes, channels, control struc-
tures, junctions and external boundary points. Internal boundary conditions are solved
at the junctions. Shown are the junctions; J-l with three channels and a structure, J-2
as a three channel open junction, J-3 with two channels and a structure and J-5 as a
two channel open junction.
The following section discusses various cases considered by CanalFlow in solving the
internal boundary conditions.
Junctions and hydraulic structures are modelled as internal boundaries. A typical open
junction where split or confluence occurs is shown as J-2 in the Figure 4.13 and 4.14.
Based on equations (4.101 and 4.106) following could be written for the channels shown
in J·2 Figure (4.14):
Qi = AH Ai + BH
Qj = Alj Aj + B lj (4.125)
Qk = AikAk + Blk
where Qi, Qj and Qk are discharges and Ai, Aj and Ak are cross-sectional areas at the
interfaces of Channel i, j and k respectively. Al and BI for i, j and k are defined in
108
4.6 Initial and Boundary Conditions
Cl Cl
C ~
C, C C, C
(J-1) k (J-2) k
• C,
.~ • • C, • •
Cl Cl
(J-3) (J-5)
Figure 4.14: Junctions J-l, J-2, J-3 and J-5 as seen in the Figure 4.13.
Equations (4.102, 4.103, 4.107, and 4.108). For the given junction following equation
can be written to conserve the mass,
(4.126)
(4.127)
(4.128)
Ali [h(Bi +0.5mi h)] -A lj [h(Bj +0.5 mj h)] - Alk [h(Bk +0.5mk h)] = Blj + Blk -BJi
(4.129)
where Bi, Bj and Bk are bottom widths and mi, mj and mk are side slopes of Channel
i, j and k respectively. The equation is further simplified as
0.5 (AJi m; -A lj mj -Alk mk) h 2 + (AJi Bi - A lj Bj -Alk Bk) h- (Blj +Blk - Bli) = °
(4.130)
109
4.6 Initial and Boundary Conditions
This is a quadratic equation and could be easily solved to obtain h and in turn obtain
Q using Equation (4.12,,). In case of rectangular channels, the side slope equals zero
and h can be written as:
h= Blj + Blk - Bli (4.131)
Ali Bi - A lj Bj - Alk Bk
In order to achieve second order accuracy following steps are involved. Combining
Equations (4.113, 4.118, and 4.126) following could be written:
0.5 (Ali mi-Alj mj-Alk mk) h 2 + (Ali Bi-Alj Bj-A lk Bk) h-2 (Blj+Blk-Bli) = 0
(4.133)
If side slope is zero,
h= 2 (Blj + Blk - Bli)
(4.134)
Ali Bi - A lj Bj - Alk Bk
A three channel junction with a structure is shown as J-1 and J-4 in the Figure 4.13.
For example in J-1 (Figure 4.14), the discharge in channel Cj is calculated based on
the empirical discharge equation of the structure which is a function of water depth in
channel Ci and Cj. Recalling Equation 4.125 and substituting the value of Qj as,
(4.135)
(4.136)
In order to calculate the discharge through a structure, scenarios shown in the Figure
4.15 and 4.16 are considered. In case of a free flow condition, discharge through the
structure is function of a gate opening and upstream water depth hI, while in case of
submerged flow it is function of a gate opening, upstream water depth hI and down-
stream water depth h2. The procedure to obtain Q and h after Equation 4.136 is similar
to the one given in the Section 4.6.4.1.
110
4.6 Initial and Boundary Conditions
~ _____ L_ 1
ho h,
1 --------lmlI
Channell Channel 2
Figure 4.15: Channels connected through a structure and h2 < ho < hi (free flow)
The two channel junction without structure is shown as J-5 in the Figure 4.13. Referring
Figure 4.14 and assuming water depths and discharges to be equal at the connecting
interface we can write,
Qi = Qj,
(4.137)
hi=hj=h
and
(4.138)
Qj = A lj Aj + Blj
Substituting values of Ai and Aj and then simplifying we get,
(4.139)
this equation could be solved to get h, if side slopes are zero then h is given by
h= BI"-B
J I •"
(4.140)
AJi Bi - A lj Bj
To apply second order correction Equation 4.139 becomes,
(4.141)
this equation could be solved to get h, if side slopes are zero then h is given by
h= 2 {Blj - B Ji )
(4.142)
AJi Bi - Alj Bj
111
4.7 Concluding Remarks
Hydraulic structure
h. _____ L_ T
h,
h,
1 --------n 1
Channel I Channel 2
Figure 4.16: Channels connected through a structure and ha < h2 < hI (submerged
flow)
Junction of two channels with a structure in between is as shown in J-3 Figure 4.13. It
involves calculating discharge through a structure first as discussed in Section 4.6.4.2.
To solve the boundary problem at the structure, discharge on the upstream side of the
structure is equal to discharge on the downstream side. Based on the discharge another
state variable can be estimated using the characteristic equations for the channels. From
the Figure 4.14, we can write,
(4.143)
(4.144)
using above equation h is calculated based on the procedure given in the Section 4.6.4.3.
The key relationships and equations which have been set out in this chapter are the
basis of CanalFlow model. It has discussed in detail the governing equations, FVM
112
4.7 Concluding Remarks
formulation, exact Riemann solver, limiters and initial and boundary conditions. The
result is a robust I-D hydraulic simulation model for canal irrigation networks. The
model takes into account different hydraulic structures and junctions encountered in the
canal network. The exact Riemann solver described in the Section 4.4.2 is capable of
handling an advance of the wave front on a wet as well as dry beds. The model can also
handle the mixed-regime How and hydraulic jump in the channel. The methodology
to solve various boundary conditions at channel ends and at junctions described in
the Section 4.6.2 plays an important role in completion of this work. The challenging
part was to integrate the domain specific concepts described in this chapter with the
00 and relational database model described in the Chapter 5. CanalFlow can be
used to analyse hydraulic conditions in canals and test an operational plan for the
network. This has been demonstrated by testing CanalFlow using the benchmark test
suite described in the Chapter 6 with results of these tests presented in Chapter 7.
113
CHAPTER 5
5.1 Introduction
The first systematic approach to software development was called the traditional method-
ology or traditional paradigm; it sometimes is also called as the structured methodology
or structured paradigm. This methodology worked well with the small-scale systems
(with approximately less than 5,000 lines of code) but not with medium-scale (ap-
proximately 50,000 lines of code) or large-scale systems (approximately 500,000 lines
of code) (Booch, 1994). The two basic blocks of software systems are the operations
performed by the system and the data on which the operations are performed. The
traditional paradigm essentially ignores the data in favour of the operations. In con-
trast, the object-oriented paradigm pays equal attention to operations and data. Thus,
object-oriented paradigm proved to be the solution and a better approach for developing
software systems. In 1980s more than 50 different object-oriented methodologies were
published (Schach, 2004). Three of the most successful methodologies were Booch's
method, Jacobson's Objectory, and Rumbaugh's Object Modelling Technique (OMT).
This chapter discusses object-oriented design and UML in brief and then moves
on to discuss the design of the CanalFlow model in detail. The software design is
114
5.2 Object-oriented Design and UML
discussed with help of class diagrams. Last section is dedicated to the database design
and discussion on the data model used in CanaIFlow.
Object-oriented analysis and design (OOAD) consists of two phases; first is analysis
which is necessarily finding and describing the objects or concepts in the problem do-
main. For example in a canal hydraulic model, some of the concepts include channel,
network, structures, and the numerical solver. The second phase is a design phase
which focuses on defining software objects and the way they collaborate to fulfil! the
requirements. For example, a channel software object may have a start node, end
node, length attribute and methods like getFroudeNumber, getWaterDepth, apply-
BoundaryCondition. Finally, these designed objects are implemented using an 00
programming language, in this case CanalFlow implementation is done using Java.
Unified Modelling Language (UML) is a visual language for specifying, constructing and
documenting the artifacts of systems (OMG, 2003). It is an universally accepted no-
tation for the requirements, analysis, and design of a software system (Schach, 2004).
In early 1990s, the most popular object-oriented analysis and design methodologies
were OMT developed by James Rumbaugh at General Electric Research and Develop-
ment Centre (Rumbaugh et al., 1991) and Grady Booch's method developed at Ratio-
nal (Booch, 1994). All object-oriented analysis techniques essentially are equivalent, as
are all object-oriented design techniques, so the differences between OMT and Booch's
method are small. In 1994, Rumbaugh joined Booch to develop a methodology that
would combine OMT and Booch's method. Thus, UML, a common notation system
for representing an object-oriented software system was developed. UML diagrams en-
able information technology professionals to communicate their ideas more quickly and
more accurately than if only verbal descriptions were used. This chapter will take help
of the class diagrams to explain the underlying architecture of the model.
115
5.3 Software Design and Development
This research followed iterative development model which involves early programming
and testing of a partial system, in repeating cycles. We relied on short quick develop-
ment steps, feedback, and adaptation to clarify the requirements and design. In this
approach, development is organised into a series of short, fixed-length cycles called iter-
ations; the outcome of each is a tested, integrated, and executable partial system. The
system grows incrementally over time, iteration by iteration, and thus this approach is
also knows as incremental development. The requirements phase was short and brief
requirements were analysed and captured in the First year report (Shende et al., 2003).
Basic classes were designed and implemented with development of a core solver class.
Slowly the code evolved with a test-driven and continuous integration approach.
The classes in CanalFlow model are organised in packages based on their functionality.
Packages are containers for classes and are used to keep the class name space compart-
mentalised (Naughton, 1996). These packages are stored in a hierarchical manner and
are explicitly imported into new class definitions. Thus, they are both a naming mech-
anism and a visibility restriction mechanism. A package consists of class(es) and/or
interface(s). A class encapsulates all concepts and defines shape and behaviour of an
object, while interface is designed to support dynamic method resolution at runtime.
An interface is jnst like a class, but it does not define instance variables and come with
methods declared without any body. Per se they do not provide any implementation
details but the classes do. Thus to implement an interface, class needs to provide
implementation of the complete set of methods from the interface and these methods
must match the type signatures of the interface method parameters exactly.
In CanalFlow, the physical concepts described using objects could be explained as,
• The canal irrigation network consists mainly of canals, reservoirs, junctions and
control structures;
116
5.3 Software Design and Development
• The canal may be divided into sections called as channels which are connected
to each other through nodes. These channels might have a controIling structure
at the start and end caIIed as node structures and these channels take some
geometrical shape;
• In order to solve a flow problem for a network, the network components are to be
assembled together with some initial data. Each of these components have got
data and methods which act on this data;
With this background, we can now look at the canalfiow package which has got sub-
packages and classes defined in it representing the concepts defined above. The overview
of canalfiow package is shown in the Figure 5.1.
117
5.3 Software Design and Development
=_1""",,1'''
=
~~ ~D
E!J ~
CIo "bUCIUtII _
I----J
., B B
= =
IM'-'-I ,I '·-·~·I
I<GUI~~'II~~Su"d.J
,,
,,
R...... ""IDSo1_
a
OnoD8_
~-.
I
I
D
1..'
I I
MIIlnF_
"--
I -
~
~ .......
, .. :,-----
-- I -
t.' 1•.'
I
... 1.."
•..
,
, I
.,
.
.', " ..
' ,
. -"-.
-,-
lS
, , ,
.
.,
~
,
Clrcul.
,
,
. ,
,, ••
, ,
, ,
, , ,
,
w••
Roc:III"II_
-. I
- •
Figure 5.1: Components of canalfiow package. Shown axe various classes and their
relations with other classes. For example in a network package, Rlemann1DSolver is a
sub-class of a Channel class additionally it implements OneDSolver interface. Also, Net-
work consists of many Channel, Node and Reservoir objects while a Node might/might
not have a Junction.
118
5.3 Software Design and Development
5. canalflow. network: This is the most important package of the model with imple-
mentation of classes like Channel, Network, Node and Riemann1DSolver.
7. canalflow.utilities: For different utilities which in turn are used by all other classes.
5.3.2 Tools
Various software tools used for the design, development, testing and documentation of
CanalFlow are listed in the Table 5.1. The Java code was implemented and compiled us-
ing Java 2 Platform Standard Edition (J2SE) 5.0 also called as JDK 1.5.0 or Tiger. J2SE
5.0 is a significant release including many new features and updates while preserving
compatibility and stability of previous versions (http://java.sun.com/j2se/1.5.0/index.jsp ).
Table 5.1: Different software tools used during various phases of the research.
119
5.4 CanalFlow Packages
5.4.1 canalfiow.channelstructures
5.4.2 canalfiow.database
120
5.4 CanalFlow Packages
121
5.4 CanalFlow Packages
~~~~JtJj~tI r,:----i
I
~~1I·;1 I
n.-.ct.lStNctw. [~-lirkn
If-,-..- ..-~-.....
- ..------~.-- ~: ;:'-:)/]1
,::;:=
, ; G !Ildx: daIbe()
#Gi!h:f>;:doo.t>Iefl'
--~~~ih
~
# .l fttwni.~l1~':dD.tlIe
;.lI frtJtrJ¥'tSlo:leSloPe:tb.tIItI
r-~or=n !.J
# ;, IstISb!oSIDpB: dc:UJIe[) ~
# ;,.~:doltie[f
;rI~:~
; rI tmrfIcUOI'IWIIh. dol.tIIII
;~sldeSlopeo:~
#' rlloBaIorrWdIh: doo.tie
122
5.4 CanalFlow Packages
Figure 5.4: DBConnection class diagram. This class is used by other classes for data-
base related functionalities.
Figure 5.5: Illustration of a Canal class diagram. Canal has one or more Channel
objects and is a part of a network package.
123
5.4 CanalFlow Packages
5.4.3 canalflow.gui
5.4.4 canalflow.network
This is the most important package in CanalFlow. It implements classes like Channel,
Junction, Network, Node and Solvers package. In the current version, only one 1-D
solver is implemented, i.e., Riemann1DSolver. In future, solvers could easily be added
to this sub-package and user can specify the option.
124
5.4 CanalFlow Packages
Class diagrams of Canal and Channel are shown in Figure 5.5 and 5.10 respectively.
The irrigation network consists of network of canals. These canals are built using
Canal class. For computational purposes, these canals are broken down into number of
channels (or reaches). These channels have x, y, and z co-ordinates and other hydraulic
details. Each channel has a parent canal object.
Channel class has all the necessary details required for the computations. Each
channel class has got its own solver. Basically, the solver class (implemented using
OneDSolver ihterface shown in Figure 5.6) takes care of solving the governing equa-
tions. In this case, an exact Riemann solver implemented using Riemann1DSolver class
(Figure 5.7), which extends Channel class so that it can use all the public and protected
variables and methods implemented in the Channel class.
Network class handles the tasks of building the network components, initialising
actual network, initiating and controlling the hydraulic simulation activity, initialising
all network-wide parameters and initiating output task. In short, this is a sort of a
controlling class which keeps track of time and decides when to output and stop the
simulation process.
Node class plays vital role in bnilding a network and channels. Each node has its
x, y, and z co-ordinates defined. Each channel has got a fromNode and toNode. A
node might have a reservoir. Similarly, few nodes in the network would be junction
nodes where multiple channels meet.
Node class also handles the process of applying boundary conditions. After every
time step boundary condition method is run on each node. The node then decides which
boundary is to be applied and then pass on the instructions to respective components.
For example in case of a wall structure, node passes control to the wall class. Wall class
then checks what channel it is connected to and in turn asks corresponding channel to
apply wall boundary. Similarly, if multiple channels meet then the control is passed
on to Junction class. Junction class then calculates the discharge and depth at the
junction and passes control to respective channel objects to apply depth-discharge
boundary. This class also does some housekeeping in terms of scanning the network,
analysing the nodes and sorting channels in different categories like channels taking
125
5.4 CanalFlow Packages
canalflow.network.soIvers':
Riemann1DSoluer
Figure 5.6: OneDSolver interface. The new solvers developed in future could implement
this interface.
~rI'_os-o'''''''
,".ao~,,,,,,,
~"'&OIOIIOfV..-.o'_
-.".orcT...a:_
.
.. " Ohop(): ......
,,~:
"'iI~:"""
... illIIUr():_
126
5.4 CanalFlow Packages
_;\10<*"'1
_:-
_'SWcI
.....,:W
....oe: _
Figure 5.8: Illustration of a Network class diagram. This is a central class which initiate
all processes pertaining to building network, simulation process and time keeping op-
erations. The concept of multithreading is used to run the simulation and visualisation
processes simultaneously as two separate threads.
127
5.4 CanalFlow Packages
Figure 5.9: Shown is the Node class diagram. Node objects are building blocks of a
network. Important role of a Node class is to initiate right kind of boundary condition
based on the available data but the actual application of boundary condition is done
by the Channel class.
128
5.4 CanalFlow Packages
off from the node, channels coming to the node, channel with structure and channel
without a structure.
Junction class specifically implements the boundary condition methods for a junc-
tion. A node initialises the junction and shares the information with it such as number
of channels connected to it and connection end of each channel (from or to channel).
Depending on the number of channels connected to it and presence of structures on the
channel, junction object decides which boundary condition method to execute.
Reservoir is one of the major network components and maintains the data and var-
ious states pertaining to reservoir. Reservoir object is initialised with some water level,
dimensions and a reservoir rule curve. Rule curve is a look-up table where the relation
between depth-discharge-submergence area is defined. Each reservoir is connected to
a channel and flow is controlled through an intermediate hydraulic structure. This
structure governs the flow in the channel. A schedule can be applied for opening and
closing of the gates.
5.4.5 canalflow.nodestructures
129
5.4 CanalFlow Packages
-,-
-.r.:H
_,lot
~:W
_:1'11
,,,,,1'11
q:<Io<.t:IIeOI'
<!SI<Irt:·dcI;lIIetJl
~:~
~:~
.'CIOIbIoIIl
~,~
_:-
~,~
"',"""*""Il
- ......
111:li0<*'ii
,._-,,'
_:1b,.tJIe[JJ]
""'Type'rc
_:CcmodIDn
..._:SUiI_
1cCo!tI>:_
-,-
~:
1cS!rD:"
130
5.4 CanalFlow Packages
~
" r8 gelFroucls(); doIJJill
... ~ getHYl'hlLllleb~:~, '. #" rI ~.. ::do,~'
, 4' ill i'IIhcte~:r.t', '
" rI ~1ltAeDer.thO: ~ ; #" rJi hScu-ca: doUIIII
, ... r1!~~rOl!Wes();doIA:IIe i #" rJl ~; doU::Ib
i ... rJ! gelHydral.llcRdJ$(); double i ., rJJ m1lc~,irI! '
... rI gelHy(i'lIUAeRadlus()) doU'JIe 4' rJl~:Vectot
. ... rI get(l(): W #" rJi acuce}dooJ:lIe{II-
" ri! getlndeXO: 'lit #" rJIlo1b;fe; No:w:t&,
.,... rI geiLeltSlCope(): ~ #" rJI W:lodeStrUdlte :'Node~
... rI QetM.r.tllls(i: doleie #" rI vellx:1),: doubien'
.... ;1 ~:dolJJIe #" rJI weterDePh; drJlde(J .
... rI getNew~:doI.tIIt 4' rJl xPol't-;'dotbI$'
... .J!o gelNewflOneO:-dOuble
: ... rI gdNew)(Po;-to; doI.tIIt
... r1! ge4RirjtSIopeO:~
... rH getSId!ISI~; double
... rI getSkmSlOP8(); doubla•
... r1! getSpedflcf'oroe()'~ ~
... rI getSlIJteVerW:Jle():cb.t'oIII,
" rJl geIS1rucb.irelmt!lxO: 1nl
... rI geii"6No!lnsiructure~o~ri
, ... rI getT~~: dooJlIe,
... r:I getVoUne(): dt.:U:lIa
~d!;ge\WIII:~rO\!pIt():dOl.ti:l
131
5.4 CanalFlow Packages
Figure 5.11: Junction class diagram. It takes care of solving boundary condition prob-
lem at the junctions.
l.Ii.':iltiim
.....
132
5.4 CanalFlow Packages
133
5.4 CanalFlow Packages
.;
.; " SETTING TYpe GATE; t't
1Fj~9--11 .; f:boIa1dIwyCQnditlon:Jrt
U tt ~:cIoUIt
,. , em: Chanr)eI
,. 'llllml,: Stmg
, ,. "nooIGtIII_~Irii'
.; 9 ~: boOIetri
.; , $eltingType,:,",
" .; 't, $\de: booIeM.
,. ',stroea~::!:lI""Ii-
.; ,'lIIfI):W
134
5.4 CanalFlow Packages
RectangularGate
~~~ ;p:;~~::;iq=;q'it;=~jFi~jF---ii:fGlfl
r:~;cq"..cfi.. ,11- =RWItSH,;1 SQlE.ception I"st,t. . .1if d
135
5.5 Database Design
• Table: Table is a collection of data in a defined structure. In the data model, the
table is displayed in a window-like manner. The table's name is displayed in the
title bar, the table columns are printed below and are indicated by an icon. A
key icon means this column is in the table's primary key list.
• Primary Key: Usually one or more columns are defined as the table's Primary
Key (PK). These columns must not contain two or more data values which are
equal. That makes it possible to clearly identify each record in the table by
the Primary Key. For example channel ID will not be duplicated in the channel
master table.
• Indices: To make the database find a specific record in the database more quickly,
it is possible to define an index on one or more columns. Indices are also used to
improve speed when two or more tables are joined together.
• Relations: Relations can only be placed between two tables. They define the rela-
tionship between the tables and can create a Foreign Key (FK) reference. Tables
can be connected by a one-to-one (channel and a channel section-type relation),
one-to-many (network-channel relation) or many-to-many (channel-output at dif-
ferent time) relation.
The Figure .5.16 shows the complete data model. Description of all tables is pro-
vided in Appendix A. The data model illustrates different tables in the database and
relationship between different tables. The whole data structure is organised in a hi-
erarchical fashion. The network is initialised first and then network components like
136
5.5 Database Design
""''''';'7,
"
, ,."
1:~:Eil~i;
I~~~·· .;,1'
'· .
I• OV',;, ""'" ,.,
. I1U1' i
I
GIII.P'Il: l1NYINT
gleJ>ilt: oo..ae(~,2)
,
919.:.....d1:t:IOLIlE(4.2)
.'
I::;j.~=<;:;)
I:~;;,;= K'~;"
,.
, , . ,
Figure 5.16: CanalFlow data modeL Shown are all tables, their attributes, data types,
relations and primary keys. Table names prefixed with 'm' are master and 't' are
transaction tables. Master tables are in which data transactions (inserts and deletes)
are limited while in transaction table they happen frequently.
137
5.6 Concluding Remarks
canals, channels, hydraulic structures and channel structures. Each of the network
component is initialised with some initial and necessary master data. As explained
before the relation between tables could be one-to-one, one-to-many or many-to-many.
For example, a network can have many canals but a canal can have only one network
associated with it. Similar relations can be derived for canals and channels. Relation-
ship between a channel and a channel structure is of type one-to-one because a channel
will have only one channel structure and a particular channel structure will have only
one channel associated with it. The Appendix A provide details of all tables in terms
of their attributes, data types, primary and foreign key, whether or not null values are
allowed, default values and auto increment keys.
In this chapter, development of the CanalFlow model has been described based on the
00 concepts, Java and MySQL described in the Chapter 3 and hydraulic modelling
theory presented in the Chapter 4. The important thing to mention regarding the
CanalFlow is the design has kept separate the domain specific logic from the data model
and GUI which is extremely important in maintaining and doing enhancements in the
software. The design implements simple 00 model to define various components in the
network. These components then interact with each other to perform tasks like building
the network, running simulation and writing output to the database. Java's multi-
tasking feature has been used to run the simulation and update data simultaneously
on the graphical interface. The data pertaining to various network components, initial
and boundary conditions and system-wide simulation parameters can be defined in a
data model presented in the Section 5.5. The network building process builds Java
objects by retrieving necessary data from the tables. The procedure to enter the data,
build network and simulate flow conditions in a network is given in the Chapter 7.
138
CHAPTER 6
6.1 Introduction
The recent years have witnessed a proliferation of software packages used in hydraulic
research and engineering applications. The advent of more capable computing platforms
has resulted in an extensive use of mathematical models in all aspects of engineering,
including hydraulics. A great deal of knowledge and skill is being transferred in this
way, from experts to end users. While it is exciting to see the products applied to
important technological issues in the society, this also places increasing demands on
the technical quality of the software. The utmost important thing for the numerical
modeller and an end user is the concept of accuracy and adaptability.
This stresses the need to examine and support the validity of the models and of
the results they produce. Testing of the model can be done by direct comparison
of numerical results with the analytical, laboratory or field data. In order to test
the unsteady flow canal model on a real canal system, the real canal system should
be complex enough to provide a good test, yet not too difficult to model. Physical
dimensions and hydraulic properties must be well defined. References in literature
regarding such systems are available. Rogers et al. (1993) have described the Central
Arizona Project and the Salt River Project for this purpose. Further, Parrish In and
Burt (1993) have described the CalPoly model for the purpose of generating data for a
139
6.2 Need for the Testing
series of predetermined flow scenarios. Attempts were made to acquire data from these
projects (Personal Communication - Burt, 2004; Personal Communication - Parrish,
2004; Personal Communication - Schuurmans, 2004) for this research work. These
attempts failed because not much research has been done in this area in recent years
and the authors seem to have lost the old data.
This chapter explains, in brief, the need for testing and shortlisted benchmark test
cases adapted for this research work. The benchmark test cases are classified under
three broad categories and detailed description of the test configuration is provided in
the Section 6.4. The results of these tests are discussed in depth in the Chapter 7.
Development of a reliable numerical model for unsteady flow in channels and channel
network is a very formidable task because of its complex formulation, non-linear equa-
tions, bi-directional and rapid response, and definition of suitable boundary conditions.
The starting point of all models is the same, namely, governing equations and the end
point of all models is the same, the predicted response patterns. But the route from
start to end may follow a wide variety of paths (Sobey, 2001). Although this may be
satisfactory, there is a need to demonstrate model credibility. The purpose of testing
a computer-program is to verify the developed computer program. Testing involves
operation of a system or application under controlled conditions and evaluating the
results.
Sequence of benchmark tests can be used to spotlight the physically and numerically
significant response patterns that are expected to be within the predictive capabilities
of steady/unsteady flow models. The purpose of these benchmark test cases is to
provide a potential user with a comprehensive picture of the capabilities and potential
blemishes of a particular code (Sobey, 2001).
140
6.3 Testing Methodology
• Capability: The capability can be assessed objectively by testing the most com-
monly required features of a software package. CanaIFlow is assessed using the
looped network, weirs and a hypothetical irrigation network test case (see Section
6.4.5 and 6.4.6).
• Reproducibility: Under this category model results are compared with the experi-
mental or real world dataset to assess the reproducibility. In this case, CanalFlow
results are compared with the analytical and laboratory dam-break data (see Sec-
tion 6.4.8).
Apart from the DEFRA/EA test cases, one more test was conducted for the sake of
completeness of this research work and to test the applicability of CanalFlow for canal
irrigation networks. The data for this test was adapted from Misra et al. (1992) and
141
6.3 Testing Methodology
the test is categorised under the Capability group (see Section 6.4.7). The darn-break
in a straight rectangular channel was also added for the sake of completeness for which
the data was adapted from Gottardi and Venutelli (2004). Table 6.1 presents the list
of test cases shortlisted for this research. Each of this test has a limited objective,
seeking to focus on crucial problems in relative isolation. They have been grouped into
numerical accuracy with identifier N, capability with identifier C, and reproducibility
with identifier R. All of the test cases carried out for this research work with identifier
Nand R are single-channel problems while test cases with identifier C are for multiple
channels or for network of channels.
Table 6.1: The benchmark test suite used to test CanalFlow. Listed are the eight test
cases classified under three broad categories.
142
6.4 Test Cases
This section provides specification of the tests listed in the Table 6.1. The test cases
discussed here tells about the test objectives, configuration, dataset and additional
notes. The details of model building process for these tests and results are provided in
the Chapter 7.
6.4.1.1 Objectives
• Present the particulars for developing and undertaking the test with CanalFlow
and discuss the associated results.
The test is undertaken in six sub-parts as defined in the Table 6.2 (Crowder et al.,
2004b,i). For each part of the test 101 cross-sections of rectangular shape are defined
at 1.0 m spacing as illustrated in the Figure 6.1. The elevation of each cross section is
defined such that the required flow regime is produced. Additional cross sections, i. e.,
interpolated are not added to the test configuration. A constant Manning's n value
of 0.03 is used throughout the channel for all parts of the test. Flow and water level
boundary conditions are defined as given in Table 6.3 below.
143
6.4 Test Cases
Table 6.2: Definition and type of boundary condition used for the sub-parts N.1.l to
N.1.6 in the sub-critical, super-critical and transitional flows test case.
CanalFlow is tested under the quasi-steady (Subsection 6.4.1.4) flow boundary con-
ditions for sub-parts N.1.1 to N.1.5. The boundary conditions, as defined for each part
of the test, are used at time t = 0 and extended through to 01:00 hour respectively.
Part N.1.6 of the test uses the same test configuration as that used for Part N.1.5,
with the exception of the boundary conditions. For Part N.1.6 of the test, the upstream
boundary is fixed at 20.0 m 3 / sec for a period of 24:00 hours. The downstream boundary
is set at 0.602 m for the first 06:00 hours and then it is increased linearly from 0.602
m to 2.0 m between 06:00 hours and 18:00 hours, after which it remains at this level
for a further 06:00 hours.
144
6.4 Test Cases
Table 6.3: Boundary condition details for sub-parts N.1.1 to N.1.5 in the sub-critical,
super-critical and transitional flows test case.
I
Downstream Upstream
Boundary Cross-Sections Boundary
0·101
6.4.1.3 Dataset
The dataset entered in the CanalFlow database to execute this test case is as given in
the Tables B.l to B.5 in Appendix B.
6.4.1.4 Notes
ii. The original test specification (Crowder et al., 2004b) defines that test part N.1.1
145
6.4 Test Cases
to N.1.5 be simulated using a steady and quasi-steady state but CanalFlow soft-
ware does not have a separate steady state solver, thus, the software is run using
an unsteady state solver and results are compared with quasi-steady solver of
other softwares (ISIS, MIKE 11 and HEC-RAS);
iii. For parts N.1.2, N.1.3, and N.1.5. discharge-depth hydrograph is used as an up-
stream boundary condition; and
iv. The results of this test case are discussed in the Section 7.4.2.
This test is an assessment of the capability of the software package to calculate the
normal sub-critical flow depth and the normal super-critical flow depth in a triangular
channel (Crowder et al., 2004d).
6.4.2.1 Objectives
• Assess the ability of CanalFlow to calculate the normal sub-critical flow depth
(Part N.2.1) and the normal supercritical flow depth (Part N.2.2) in a triangular
channel under quasi-steady (Subsection 6.4.2.4) boundary conditions; and
• Present the particulars for developing and undertaking the test with CanalFlow
and discuss the associated results.
The test configuration is illustrated schematically in Figure 0.2. For Part N.2.1, the
triangular channel is defined by eleven cross-sections that are placed 300 m apart each
with a side slope of 1:2, a constant Manning's n value of 0.035 and a constant bed slope
of 0.001.
146
6.4 Test Cases
I
Downstream Upstream
Boundary Cross~Sections Boundary
Chainage = 0.0 m 0-11 Chainage = 3000.0 m (N.2.1)
Slope = 0.001 (N.2.1) Chainage = 150.0 m (N.2.2)
Slope = 0.02 (N.2.2)
Figure 6.2: Schematic illustration of the triangular channel test case (N .2). The channel
of different length and slope but same number of cross-sections is used to conduct two
sub-parts of the test. The flow direction is from right to left (higher elevation to lower
elevation)
For Part N .2.2, the triangular channel is defined by eleven cross-sections that are
placed 15 m apart each with a side slope of 1:2, a constant Manning's n value of 0.035
and a constant bed slope of 0.02. Interpolated cross sections are not used in the test.
For Parts N.2.1 and N.2.2 of the test, CanalFlow is tested under the quasi-state
boundary conditions as defined in the Table 6.4. The boundary conditions defined for
each part of the test, are used at time 00:00 hour and extended through to 01:00 hour.
Table 6.4: Boundary condition data for parts N.2.1 and N.2.2 of the triangular channel
test case.
6.4.2.3 Dataset
The dataset entered in the CanalFlow database is as given in the Table B.G in Appendix
B.
147
6.4 Test Cases
6.4.2.4 Notes
ii. CanalFlow is run in a full unsteady mode and for the super-critical flow test
case (Part N.2.2), discharge-depth hydrograph is used as an upstream boundary
condition; and
iii. The results of this test case are discussed in the Section 7.5.2.
These assumptions are acceptable when the ratio of wave amplitude to the mean water
depth is relatively small. The analytical solution calculates the water elevation in a
channel of finite length with a closed end by superimposing the incident wave with the
reflecting wave, as schematically illustrated in Figure 6.3. The test set-up is adapted
from Crowder et al. (2004f).
6.4.3.1 Objectives
148
6.4 Test Cases
HWL-21.0m
~ ---:...:::..::.- ---- - - - - - - -- - -- - -- - - - ~WL=20.0m
25.0 m
________________________ ~ ~WL= 19.0m
1-----.. Reflecting wave
"~---------------------~)I
Closed end 100 Km Open end
Figure 6.3: Schematic diagram of a tidal wave in the idealised channel. The upstream
side of the channel is closed while depth is varied on the downstream side, thus, creating
an incident wave which then is reflected after hitting the closed end of the channel.
• Present the particulars for developing and undertaking the test with CanalFlow
and discuss the associated results.
The model sets up a rectangular channel of 100 km length and uniform elevation (no
bed slope) using 201 cross-sections of bed width of 1000 m and a side wall height of
25 m at an equal spacing of 500 m. Constant Manning's n roughness value of 0.025 is
used along the complete length of the channel.
Cross-sections are defined with a physical geometry and not interpolated cross-
sections. A flow/time boundary of zero inflow for the duration of the simulation is
defined at the upstream side of the channel so as to represent the closed boundary. In
CanalFlow this is achieved by specifying a waIl structure at the upstream end. The
open boundary at the downstream end of the channel is specified with a head/time
boundary (stage hydrograph), simulating a sinusoidal wave profile of time period 2.0
hours and amplitude 1.0 m.
149
6.4 Test Cases
6.4.3.3 Dataset
The initial data entered in the CanalFlow database is as given in the Table B.7 in
Appendix B.
6.4.3.4 Notes
ii. The results of this test case are discussed in the Section 7.6.2.
This test is an assessment of the ability of the software package to recreate the special
case of an unsteady flow, known as the monoclinal rising wave, as a typical case of
uniformly progressive flow (Crowder et al., 2004g).
The Monoclinal rising wave described by Chow (1973), is an example of uniformly
progressive flow. The wave is regarded as a special case of unsteady flow as it describes
the transition from one set of steady state condition to another. The theoretical analysis
of the wave is such that for a given channel of uniform cross-section and slope, the
flow at the upstream end increases uniformly over a given length of time, from an
initial constant value to a subsequent constant value. The resultant wave is one which
possesses a stable longitudinal profile, of the form shown in Figure 6.4, the shape of
which remains unaltered as the wave travels through the channel.
6.4.4.1 Objectives
• Determine whether CanalFlow is able to recreate the special case of unsteady flow
(monoc1inal rising wave) described by Chow (1973) as a typical case of uniformly
progressive flow; and
150
6.4 Test Cases
~ ~,•••••• fo-_U-"w_t-.
•••••••••••••••••• •••
Downstream stage
Figure 6.4: The Monoclinal rising wave profile. It is a translatory wave of stable form,
i.e., a uniformly progressive wave, that travels down the channel at constant velocity U w
from an upstream region of uniform flow of depth hI and velocity Ul to a downstream
region of uniform flow of depth h2 and velocity U2. During a time interval f!.t, the
wavefront moves forward a distance equal to Uwt.
• Present the particulars for developing and undertaking the test using CanalFlow
and present the associated results.
151
6.4 Test Cases
Upstream Downstream
Boundary Cross-Sections Boundary
0-201
Figure 6.5: Schematic illustration of MonocJinal wave test case (N.4). A rectangular
channel of 1000 km is defined to perform the test with a discharge hydrograph and a
stage hydrograph specified at the upstream and downstream boundary, respectively.
The flow direction if from left to right.
6.4.4.3 Dataset
Initial conditions for stage and flow at every cross section are entered into the model
based upon Chow's analytical solution and adapted from the data sheet used in DE-
FRA/EA benchmark test suite (Personal Communication - Crowder, 2005). The data
entered in CanalFlow is specified in the Table n.8 in Appendix n.
6.4.4.4 Notes
i. Apart from the above-mentioned test configuration this test case is run using
different limiters, namely, minmod, monotonised centered, superbee and van Leer;
and
ii. The results of this test case are discussed in the Section 7.7.2.
152
6.4 Test Cases
6.4.5.1 Objectives
• Assess the ability of CanalFlow to calculate a split and confluence flow scenario
(a looped system);
• Assess the solution criteria and the stability of the software package under quasi-
steady (see Subsection 6.4.5.4) boundary conditions; and
• Present the particulars for developing and undertaking the test using CanalFlow
and present the associated results.
The test configuration is illustrated schematically in the Figure 6.u with a channel
geometry details as given in the Table 6.5. A constant Manning's n value of 0.012,
0.0125, 0.013 and 0.0135 is used for Channels A through to D respectively. Channel
B is three times as long as Channel C with the same vertical drop between the two
junctions.
Figure 6.6: Schematic illustration of the looped system (C.1) test. There are four
channels: Channel A, Channel B (B1 and B2 combined), Channel C and Channel D.
At the downstream end of Channel A the system diverges (splits) into Channels B and
C, which then converge to Channel D at their downstream sections. In CanalFlow test
setup, Channel B is divided into 2 sub-channels of equal length.
Interpolated cross sections spaced at 5.0 m distance are used for this test. The test
configuration can be summarised in the Table 6.5. The software package is tested with
153
6.4 Test Cases
Table 6.5: Channel geometry details for the looped system (C. 1) test case.
6.4.5.3 Dataset
The initial dataset entered in the CanaIFlow database is as given in the Table B.9 in
AppendixB.
6.4.5.4 Notes
iii. In CanalFlow, the test is carried out using full unsteady flow simulation with 5.0
m as the distance step for all channels (interpolated cross-sections); and
iv. The results of this test case are discussed in the Section 7.8.2.
154
6.4 Test Cases
6.4.6.1 Objectives
• Assess the ability of CanalFlow to model a Broad Crested weir under steady and
unsteady boundary conditions;
• Assess water level and head loss results for the weir under both free flow and
drowned flow conditions; and
• Present the particulars for developing and undertaking the test using CanalFlow
and present the associated results.
The test configuration is illustrated schematically in the Figure 6.7 and the parameters
defining the Broad Crested weir are given in the Table H.6.
r - Weir
r - Channel A EfiI._____-'/:. ___ C_h_ann_e_I_B_-I,
u...
p'-"e-.-m---''---------I
Boundary
m
---+- 2.0 m+--
Downstream
Boundary
100.0 m tOO.Om
I~ )1
Figure 6.7: Schematic illustration of the weir (C.2) test case configuration. A weir is
placed between Channel A and Channel B, both of which have two cross sections at a
spacing of 100 m, a constant Manning's n of 0.014 and a bed slope of 0.005 and 0.001,
respectively. The flow is from left to right.
155
6.4 Test Cases
Parameter Value
The software is also tested with an unsteady boundary condition. The upstream bound-
ary is fixed at 0.15 m 3 / sec for a period of 24:00 hours. The downstream boundary is
set at 0.3 m for the first 06:00 hours and then it is linearly increased from 0.3 m to 0.8
m between 06:00 hours and 18:00 hours, after which it remains at this level for further
06:00 hours.
6.4.6.3 Dataset
The dataset entered in the CanalFlow database is as given in the section above and in
the Table 6.6.
6.4.6.4 Notes
156
6.4 Test Cases
ii. In CanalFlow, the test is carried out using full unsteady flow simulation with 2.5
m as the distance step (~x) used for both channels; and
iii. The results of this test case are discussed in the Section 7.9.2.
6.4.7.1 Objectives
• Assess the ability of CanalFlow to handle the structures and ability to operate
these structures as per the gate operation schedule;
• Compare discharge and depth hydrographs at the upstream and downstream side
of the channels and carry out mass balance at selected junctions; and
• Present the particulars for developing and undertaking the test using CanalFlow
and present the associated results.
The test configuration is as shown in the Figure 6.8. The network consists of 41
channels, 9 structures, 20 junctions and a reservoir. Total length of the network is 43.5
km. Details of the channels are as given in the Table B.1O in Appendix B while details
of the control structures are presented in the Table 6.7. Constant depth of 5.0m is
malntained in the reservoir and structures are operated according to the gate schedule
given in the Table 6.8. Total simulation time is 6 hours.
157
6.4 Test Cases
oControl structure
~ Transmissive boundary
C8
CI8 CJ7
Cl' C7 C35
CI5
S-9
CI1
C33
J2
S-I
Cl C4
025
Figure 6.8: Schematic illustration of the Irrigation network test (C.3) case. The chan-
nels are numbered with prefix C, control structures with S and junctions with J. Rec-
tangular gates are used as control structures and they are operated as per the given
schedule. Constant depth of water is maintained in the reservoir and water flows
through main channels Cl to C4 and then through other channels connected to one of
these channels.
158
6.4 Test Cases
Table 6.7: Details of the rectangular gates used in the irrigation network (0.3) test
case. Cd is the coefficient of discharge used for the free and submerged flow condition.
Table 6.8: Gate operation schedule for the rectangular gates used in the irrigation
network test case (Refer Figure 6.8). The gate number 1, 2, 5 and 8 are opened in
first hour, kept at the same opening during second hour and then closed in third hour.
Between 3 and 6 hours they are kept at the same gate opening. Gates 3, 4, 6 and 7 are
kept at constant gate opening and are not operated at all.
159
6.4 Test Cases
6.4.7.3 Dataset
The dataset is entered in the CanalFlow database as given in the Table B.lD in Ap-
pendix B and in Table 6.7.
6.4.7.4 Notes
i. The test is carried out using full unsteady flow simulation and LJ.x = 2.5 m for
all channels; and
ii. The results of this test case are discussed in the Section 7.10.2.
This is an assessment of the ability of the software package to replicate the behaviour of
a surge wave, caused by the sudden collapse of a large body of water, in a rectangular
channel (Gottardi and Venutelli, 2004) and in a channel with a local constriction and
expansion and compare the numerical results against laboratory results (Crowder et al.,
2004h).
6.4.8.1 Objectives
• Assess the ability of CanaiFlow to replicate the behaviour of a surge wave, caused
by the sudden collapse of a large body of water, in a channel with a local con-
striction and expansion;
• Benchmark the numerical results with laboratory results obtained in the Euro-
pean Commission's The European Concerted Action on Dam-Break Modelling
(CADAM) project; and
160
6.4 Test Cases
• Present the particulars for developing and undertaking the test using CanalFlow
and present the associated results.
In the first part of this test CanalFlow is setup to simulate an idealised dam-break in
a rectangular, unit width, horizontal and frictionless channel, 200 m long. The initial
conditions were hI = 10 m for the water depth upstream of the dam, which is located
in the middle of the channel, and h2 = 1 m is the water depth downstream of the
dam. The distance step dx equals 0.5 m. The results of this simulation are plotted,
as depth along the channel, at time t = 6 seconds. The data for this test is adapted
from Gottardi and Venutelli (2004).
A schematic diagram of the channel setup with contraction and expansion is given
in the Figure G.g. The model comprises of a horizontal channel of uniform 0.5 m wide
rectangular cross-section. The overall length of the channel is set up at 19.30 m, with
the first 6.10 m of the channel at the upstream end specified for the reservoir. The
actual CADAM model has a removable sluice gate built into the channel to retain
the water within the reservoir but which is not included in the CanalFlow model. A
constriction is located 7.70 m downstream of the gate, with an overall length of lA
m, the first and last 0.2 m of the constriction being tapered at 45 deg angles to the
channel walls. The initial conditions for the test are set at 0.3 m depth of water in the
reservoir upstream of the gate, and 0.003 m in the channel downstream of the gate.
The experimental data for the test was obtained by removing the gate, and measuring
the depth and velocity of flow at the benchmarking stations 81 to 84 (Figure 6.9). 81
is located 1.0 m upstream of the dam gate, 82, 6.10 m downstream of the dam gate,
83, 8.80 m downstream of the dam gate, and 84, 10.50 m downstream of the dam gate.
The measurements were taken every 0.04 seconds, and recorded up to a simulation time
of 10 seconds (Crowder et al., 2004h).
CanalFlow is set without a controlling sluice gate and the initial water levels and
flows specified are as given in the Appendix B. By setting up the model in this way
any stability/operational problems that may have been caused by the operation of the
161
6.4 Test Cases
~..•.•. :n
';" "
:;: rh 0.2 m
Io.5m
Upstream Reservoir
I·
0.2 m
81
+ + ,I
Gate
I.Om Q.2m
">. . . . s:~ ~ .:~.~ :
82 /
:
83
....
\
".
84
Ij>· 1;0;1
\.) ()
1 1 1 1
:j>. . . . ........ 1 . ·. ···<I
6.1 m
'I' 7.7 rn l.4m
t 19.3 m .1
Figure 6.9: Schematic illustration of the Contractions and expansions test (R.l). The
laboratory model consists of a 19.30 m long rectangular channel with a removable sluice
gate at 6.10 m and a constriction at 13.8 m. The constriction is 0.1 m wide and 1.0 m
long. Also shown are the gauging stations SI to S4 .
.sluice gate structure would be neglected. Manning's n throughout the channel is 0.01.
Interpolated cross-sections are not used for this test.
6.4.8.3 Dataset
The dataset used by CanalFlow to execute this test case is as given in the Table Rn
in Appendix R
6.4.8.4 Notes
i. Sluice gate is not modeled in this test case. Thus, instantaneous dambreak occurs
in a rectangular channel;
iii. A part from the above-mentioned test configuration this test case is run using
162
6.4 Test Cases
different limiters, namely, minmod, monotonised centered, superbee and van Leer
and without any limiter; and
iv. The results of this test case are discussed in the Section 7.11.2.
163
CHAPTER 7
7.1 Introduction
This chapter presents the results from applying CanalFlow to the test cases described in
the Chapter 6. Section 7.2 provides details of the testing platform, Section 7.3 presents
a general procedure to configure a test case using CanalFlow and the results of all test
cases listed in the Table 6.1, page 142 are discussed from Section 7.4 onwards. The
presentation of results and analysis carried out here is as per the guidelines given in the
benchmark test suite of the Environment Agency. Each test case is discussed in terms
of the model building process with CanalFlow, the simulation parameters and output
results compared with analytical and/or results of other softwares for the same test.
The time required to run these tests was recorded in order to get some idea about the
speed of the algorithm. The final conclusions and future development suggestions are
presented in the Chapter 8.
This section presents details of the hardware used to run the tests. The benchmark
tests were carried out on the laptop running on a Windows 2000 (Service Pack 4)
operating system and a Mobile lntel Pentium(R) 4, 3.06 GHz processor with 1 GB of
164
7.3 General Test Configuration Procedure
RAM. The Java code was compiled and run using JDK 1.5.0. The test data was stored
in a MySQL storage engine running on the same laptop with MySQL Connector / J
3.0.11 as a JDBC driver (see Section 504.2). The list of other tools used for the testing
is provided in the Table 5.1, page 119, while details of the database tables (metadata)
are provided in the Appendix A.
We will discuss the general procedure to configure a test case in CanalFlow. Before
discussing the general steps involved in building the model we refer to the main screen of
CanalFlow shown in the Figure 7.1. It is noted that this is a basic GUI of CanalFlow
and is developed just to enable the user to build a network and run the simulation
process. At this moment, CanalFlow does not provide any GUI for the data entry,
instead, it uses a front-end of MySQL storage engine to enter the data (see Figure i.2).
The whole task of running a test case could be divided in the following steps:
165
7.3 General Test Configuration Procedure
I- ~a~·':'~"l,,:,~~.;~~~",.:, ';;j~" n;';2"_'':D,1O:L>i..:JiLL''';; ~>(, ;....~:,:;:/;;;"'d,_; ;;I,';;:;';";''''''; ;:\.;/<li.~;;:;;;';;~;;'4. ,:-;_;';;' ,-I'.::,;;;:',,':", ,.;~;;'~,jw<O:;",,"}!:,;Jki;ii:.d
i';~,,~,c, ~C,G","G"G;',",;,"''';J'~ ;;,r;",;,;"" ~''''', er, ";'''"';,'''r''';",,, ";. ,;,.;",;";,:,.,:'";,,,,,,,',u" !
I
I N.twork $Chematlc:,Test N.1.1
I
i
I
G-1
.
.
• •
rmu1atkln SlIltUs
St. . TIrne: fNJ .... 1400~OOBST: ' Pl'ogr_ r 'Z -i" L ',' e
IrTWlaIlon messagasl
Im\llaUOn~!'Oee" started,
~.rung uulpu11tl ftJe OB. Please waIf:"
ca,' I O,I02726016l376058~ ImulaUon over.
F'IOU1S00D);OOSST 2( 181&1'(1 100k41Tll-nuie(s),anll'j,sacond(S)1I! ~cute 'ltlB iasi
Cl.nenlTI_ rrl . . . tSOO,OOIXlBST2I,
r-;;;:"","="""""=""",,.,
c.n __..:1
I ,l1li111 .....,
llmeSlep:: I O.1SJ6266S50616268 1018.11619:lJ:l12733
Figure 7.1: CanalFlow main screen. This is a basic GUr developed to help build the
network, run simulation and visualise the flow in a channel. It also shows run-time
simulation parameters such as current time, time step (At), Courant number, initial
and final mass of the system and overall progress. After the simulation is over, the
screen on the right hand side bottom corner displays the total time taken to execute
the task.
166
1.3 General Test Configuration Procedure
1111"._
111 .. b","".
8 r canaffloYi
I!I mrn_cM'l~
00 gm] m_t:hamef
ffi [[lJj m_ciTcular
,lEM .m_clivert
!E M m..netw6tk'
re illnl nt.~l~fu~~~
re I!!!Jim.JlOde.:
EH @jj m..,parahQla
m
,00. .I1'J.PlIr,~m
III I!!!Ji mJ>Uffip
[!( ffIDj ,mj:,U(nIi_cr.'/
!!l I!!!Ji "'Jad..o&e
ar l!H1l mJect_gate
8:1llm] m_~5:_curve'
!!l ~ m_r,eservoll"
!±I !!mi m.;.staQe'_dis
I±I mm m.;.structure,
Ii1 ~ mJr~pel
!!l1!!!Jirn..-
ltJ lffm tJnUiata
!!J ~ t_node~ary
!!l I!!!Jit_output.
!!l I!!ll '_«h'dui.
l!! 1!!!Ji'_voI
Figure 7.2: Front-end of the MySQL storage engine (refer Table 5.1, page j 19). The
screen is divided in two parts; the left hand part shows list of available tables and
the right hand part shows columns of the selected table in which data can be entered.
Shown is the example of network parameter table and corresponding data. The users
aware of SQL can also use SQL editor to insert/modify/delete the data. SQL scripts
are generated for building network and retrieving output data (which is then exported
to Excel sheet using MySQL Front) for each test case.
167
7.3 General Test Configuration Procedure
• Define network parameters like desired Courant number, time step and Iim-
iter to be used, start time and total simulation duration. The table structure
is as given in the Table A.5, page 231, while the typical data for a network
is shown in the Table 7.1;
• Initialise structures, if any, using the Table A.15, page 235. This will just
tell CanalFlow the number of structures and their location in a network
while the details of structures are to be provided in a table depending on a
structure type. For example the rectangular gate details are provided in the
Table A.H, page 2:33;
• Provide initial state in a channel by using the Table A.18, page 236;
iii. Execute the simulation process using Start sub-menu under the Simulation menu
(Figure 7.1);
iv. After the simulation is over user can retrieve the output data from the database
using SQL editor. In future, this process can be integrated in the GUI to shield
the user from nitty-gritties of the SQL. The data then can be exported to Excel
format from MySQL Front.
168
7.4 Test N.1: Sub-critical, Super-critical and Transitional Flows
As per the specifications in the Section 6.4.1, this test is an assessment of the ability of
CanalFlow to calculate sub-critical, super-critical and transitional flows and assess the
- -- -
numerical accuracy of the software package with reference to analytical results (Crowder
et al., 2004i). The analytical results for this test were obtained from the work done
by MacDonald (1994).
The analysis of results for the sub-parts N.1.1 to N.1.5 of the test case has been
done in terms of:
• The percentage difference between CanaIFlow and analytical solution for longi-
tudinal water surface profile; and
• Root Mean Square (RMS) error over the length of the reach for each part of the
test. CanalFlow RMS error is compared with RMS values of ISIS, HEC-RAS and
MIKE 11 run for the same test. RMS error is given by:
i=l
RMS= (7.1)
N
The analysis of results for the sub-part N.1.6 of the test has been done in terms of:
• The water surface profile generated by CanalFlow at the end of 6, 12, 18 and 24
hour.
This test has been undertaken in six separate sub-parts as per the specifications in
Crowder et al. (2004b). CanaIFlow is run in an unsteady mode with the quasi-steady
169
1.4 Test N.l: Sub-critical, Super-critical and 'Transitional Flows
boundary conditions (see Section 6.4.1.4). Initial conditions in terms of depth, dis-
charge, and elevations are specified and they are provided in Tables B.l to B.5 in
Appendix B. For all sub-parts of the test a variable time step is used which is based
on the desired Courant number of 0.80. Interpolated cross-sections are not used and
for all parts of the test dx = 1.0 m. The diagnostics file provided no warnings for any
parts of the test.
CanalFlow uses MoC to apply the boundary conditions. In case of a super-critical
flow, two characteristic curves have positive slope so that two boundary conditions
must be given at the upstream end (Garcia-Navarro and Saviron, 1992a). Thus, for
parts N.1.2, N.1.3 and N.1.5, discharge-depth hydrograph is specified at the upstream
end. Initial depth is provided along with upstream boundary discharge and extended
through to 01:00 hour respectively. The default network parameters set for all sub-parts
of the test N.1 are given in the Table 7.1.
Table 7.1: Network parameters set for the sub-critical, super-critical and transitional
flows (N. 1) test case.
Parameter Value
170
7.4 Test N .1: Sub-critical, Super-critical and Transitional Flows
Figures 7.3 to 7.5 present the results for all six sub-parts of the test N.1. The flow of
water is from higher elevation to lower elevation, i.e., from right to left in the figures.
i. The results for the sub-critical part (N.l.1) are illustrated in the Figure 7.3(A}.
This figure shows the water level profile and percentage error when compared
with the analytical solution. It is visible from this figure that the water levels
produced are smooth with no visual fluctuations and CanalFlow calculated water
level closely follow the analytical solution. The percentage error in water depth
when compared to the analytical solution, as illustrated in Figure 7.3(A}, is within
0% and -0.5 % and the RMS error is 0.002. CanalFlow took 4 minutes and 3
seconds to execute the test in a variable time step mode for simulation time of
24 hours. This time includes writing output to the database.
ii. The results for the super-critical part (N.1.2) are illustrated in the Figure 7.3(B}.
The results closely follow the analytical water surface profile. The percentage er-
ror in water depth when compared with the analytical solution is between -0.2 %
and -1.0 % with a slight fluctuation on the downstream side and the RM8 error
for this part is 0.012. Throughout the channel, CanalFlow has underestimated
the water levels by approximately -0.5 %. CanalFlow took 4 minutes and 40
seconds to execute the test in a variable time step mode for simulation time of
24 hours.
iii. For the super-critical to sub-critical flow (N.1.3) test CanalFlow clearly follows
the analytical water surface profile with smooth levels and no major fluctuations
as seen in the Figure 7.4(A}. The error in analytical and calculated solution is be-
tween ±0.2 %. This is far less than other commercially available softwares tested
in the EA benchmark. For example, 1818 and MIKE 11 under predicted water
level by up to approximately 8 % over the complete length of the channel (Crow-
der et al., 2004b). As seen in the Figure 7.6, CanalFlow RM8 error for this test is
approximately 6 times less than other softwares. CanalFlow took 5 minutes and
171
7.4 Test N.1: Sub-critical, Super-critical and 'fransitional Flows
~
"
iii 0.5 -1.0
0.0 -2.0
0 10 20 30 40 50 60 70 80 90
Chalnage (m)
(A)
4.0 2.0
:§: 3.0
-
1.0
j
c
:8
w 1.0
2.0
.~.
---=== 0.0
-1.0
C
~
0
t:
UJ
0.0 -2.0
o 10 20 30 40 50 60 70 80 90
Chalnage (m)
(B)
Figure 7.3: The figure shows bed elevation and comparison of longitudinal water surface
profile generated by CanalFlow and as given by analytical solution plotted on the
primary y-axis. It also shows the percentage error between CanalFlow and analytical
solution plotted on the secondary y-axis for test (A) sub-critical flow (N. 1.1 ) (B)
super-critical flow (N.1.2)
172
7.4 Test N.l: Sub-critical, Super-critical and Transitional Flows
2 . 5 . , - - - - - - - - - - - - - - - - - - - - - , - 2.0
0.5 t--=lIIfI~
0.0 -2.0
o 10 20 30 40 50 60 70 80 90
Chalnoge (m)
(A)
4.0 r - - - - - - - - - - - - - - - - - - . 10.0
:§: 3.0 ~=====::t--_::;;;--"':O::==::======+ 5.0
....
,gc 2.0 -jtfIfI_ _ _ _ ~
0.0 ""
"0
~ ~
w
W 1.0+--------.~----------------------------__+ -5.0
(B)
Figure 7.4: The figure shows bed elevation and comparison of longitudinal water surface
profile generated by CanalFlow and as given by analytical solution plotted on the
primary y-axis. It also shows the percentage error between CanalFlow and analytical
solution plotted on the secondary y-axis for test (A) super-critical to sub-critical flow
(N.1.3) and (B) sub-critical to super-critical to sub-critical flow (N.1.4)
173
7.4 Test N.l: Sub-critical, Super-critical and Transitional Flows
51 seconds to execute the test in a variable time step mode for simulation time
of 24 hours.
iv. The results for sub-critical to super-critical to sub-critical flow (N.1.4) test are
shown in the Figure 7.4(B). The Figure shows that the CanalFlow has captured
the hydraulic jump and the CanalFlow error as compared to analytical solution
is almost zero percent except at the point of an overestimation of water depth by
approximately 5 % near the area of jump. The HMS error for this part is 0.017
which is again approximately 6 times less than other softwares tested in (Crowder
et al., 2004b). CanalFlow took 4 minutes and 46 seconds to execute the test in a
variable time step mode for simulation time of 24 hours.
v. Results of the test N.l.4 are repeated in the super-critical to sub-critical to super-
critical flow test (N.l.5). From the Figure 7.5(A) it can be clearly seen that
CanalFlow has captured the jump reasonably well. The figure also shows that
CanalFlow has overestimated the water depth at a point upstream of the hydraulic
jump by roughly 10 %. Otherwise, the percentage error is in the range of 0.5 %
and -3.5 %. The HMS error for this sub-part is 0.027. CanalFlow took 4 minutes
and 39 seconds to execute the test in a variable time step mode.
vi. No analytical solution is available for the transitional flow (N.l.6) test case. Thus,
the longitudinal water profiles after 6, 12, 18 and 24 hour are shown in the Figure
7.5(B). The figure shows that CanalFlow is able to run the transitional flow test
and no instabilities are observed as the flow changes from the super-critical flow to
sub-critical flow along the upper section of the channel. CanalFlow took 5 minutes
and 4 seconds to execute the test in a variable time step mode for simulation time
of 24 hours.
174
7.4 Test N.l: Sub-critical, Super-critical and Transitional Flows
4.0 10.0
]: 3.0 5.0
~~
~ C
--
~ 2.0 0.0 ~
g
~ w
iii 1.0 -5.0
0.0 -10.0
o 10 20 30 40 50 60 70 80 90
Chalnago (m)
(A)
:8m 1.5
~ 1.0
../'
w ~
0.5
0.0
0 10 20 30 40 50 60 70 80 90
Chainago (m)
(B)
Figure 7.5: (A) Bed elevation and comparison of longitudinal water surface profile
generated by CanalFlow and as given by analytical solution plotted on the primary
y-axis and percentage error between CanalFlow and analytical solution plotted on
the secondary y-axis for super-critical to sub-critical to super-critical flow test (N.1.5)
(B) Longitudinal water surface profile at the end of 6, 12, 18 and 24 hour. The flow
boundary at the upstream side (on the right) is kept constant while at the downstream
side (on the left) the depth boundary is kept constant for first 6 hours and then increased
linearly from 0.602 m to 2.0 m till 18:00 hour and then kept at that level for another
6 hours.
175
1.4 Test N.l: Sub-critical, Super-critical and Transitional Flows
~=
~
0
,fi 0.080
'i
!!,l '>
a:
0.060
0.040
" - ",- " I,: _
0.020
0.000
N.1.1
".
IT:'
N.1.2
.
N.1.3 - N.1.4
L '1=
c·. .
N.1.5
Sub-Parts of Test Case N.1
Figure 7.6: Comparison of RMS errors for five sub-parts of the sub-critical, super-
critical and transitional flow (N. 1) test case. It is evident from the bar graph that
amongst ISIS, HEC-RAS, MIKE 11, and CanalFlow, CanalFlow has the lowest RMS
error for all parts of the test. The average RMS error is approximately 4 times less
than that of ISIS and HEC-RAS and by approximately 6 times less than MIKE 11.
116
7.5 N.02: Thiangular Channels
As per the specifications in the Section 6.4.2, this test is an assessment of the capability
of CanalFlow to calculate the normal sub-critical flow depth and the normal super-
critical flow depth in a triangular channeL As discussed earlier in the Section 7.4.1, the
depth-discharge hydrograph boundary condition is used on the upstream side of the
channel in a super-critical flow test.
The analysis of results for parts N.2.1 (sub-critical fiow) and N.2.2 (super-critical
flow) of the test has been done in terms of:
• Comparison of longitudinal water surface profile with the normal depth. The
normal depth is determined based on the Manning's equation:
3
• Longitudinal profile showing the percentage difference between normal depth and
longitudinal water surface profile.
To execute the test N.2, CanaIFlow was run in an unsteady mode with the quasi-
steady boundary conditions (see Section 6.4.2.4). Initial conditions in terms of depth,
discharge, and elevations are as provided in the Table B.6 in Appendix B. Interpolated
cross-sections were not used and the diagnostics file provided no warnings for any parts
of the test. All other simulation parameters were as listed in the Table 7.1.
The results for both the sub-critical and super-critical parts of the test N.2 are presented
in Tables 7.2 and 7.3 and in the Figure 7.7.
177
7.5 N.02: 'Iriangular Channels
i. The Table 7.2 presents the water level data generated by CanaIFlow and other
softwares used by EA benchmark suite. From the Table 7.2 it is clear that
CanalFlow has produced an exact match to the normal depth from chainage
900 m and above. The Figure 7.7(A) compares the normal depth determined by
Manning's formula (Equation 7.2) and the one generated by CanalFlow. The per-
centage error seen from the figure is in range of -0.01 % to -0.06 %. CanalFlow
was run in a variable time step mode and it was observed that the approximate
time step for desired Courant number to be 0.80 was 48 seconds and it took 1 sec-
ond for CanalFlow to execute the task for the total simulation time of 24 hours.
1~1~#±:Sgl
o 15 30 45 60 75 90 105 120 135 150
Chlllnaga (m)
(A) (D)
i~~~Ia;{::~:
.2.00~ 1
:~::~:
2
::13 4 5 6 7
ChaJnaga (m)
8 9 10 11
i J1~f$:: 1 2 3 4
:::::::1
5 6 7
Chalnage (m)
8 9 10 11
(C) (D)
Figure 7.7: Comparison of normal depth and CanalFlow depth for (A) the sub-critical
and (B) the super-critical part of the test (C) shows the percentage error in depths
generated by CanaIFlow, ISIS, HEC-RAS, MIKEll and the normal depth for the sub-
critical flow, (D) shows the similar percentage error for the super-critical flow test case.
ii. Table 7.3 shows that CanalFlow results are almost identical to the normal depths
178
7.5 N.02: Triangular Channels
179
7.5 N.02: Triangular Channels
Table 7.2: Triangular channel test (N.2): Comparison of water levels generated by
CanalFlow, ISIS, REC-RAS and MIKE 11 run in a resistance radius (RR) and hydraulic
radius (RR) mode for the sub-critical flow test case (N .2.1). The results for softwares
other than CanalFlow presented in this table and Table 7.3 are adapted from Crowder
et al. (2004d).
Water Level (m)
Chalnage Normal ISIS HEC-RAS MIKE 11 MIKE 11 CanalFlow
(m) Depth (RR) (HR)
(m)
0 13.012 13.000 13.000 13.000 13.000 13.004
300 13.312 13.306 13.305 13.201 13.305 13.310
600 13.612 13.608 13.601 13.436 13.608 13.613
.00 13.912 13.910 13.908 13.698 13.910 13.912
1200 14.212 14.211 14.209 13.977 14.211 14.212
1500 14.512 14.511 14.509 14.266 14.511 14.512
1800 14.812 14.812 14.810 14.560 14.812 14.812
2100 15.112 15.112 15.110 14.857 15.112 15.112
2400 15.412 15.412 15.410 IS. IS!') 15.412 15.412
2700 15.712 15.712 15.710 15.455 15.712 15.712
3000 16.012 16.012 16.010 15.754 16.012 16.012
Table 7.3: Triangular channel test (N.2): Comparison of water levels generated by
CanalFlow and other softwares used in the EA benchmark suite for the super-critical
flow test case (N.2.2).
Water Level (m)
Chalnage Normal IS IS HEC-RAS MIKE 11 MIKE 11 CanalFlow
(m) Depth (RR) (HR)
(m)
0 11.718 11.7 11.7 11.7 11.7 11.731
15 12.018 12.01 12.009 11.935 12.008 12.018
30 12.318 12.314 12.312 12.2 12.312 12.318
45 12.618 12.615 12.613 12.483 12.614 12.618
60 12.918 12.916 12.913 12.775 12.915 12.918
75 13.218 13.216 13.214 13.071 13.216 13.218
.0 13.518 13.516 13.514 13.369 13.516 13.518
105 13.818 13.816 13.814 13.668 13.816 13.818
120 14.118 14.116 14.114 13.968 14.116 14.117
135 14.418 14.416 14.414 14.268 14.416 14.413
150 14.718 14.716 14.714 14.568 14.716 14.710
180
7.6 N.3: Ippen Wave
As per the specifications in the Section 6.4.3, this test is a comparison of results gener-
ated by the CanalFlow with an analytical solution based on the hydrodynamic theory
of tidal wave propagation in a horizontal channel of uniform cross-section and finite
length (see Section 6.4.3). The analysis of test results is done in terms of:
• Comparison of calculated water levels and velocities against the analytical values
over the duration of the simulation at the closed boundary and at locations of 25
km, 50 km and 75 km downstream of the closed boundary; and
• Comparison of the calculated longitudinal water level and velocity profiles with
the analytical solutions over the length of the channel at the fixed simulation
times of 2.0, 2.5, 3.0, 3.5 and 4.0 hours.
The simulation was run for 4 hours (Crowder et al., 2004f) and the CanalFlow simula-
tion parameters for this test were as given in the Table 7.1. The initial channel data
used is presented in the Table B.7, in Appendix B. Approximate time step for the
desired Courant number, 0.8, was 26 seconds and CanalFlow took 3 seconds to execute
the test. Initial volume of the system shown by CanaIFlow was 2.0E + 9 m 3 and at the
end of the simulation it was 1.9994E + 9 m 3 . This shows thats the mass is conserved
in a closed system with no inflow boundary.
The results for Ippen wave test are presented in terms of water level and velocity time
series and profiles in Figures 7.8 to 7.10. The analytical solution for this test was
obtained from the work done by Ippen (1966).
i. Analysis of water level and velocity time series is done in the Figure 7.8. It shows
that the CanalFlow calculated water levels are generally very close to the analyti-
cal values, and thus follow the expected trend, but it does not follow the analytical
181
7.6 N.3: Ippen Wave
values exactly at any location for the duration of the test. Variation in the an-
alytical and calculated water levels is observed after 2.5 hours of the simulation
at 25, 50 and 75 km (Figure 7.8(C,E and G)). The reason for this discrepancy
might be due to the fact that the actual length of the wave is slightly longer
than the length of the channel (by 800 m). Therefore, as the wave is supposed
to travel beyond this point within the first 2 hours, it would be expected that
the water levels would be slightly different even without any reflection. However,
as the wave will have reached the closed boundary slightly before 2 hours of the
simulation time, some reflection will have already occurred despite the fact that
none should have, and the resultant superimposition will have forced the water
levels to rise even further (Crowder et al., 2004f).
ii. Analysis of water level and velocity profiles is done in the Figure 7.9 and 7.10.
The standing wave of Ippen is the one which would result from the total super-
imposition of the reflected wave with the incident wave after one time period of
the reflected wave and from Ippen's mathematical solution it can be seen that
the harmonic function for the wave amplitude applied at the open end of the
channel follows the cosine relationship (Crowder et al., 2004f). In comparison,
the harmonic function applied at the downstream end of the model followed a
sine curve relationship, as confirmed by the shape of the curve representing the
water levels along the length of the channel at only 2 hours. For this reason, the
generated wave formulation is out of phase with Ippen's theory for tidal waves
entering channels of finite length (Ippen, 1966).
As the applied water levels at the open end of the channel follow a sine wave
relationship as opposed to a cosine relationship, the resultant phase difference is
a 90 deg time angle, equivalent to a time lag of 30 minutes. Therefore, in order
to achieve the final standing wave that Ippen refers to, it would be necessary to
run the simulation for 4.5 hours as opposed to 4.0 hours. Total superimposition
of the reflected wave with the incident wave would still be achieved after only 4.0
hours. However, no greater effect of further reflection and superimposition should
182
7.6 N.3: Ippen Wave
i~F:'s2d
1~. •
Jm·····,,·····,······t
,.• ,.• _ow, -". , ,.• ,
" " " r-lI1rs) " •
(A) (B)
/:
I__ ~ -e-c.n.IF;;] ' _ _ Ano/yIICIII-e-c.n..flow I
f~F·
.. .•
,.• , ,.
,..",.(11..)
. . 1-
irl ?§d..
,.• ••
I~l.....
fM
> -0.10
-UD
,.• ,.• ,.•
T1""'(II")
,.•
(C) (D)
,.• 2.5
-On,
(E)
3.0 3.5 •.• ,.• 2.5 3.0
TI_tll ..)
(F)
3.8
!
••
f~@?f":
.•
, .•
-On,
~+g H~~;zP4
.. ,. . .•
"
,
"
-,~
,
"
TI_(H ..)
, ••
(C) (H)
Figure 7.8: Ippen wave test: Water level and velocity time series at (A) and (B) 0 km,
(C) and (D) 25 km, (E) and (F) 50 km, and (G) and (H) 75 km. The Figure (C) and
(E) shows that the calculated water levels follow analytical solution closely between
2.0 and 2.5 hours and after that it break from the analytical solution. At 75 km (G)
calculated solution follows the analytical solution between 3.0 and 3.5 hours and after
that variations are seen. The Figure (B,D,F and H) shows that the velocities produced
by CanalFlow follow the analytical solutions fairly closely. In (D) and (F) it is seen
that the velocities are identical to analytical solution till 2.4 hours and after that they
start breaking up. While in (H) they are identical between 2.5 hours and 3.5 hours.
183
7.6 N.3: Ippen Wave
occur at anytime thereafter as the reflected wave simply travels out of the open
boundary after 4.0 hours.
The Figure 7.9{E) and 7.1O{A), provide further evidence of the correct superim-
position of the reflected wave with the incident wave. At 3.0 hours the reflected
wave will have travelled half the length of the channel, but would be out of phase
184
7.6 N.3: Ippen Wave
I ~:~'_20~ ~
i~t ~
t====:::;;;---~.~t;;;;;~
0 20
~40 60 80 100
f;:~ ~-----=-=-I
>-0.80 _
_1.20.1.-_ _ _ _ _ _ _ _ _ _ _ _- '
o 20 40 60 80 100
Chainage (m) Chalnsg. (m)
(A) (B)
1
~'().80t========:::-==~===j
-1.20
o 20 40 60 80 100 o 20 40 60 80 100
Chalnsge (m) Chainsge (m)
(C) (D)
(E) (F)
Figure 7.9: Ippen wave test: Water level and velocity profiles after (A) and (B) 2 hours,
(C) and (D) 2.5 hours, (E) and (F) 3 hours.
185
7.6 N.3: Ippen Wave
o 20 40 60 80 100 o 20 40 60 80 100
Chalnage (m) Chalnage (m)
(A) (B)
I~ r-l--....~1iijii::1
o 20 40 60 80 100 o 20 40 60 80 100
Chainag. (m) Chalneg. (m)
(C) (D)
Figure 7.10: Ippen wave test: Water level and velocity profiles after (A) and (B) 3.5
hours, (C) and (D) 4 hours.
with the incident wave by 180 deg. As such, it would be expected that the two
waves would cancel each other out up to channel chainage 50 km, whilst the
water levels between 50 km and the open boundary would be representative of
the remainder of the incident wave. These observations are confirmed in the Fig-
ure 7.9(E). The water levels are approximately horizontal up to chainage 50 km,
showing complete cancellation of the two waves. After chainage 50 km, the water
level then returns to the sine wave formation of the incident wave as expected.
The same inference is demonstrated in the Figure 7.1O(C), showing the calculated
water levels at 4.0 hours. As the reflected wave will have returned to the open
end of the channel at this time, though again be out of phase with the incident
wave by 180 deg, it would be expected that both waves would cancel each other
out over the entire length of the channel. It is clear from the graph that this has
happened to a certain extent, the slight remaining water levels being due again to
186
7.6 N.3: Ippen Wave
the damping effect of the channel bed. Essentially, the amplitude of the returning
wave would be less than that of the incident wave as it has travelled twice the
length of the channel by the time it reaches the open boundary.
The wave that Jppen describes is the result of an initial cosine wave, as opposed
to the sine wave applied to this model, it would be expected that the points of
maximum velocity, at distance of L/4 and 3L/4, would occur at the simulation
times of 2.5 and 3.5 hours. Furthermore, this would be consistent with the ob-
servations made earlier that there was no amplitude of the water surface about
the mean water level at these corresponding time and locations. However, Fig-
ures 7.9(B,D and F) and 7.10(B and D) show that maximum positive or negative
velocities are achieved at these channel chainages, but not at the same times zero
amplitudes are recorded at these locations. Instead, maximum velocities are gen-
erated at the approximate locations after 2.0, 3.0 and 4.0 hours. Zero velocities
are also recorded at these times at the closed boundary and at multiples of L /2
chainage in the downstream direction, which correspond with the locations of
zero amplitude.
All calculated and analytical velocities, at any given time and location, can be
observed to switch between positive and negative values. Positive velocities refer
to flow travelling towards the open boundary, and negative velocities towards the
closed boundary. It can therefore be deduced from the Figures 7.9(B,D and F)
and 7.10(B and D), that positive values of velocity correspond with low tide water
levels, negative velocities with high tide levels, and zero velocity when the water
level returns to mean water level.
The deviations seen in the calculated and analytical velocities might again be
due to the marginal difference in the length of the wave and the length of the
channel. As there is some premature reflection at the closed boundary, forcing the
water levels slightly above the expected values at the start of the simulation, the
resultant imbalance of water within the system could give rise to slightly varied
velocities.
187
7.7 N.4: Monoclinal Wave
As discussed in the Section 6.4.4, this test is an assessment of the ability of CanalFlow
to recreate the special case of an unsteady flow, known as the monoclinal rising wave,
as a typical case of uniformly progressive flow (Crowder et al., 2004g).
The results of Monoclinal test case are analysed in terms of:
The results are compared with MIKE 11 as MIKE 11 is the only other software cited
in Crowder et al. (2004g) which uses Chezy's relation as specified by the test. Other
softwares namely ISIS and HEC-RAS have used equivalent Manning's relation.
The data used for this test case is presented in the Table B.8, Appendix B while the
simulation parameters used are shown in the Table 7.4. The change in these parameters
over the default parameters is the time step, which is kept fixed at 60 seconds and use
of Chezy's resistance law instead of Manning's (as per the test specification in Crowder
et al. (2004g)). Additionally this test case was run using different limiters, namely,
Minmod, Monotonised centered, Superbee and van Leer. For fixed time step of 60
seconds, it took 3 seconds for CanalFlow to execute the test case.
The results of the simulation are presented in the Figures 7.13 to 7.1S. The figures show
the analytical water surface profile of the wave entered at the start of the simulation,
together with successive profiles of the CanalFlow calculated water levels at 5, 10, 15
and 20 hours. It also shows the percentage error calculated by comparing CanalFlow
and MIKE 11 output with that of analytical solution plotted on the secondary y-axis.
188
7.7 N.4: Monoclinal Wave
Table 7.4: Network parameters set for the Monoclinal wave (N.4) test case. The limiters
used are Minmod (MM), Monotonised centered (MC), Superbee (SB) and van Leer
(VL).
Parameter Value
The minimum and maximum errors as compared to analytical solution using various
limiters are presented in the Table 7.5.
189
7.7 N.4: Monoclinal Wave
]: 6.0 \ C
=go
c
5.0 ~.-~~-~------~
~,- \ j'
~-----.,~~~'~~------~ -2,0
0.0 ~
g
w
4.0 ~
3.0 t--------""'--=========t -4,0
2.0 " -_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _-1.
-6,0
1,OE+06 9,OE+05 8.0E+05 7.0E+05 6,OE+05 5,OE+05
Chainago (m)
(A)
2.0 -6.0
1.0E+06 9,OE+05 6.0E+05 7.0E+05 6.0E+05 5.0E+05
Chalnago (m)
(D)
190
7.7 N.4: Monoclinal Wave
(A)
8.0 4.0
7.0
" f'V... 2.0
:§: 6.0
:5
.Yf\.. 0.0
~
g
c
go 5.0 ~\ -2.0 w
4.0
3.0
\. -4.0
2.0 ..s.0
1.0E+06 9.0E+05 8.0E+05 7.0E+05 6.0E+05 5.0E+05
Chalnago (m)
(B)
191
1.1 N.4: Monoclinal Wave
\ 1\
2.0
0.0
-g
~
c•
S.O
""\1. ·2. o w
4.0
3.0
\. -4. o
2.0 ·s.o
1.0E+06 9.0E+OS 8.0E+OS 7.0E+05 6.0E+OS S.OE+OS
Chainage (m)
(A)
2.0 ·6.0
1.0E+06 9.0E+OS 8.0E+OS 7.0E+OS 6.0E+OS S.OE+OS
Chainage (m)
(BJ
192
7.7 N.4: Monoclinal Wave
2.0 ·6. o
1.0E+06 9.0E+OS 8.0E+05 7.0E+OS 6.0E+OS S.OE+OS
Chalnoge (m)
(A)
8.0 i==========~-__
~=------l- 4.0
_ 7.0 I-----------~__:~ I"<.~\__+ 2.0
.5. 60
;;
g.S.O
c
t=====~~-~~\-7~ 0.0 0t:~
.~
f------------------=~\L--+ .2.0 w
. . .Y\.
. c
4.0
3.0 } - - - - - - - - - - - - - - - - - - - - - " - = = t -4.0
2.0 ·6.0
1.0E+06 9.0E+OS 8.0E+05 7.0E+05 6.0E+OS S.OE+OS
Chalnoge (m)
(B)
193
7.7 N.4: Monoclinal Wave
9.0 B.O
8.0 4.0
:5
7.0
:§: B.O
go 5.0
c
"'"
~
\ "'-
2.0
0.0
-2.0 w
~
g
4.0
3.0
\\J -4.0
2.0 -B.O
1.0E+OB 9.0E+05 B.OE+05 7.0E+05 B.OE+05 5.0E+05
Ch_In_go (m)
(A)
2.0 -B.O
1.0E+OB 9.0E+05 B.OE+05 7.0E+05 B.OE+05 5.0E+05
Chainage (m)
(B)
194
7.7 N.4: Monoclinal Wave
-4. o
3.0
2.0 -6. o
1.0E+06 9.0E+OS 8.0E+05 7.0E+OS 6.0E+OS S.OE+OS
Chalnago (m)
(A)
195
7.7 N.4: Monoclinal Wave
M 4~
2 . 0 . 1 - - - - - - - - - - - - - - - - - - - - ' - -6.0
1.0E+06 9.0E+05 B.OE+05 7.0E+05 6.0E+05 S.OE+OS
Chainoge (m)
(A)
196
7.7 N.4: Monoclinal Wave
g
\ "<
Q f---------------~----_+
w
·2.0
4.0
3.0 t--------------'--====t -4.0
\
2.0 L.._ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _...L
·6.0
1.0E+06 9.0E+05 8.0E+05 7.0E+05 6.0E+05 5.0E+05
Chalnage (m)
(A)
2.0 -6.0
1.0E+06 9.0E+05 8.0E+05 7.0E+05 6.0E+05 5.0E+05
Chainage (m)
(B)
197
7.8 C.1: Looped System
Table 7.5: Minimum and maximum errors between CanalFlow and analytical solution
using various limiters. The limiters used are Minmod (MM), Monotonised centered
(MC), Superbee (SB) and van Leer (VL). The data clearly shows that the van Leer
limiter has a minimum range between minimum and maximum error while Superbee
has got maximum error on the negative side and Minmod has got maximum error
on the positive side with an exception at 10 Hr where Monotonised centered has got
maximum positive error. The last row in the table shows MIKE 11 error for the same
test case.
Limiter 5 Hr 10 Hr 15 Hr 20 Hr
Min Max Min Max Min Max Min Max
MM -2.80 1.52 -0.44 3.94 -0.52 2.84 -1.89 2.46
MC -2.17 0.17 -0.42 4.03 -1.11 1.10 -1.89 0.15
SB -3.86 0.19 -2.61 3.80 -3.52 0.58 -4.90 0.27
VL -2.34 0.40 -0.42 3.78 -0.49 0.91 -1.92 0.44
MIKE 11 Results (Personal Communication - Crowder, 2005)
MIKE 11 -1.93 1.26 -0.80 5.57 -0.87 4.24 -0.94 2.94
As described in the Section 6.4.5, this test is an assessment of the capability of the
software package to calculate a diverging and converging flow scenario, i. e., a looped
system (Crowder et aI., 2004c). The test is carried out for two boundary conditions,
QS1 and QS2 (see Section 6.4 ..5). For both QS1 and QS2 simulations, results for
discharge and water level results have been reported upon in tabular form at each of
the cross sections defined by the dataset. In addition, the discharge and water level
results in each channel are presented in a graphical form.
198
1.8 C.l: Looped System
CanalFlow was set to run this test for the network as shown in Figure 6.6, page 153.
Channel B was divided into Bl and B2 of equal length, 750 m. This was done due
to the way CanalFlow handles the channel geometry as the current version does not
handle curved channels. The simulation parameters are same as presented in the Table
7.1. For the desired Courant number to be 0.8, dt for this test case was approximately
equal to 0.53 and it took 5 minutes and 23 seconds for CanalFlow to execute the task
for QS1 and 5 minutes and 25 seconds for QS2 with dx = 5.0 m for both QS1 and QS2
and for all four channels. The initial conditions for all channels are specified as given
in the Table B.9, Appendix B. Junctions are specified to model the convergence and
divergence of flow (see Section 4.6.2, page 98).
The Figures 7.19 provide results generated by CanalFlow and other softwares in the
EA benchmark study for the QS1 boundary condition. The summary of results is
also provided in the Table 7.6. As seen in the Figure 7.19, almost all softwares have
predicted different flow parameters in the network. CanalFlow results in predicting
depths and discharges in Channel B and Channel C are close to REC-RAS and upto
certain extent to ISIS. Slight differences were observed in Channel A and Channel D
and the variations were in the range of 0.02 % to 0.03 % as compared to the flows
predicted by other softwares in these reaches.
The Figure 7.20 and Table 7.7 provide results of CanalFlow and other softwares
used in EA benchmark study for the QS2 boundary condition. Variations are seen in
depths and discharges in all channels as compared to other softwares. The variations in
Channel A and D are smaller as compared to QS1. The range of variation in discharge
in A and D is -0.008 % to 0.02 %. As compared to other softwares, CanalFlow has
predicted higher water levels in all the reaches. In terms of discharge, it has assigned
less discharge to Channel B and more to Channel C as compared to other softwares.
The summary of results is provided in Table 7.7.
199
7.8 C.l: Looped System
3 4 , 3
Chain_g. (mj
4 ,
Chalnage (m)
(A) (B)
-B-ISIS --+-HEC-RAS -M-MIKE11 (RR) -B-ISIS -+-HEC-RAS -M-MIKE11 (RR)
-6-MIKE11 (HR) _CanelFlow -6-MIKE11 (HR) _CanalFlow
(C) (D)
-B-ISIS --+-HEC-RAS -M-MIKE'1 (RR) -6-ISIS --+-HEC-AAS -M-MIKE11 (RR)
-6-M1KE11 (HR) ___ canslRow -iL-MIKE11 (HR) _CaI\llIF!ow
f>
500 375 250 125 o 500 375 250 125 o
Chalnage (m) ChalnlgB (m)
(E) (F)
3.08
-
260"1t - - - - - - - - - - - - - - - - - 1
250.3
250.1t~.~::::i~::::~i::::~i~::~It=~
iI 249.91
:[3.06
•
iil• 3 .04
3.02
Ii: 249.7 t·- - - - - - - - - - - - - - - - - 1
3.00 ~
249.5 +---_--_--_--_---<
100 75 50 25 0 375 250 125 o
Ch,ln8ge (m) '00
Chalnage (m)
(G) (11)
Figure 7.19: Looped system (C.1) test: Flow and stage results for the QSl boundary
condition. Shown are the water level (stage) and flow results (A) and (B) for Channel
A, (C) and (D) for Channel B, (E) and (F) for Channel C, (G) and (H) for Channel D.
200
7.S C.I: Looped System
(A) (B)
~~-·~T~O~.~~~~;~I
-lr-MIKE11 (HR) ___ CBnalFlow
2."
2.30 ...
-
:[2.10 li§:2-+-1
',.90
1.70
115.0 . .
1.50
,.,. 1125
ChliJig. (m) 375
• 1500 1125 750
Chatne". Cm)
375
•
(C) (D)
-e-ISIS -40-HEC-RAS -M-MIKEl1 (RR)I -B-ISIS -+-HEC-RAS ~MIKE11 (RR)
---tr-MIKE11 (HR) -+-CBnBIFIDw -6-MIKE11 (HR) -+-CanalFlow
:±;--=3--::~:~·;~~~d;
2."
2.30
Ill.
-- I§::bl;:E-
:[2.10
k,·90
1.70
1."
129.0
13 ~ ~ ~ ~ I
'
500 ", 2"
Chalnage (m)
125
• 500 375 250
ChIoINlp(m)
125
•
(E) (F)
f'" .
249.7 +-~~~~~~~~~~~-~~~-I ~"'lK
249.5+-~~~~~~~~~~~~~~~
100 75 50
Chaln.ge (m)
25 500 375 250
CtMIh,.ge(m)
12. •
(G) (H)
Figure 7.20: Looped system (C.1) test: Flow and stage results for the QS2 boundary
condition. Shown are the water level (stage) and flow results (A) and (B) for Channel
A, (C) and (D) for Channel B, (E) and (F) for Channel C, (G) and (H) for Channel D.
201
7.8 C.1: Looped System
Table 7.6: Looped system (C.l) test: Flow and stage results for the QSl boundary
condition. Flow is expressed in m 3 / sec while chainage and stage is in m. HEC-RAS,
ISIS, MIKE 11 results are adapted from Crowder et al. (2004c).
HEO-RAS IBIS MIKE 11 (RR) MIKE 11 (HR) CanalFlow
Channel Chalnage Flow Stage FI~ Stage Flow Stage FI~ Stage Flow Stage
A 100 250.00 3.065 250.00 3.066 250.00 3.072 250.00 3.100 250.07 3.066
75 250.00 3.062 250.00 3.063 250.00 3.069 250.00 3.096 250.07 3.062
50 250.00 3.059 250.00 3.060 250.00 3.066 250.00 3.093 250.07 3.059
25 250.00 3.056 250.00 3.051 250.00 3.063 250.00 3.089 250.07 3.056
0 250.00 3.052 250.00 3.053 250.00 3.060 250.00 3.086 250.10 3.052
B 1500 78.15 3.052 80.99 3.053 92.95 3.060 102.70 3 ..086 77.68 3.053
1125 78.15 3.039 80.99 3.041 92.95 3.044 102.70 3.{}63 77.68 3.040
750 78.15 3.028 80.99 3.030 92.95 3.031 102.70 3.043 77.69 3.029
375 18.15 3.019 80.99 3.022 92.95 3.020 102.10 3.027 11.68 3.020
0 78.15 3.012 80.99 3.015 92.95 3.011 102.70 3.-014 11.71 3.013
C 500 171.85 3.052 169.01 3.053 157.05 3.060 147.30 3.086 112.34 3.052
375 111.85 3.036 169.01 3.038 151.05 3.038 147.30 3.058 112.35 3.036
250 171.85 3.025 169.01 3.027 157.05 3.024 147.30 3.039 112.31 3.025
125 171.85 3.011 169.01 3.020 151.05 3.016 141.30 3.025 112.38 3.018
0 171.85 3.012 169.01 3.015 157.05 3.011 147.30 3.014 172.52 3.011
D 100 250.00 3.012 250.00 3.015 250.00 3.011 250.00 3.014 250.06 3.013
75 250.00 3.009 250.00 3.011 250.00 3.009 250.00 3.010 250.06 3.009
50 250.00 3.006 250.00 3.007 250.00 3.006 250.00 3.007 250.06 3.007
25 250.00 3.003 250.00 3.004 250.00 3.003 250.00 3.003 250.07 3.003
0 250.00 3.000 250.00 3.000 250.00 3.000 250.00 3.000 250.09 3.000
202
7.8 C.l: Looped System
Table 7.7: Looped system (C.1) test: Flow and stage results for the QS2 boundary
condition. Flow is expressed in m 3 / sec while chainage and stage is in m. HEC-RAS,
ISIS, MIKE 11 results are adapted from Crowder et al. (2004c).
HEC-RAS ISIS MIKE 11 (HR) MIKE 11 (HR) CanalFlow
Channel Chainage FI~ Stage FI~ Stage FI~ Stoge FI~ Stage Flow Stage
A 100 250.00 2.420 250.00 2.407 250.00 2.370 250.00 2.412 249.99 2.431
15 250.00 2.405 250.00 2.389 250.00 2.354 250.00 2.396 249.99 2.416
50 250.00 2.389 250.00 2.371 250.00 2.338 250.00 2.381 249.99 2.399
25 250.00 2.372 250.00 2.354 250.00 2.321 250.00 2.365 249.99 2.383
0 250.00 2.355 250.00 2.334 250.00 2.305 250.00 2.349 249.98 2.370
B 1500 117.61 2.355 120.50 2.334 116.85 2.305 117.80 2.349 115.80 2.367
1125 117.61 2.208 120.50 2.185 116.85 2.159 117.80 2.198 115.80 2.234
150 117.61 2.062 120.50 2.037 116.85 2.014 117.80 2.046 115.81 2.108
315 117.61 1.918 120.50 1.889 116.85 1.873 117.80 1.894 115.81 1.989
0 117.61 1.778 120.50 1.743 116.85 1.734 117.80 1.740 115.82 1.882
C 500 132.39 2.355 129.50 2.334 133.15 2.305 132.20 2.349 134.20 2.365
315 132.39 2.207 129.50 2.186 133.15 2.159 132.20 2.197 134.20 2.218
250 132.39 2.060 129.50 2.037 133.15 2.015 132.20 2.046 134.19 2.080
125 132.39 1.915 129.50 1.890 133.15 1.873 132.20 1.893 134.18 1.963
0 132.39 1.778 129.50 1. 743 133.15 1.734 132.20 1. 740 134.23 1.881
D 100 250.00 1.778 250.00 1.743 250.00 1.734 250.00 1. 740 250.00 1.878
250.02 1.835
"'
50
250.00
250.00
1. 733
1.687
250.00
250.00
1. 709
1.675
250.00
250.00
1.702
1.669
250.00
250.00
1. 707
1.672 250.05 1.783
25 250.00 1.642 250.00 1.638 250.00 1.635 250.00 1.636 250.13 1.715
0 250.00 1.600 250.00 1.600 250.00 1.600 250.00 1.600 249.80 1.638
203
1.9 C.2: Weirs
As discussed in the Section 6.4.6, this test is an assessment of the capability of the
software package to model a broad crested weir (Crowder et al., 2004e). The test is
run with a steady state and an unsteady state boundary conditions.
The analysis of results of the test has been limited to a
• Comparison of stage (water level) at the cross sections immediately upstream and
downstream of the weir structure and head loss over the structure for a steady
state boundary condition; and
• Comparison of stage (water level) verses time at the cross sections immediately
upstream and downstream of the weir structure for an unsteady case.
CanalFlow was set for the test configuration as shown in the Figure 6.7, page 155 and
run using the default simulation parameters provided in the Table 7.1. In a variable
time step mode, CanalFlow took 1 second and 3 seconds to execute the test for steady
and unsteady state case, respectively.
The steady state results from the Broad Crested weir test for free (SSl) and drowned
(SS2) flow conditions are presented in the Table 7.8. The unsteady results (transient
conditions) are presented in the Figure 7.21. These results show that the CanalFlow
can handle weir with steady and unsteady downstream boundary conditions.
204
7.9 C.2: Weirs
Table 7.8: Weir flow (C.2) test case: Steady state results for free and submerged flow
condition. The output shows that the head loss in SSl is equal to the ISIS and for SS2
it is closer to ISIS than MIKE 11 and HEC RAS.
0.'
EO.a
i
...I
0.7
0.6
./' "'"
~ 0.5 /
;: 0.4 ./'
0.3
5 9 13
TIme (hrs)
17 21 5 • 13
Time (hrs)
17 21
(A) (B)
Figure 7.21: Weir flow (C.2) test case for an unsteady boundary condition: Figure
shows water level time series at (A) downstream side of the Channel A and upstream
side of the broad crested weir and (B) upstream side of the Channel B and downstream
side of the broad crested weir.
205
7.10 C.3: Irrigation Network
• Discharge and depth hydrographs at upstream side of the channel C-1 and down-
stream side of the channel 0-4; and
• Discharge hydrograph and mass balance at junctions J-1, J-2 and J-3.
The canal irrigation network for this test case was configured as shown in the Figure
6.8, page 158. The channel geometry data is as shown in the Table RIO, Appendix
B while the network parameters are as shown in the Table 7.1. The channel section
details and other test specifications are provided in the Section 6.4.7. For the desired
Courant number, 0.8, it took 51 minutes and 16 seconds for CanaIFlow to execute the
simulation task for entire network with dx = 2.5 m. The simulation was run for 6
hours.
The initial and final flow profile in Channels 0-1 to 0-4 is shown in the Figure 7.22(A)
while 7.22(B) shows the discharge and depth hydrographs for 0-1 and C-4. Discharge
hydrograph and mass balance at junctions J-1 to J-3 is presented in the Figure 7.23.
Percentage error in discharge in C-1 and sum of discharges in C-2 and C-5 is shown
in the Figure 7.23(A). The error range is between -0.003 % and 0.28 %. The junction
clearly models the gate opening behaviour in the first hour and then closing in the
third hour of the simulation. Similar analysis is carried out for junctions J-2 and J-3.
The mass balance error for J-2 is in between -0.013 % and 0.13 % while for J-3 it is in
between -0.02 % and 0.18 %.
206
7.10 C.3: Irrigation Network
5.0
4.0
~IHr
]: Initial 6Hr
c
~ 3.0
•
>
•
iil
Bed
2.0
1.0
1 501 1001 1501 2001 2501 3001
Chalnage (m)
(A)
30 5
'U 24
I~
• ~
!
E
18
/'
•0>
I~ ~
!1i 12
~
~
c 6
/
o o
15 60 105 150 195 240 285 330
Time (mln)
(B)
Figure 7.22: Irrigation network (C.3) test case: (A) Initial flow profile and a profile
after 1 and 6 hour for the Channels C-l to C-4, (B) Discharge hydrograph plotted
on primary y-axis and depth hydrograph on secondary y-axis for the upstream side
of Channel C-l and downstream side of the Channel C-4. The behaviour seen in the
discharge hydrograph is obvious because of the increase in the gate opening during the
first hour and then it was kept constant during second hour then reduced in third hour
and then kept constant for rest of the period.
207
7.10 C.3: Irrigation Network
(A)
I-e-C2 _ _ C3 -e-C20 ___ Error I
(B)
i 15 0.3
I't~;;!::::::::::::~E f
15 90 165
Time (mln)
240 315
(C)
Figure 7.23: Irrigation network (C.3) test case: Discharge hydrograph and mass balance
at junctions at (A) J·l where channel C-l is split into C-2 and C-5 with a rectangular
gate on C-5 (E) J-2 where channel C-2 is split into C-3 and C-20 with a rectangular gate
on C-20 and (C) J-3 where channel C-3 is split into C-4 and C-33 with a rectangular
gate on C-33. The percentage error at each junction is a error between discharge of the
upstream channel and sum of the discharges of split channels.
208
7.11 R.l: Dam-Break in a Rectangular Channel and in a Channel with
Contraction and Expansion
As described in the Section 6.4.8, this test is an assessment of the ability of CanalFlow
to replicate the behaviour of a surge wave, caused by the sudden collapse of a large
body of water, in a straight rectangular channel (Gottardi and Venutelli, 2004) and in
a channel with a local constriction and expansion and compare the numerical results
against laboratory results (Crowder et al., 2004h).
The test setup is as shown in the Figure 6.9, page 162. The test data used is presented
in the Table 13.11, Appendix B. No special program to develop initial conditions was
used as opposed to various softwares used in Crowder et al. (2004h). As given in the
initial data table, cross-sections were defined at 0.1 m and time step used for this
test was 0.01 seconds and output time step was 0.04 seconds. It took 14 seconds for
CanalFlow to execute this test. This test case was run using Minmod, Monotonised
Centered, Superbee, van Leer and no limiter option. The results are analysed in terrns
of depth time series at locations SI, S2, S3 and S4 (see Figure 6.9).
For the first part of the test, the analysis is limited to the stage versus distance at
t = 6 seconds. The space domain is discretised in 400 cells. The figure 7.24 shows good
agreement between CanalFlow and analytical solution. The analytical solution for this
part of the test is given in Stoker (1957).
For the second part of this test case, the analysis is limited to the stage versus time
comparison at SI, S2, S3 and S4. The water level time series at these locations is shown
in Figures 7.26 to 7.28. As seen in all figures, the water level predicted by CanalFlow at
SI, starts dropping at approximately 0.32 seconds as opposed to experimental set-up
at 0.40 seconds. The drop in water level in CanalFlow is more prominent after 6.2
209
7.11 R.1: Dam-Break in a Rectangular Channel and in a Channel with
Contraction and Expansion
I-Analytical -CanalFlow I
12.0
10.0
8.0
I
~
c
.. 6.0
4.0
2.0
0.0
0 20 40 60 80 100 120 140 160 180 200
Chalnago (m)
Figure 7.24: Dam-Break in a Rectangular Channel: The figure shows depth along the
length of the rectangular channel with unit width at 6 seconds.
seconds. This is because, in the experimental set-up, sluice gate is used which opens
fully in 2 seconds while in CanalFlow sluice gate is not used (instant dam-break).
The calculated water level histories at the remaining three benchmarking stations
82, 83, and 84, show different behaviour in the change in water level to that of the
experimental data. It can be seen from the figures that the experimental water level
recorded at station 82 changes from 0.003 m after 3.0 seconds of the simulation. This
change in water level corresponds to the arrival of the surge wavefront at station 82.
Therefore, it can be determined that the surge wave has travelled 6.1 m in 3.0 seconds,
giving rise to an approximate travelling velocity of 2.0 m/sec. It can then be seen from
the plot of the experimental water level history that the water level only rises by a very
small amount between 3.5 seconds and 7.3 seconds, after which time the water level
increases rapidly, settling at an approximate level of 0.23 m after 8.5 seconds. The level
remains approximately constant thereafter for the remainder of the simulation. This
recorded rise in the water level after 7.3 seconds is due to the superimposition of the
resultant reflected wave caused by the angled walls at the front of the constriction on the
210
7.11 R.1: Darn-Break in a Rectangular Channel and in a Channel with
Contraction and Expansion
travelling surge wave. CanalFlow water level histories from Figures 7.25 to 7.28 show
that it does not follow this same process of change as that observed experimentally. The
water levels rise gradually after about 3.24 seconds and settle at approximately 0.09
m. However, in no limiter case, Figure 7.29(B), the behaviour at S2 is different than
other limiters. In this case, CanaIFlow follows the experimental trend where a sharp
jump has been seen after 7.3 seconds. The wave rises to approximately 0.24 m and then
settles at about 0.23 after 8.5 seconds. The discrepancies between the results in limiter
and non-limiter case at S2 are likely due to the one-dimensional solution generated by
the software. The test is not entirely appropriate for one-dimensional analysis as the
water surface becomes very irregular when the surge wave hits the constriction and it
would be appropriate to compare the CanalFlow results with a 2-D software package.
However, CanalFlow results are much better than ISIS, HEC-RAS and Mike 11 and it
did not create any instabilities in the solution.
It can be seen from the figures presenting the results at station S3, that the
CanalFlow calculated results fiuctuate around the experimental data after 4.68 sec-
onds. The final water level produced by CanalFlow at station S4 is almost same as
produced in the laboratory test. Again, the no limiter case is an exception where
CanalFlow result do not follow the experimental trend. The surge wave produced in
the physical model can be observed to arrive at station S2, S3 and S4 at 3.0, 4.32 and
5.52 seconds respectively. The same wave in CanalFlow arrives at 2.96, 4.48 and 5.2
seconds respectively.
211
7.11 R.I: Dam-Break in a Rectangular Channel and in a Channel with
Contraction and Expansion
1-Experimental -CanalFlow I
j
I-Experimental -CanalFlow I
0.4
E 0.3
i ~:+--_=~~ g ____
o 2 4 6 8
-I
10
Time (secs) Tlme (secs)
(A) (B)
~:: - - - - - - - - - - - - - 1
0.4
~03
t===:;;;;;;;;:;:.;~
.. 0.21-.1
0. 0.2
.!i f-----r
0.1 _ _ _ _ _..JL_ _-
0.0
1
_ _--1
c
O.'+--=::::::O~~
0.0..1- c
o 2 4 6 8 10 o 2 4 6 8 10
Time (sGca) TIme (_.ca)
(C) (D)
Figure 7.25: Dam break in Contractions and expansions test case (R.1): CanalFlow
predicted water depth time series using Minmod limiter at (A) 81, (B) 82, (C) 83 and
(D) 84.
212
7.11 R.l: Dam-Break in a Rectangular Channel and in a Channel with
Contraction and Expansion
I~:t-I=::::tt=~-:=~g:
0'4~
E 0.3 +"'~------------I
!:~ -
0.0 +----~_-_--_-_-_l
o 2 4 6 6 10
Tlme (sees)
2
•
Time (secs)
• "
(A) (B)
0.4
~0.3
0. 0.2 t.---------------1
t! ~:: t=====:::£O~A;¥;;'";.;.,,;
1 . ;;;.;;~
o •
Time (secs)
"
o 2 4 6
Time (secs)
8 10
(C) (D)
Figure 7.26: Dam break in Contractions and expansions test case (R.1): CanalFlow
predicted depth time series using Monotonised Centered limiter at (A) 81, (B) 82, (C)
83 and (D) 84.
213
7.11 R.l: Dam-Break in a Rectangular Channel and in a Channel with
Contraction and Expansion
TImlt(IIIICS)
IEI
c 0.1
0.0
0 2
r
, ,
TIme (secs)
a
, 10
(A) (B)
I
E0.3
iID;
0·'1
0.21-.-------------1 100·'.211-.- - - - - - - - - - - - - j
Q 0.1 +---------1:"""'------1
C:'
O.o.i-_ _ _..AL--_-_--l Q 0.1 t===::c;;;;;;;;;;;;;;;;;;;!
c .; ,
o 2 , 10 , , 0.0
o 2 4 6 10
Time (secs) Time (secs)
(C) (D)
Figure 7.27: Dam break in Contractions and expansions test case (R.l): Cana1Flow
predicted water depth time series using 8uperbee limiter at (A) 81, (B) 82, (C) 83 and
(D) 84.
214
7.11 R.l: Dam-Break in a Rectangular Channel and in a Channel with
Contraction and Expansion
o 2 4 , • 10
Time (SDca) Tlma (secs)
(A) (B)
~:: t-=====;;;;;;;;;;;;;;;;;~
iC01: F:
I::0.11-.
1--------------1
~:::C;;;;;~
0.0 l1 -_ _ _... '--~ _ _- - l -l-- c
0.0 1
o 2 4 6
Tlma (sacs)
• 10 o 2 4
Time (sacs)
,
• 10
(C) (D)
Figure 7.28: Dam break in Contractions and expansions test case (R.l): CanaIFlow
predicted water depth time series using van Leer limiter at (A) 81, (B) 82, (C) 83 and
(D) 84.
215
7.11 R.1: Dam-Break in a Rectangular Channel and in a Channel with
Contraction and Expansion
TIme (secs)
i:~
!!: "I
0.1
0.0
0 2
•
4 6
Tlme (secs)
d
8 10
(A) (B)
0.4 0.4
i:: t------,---------j
-----t i:: 1-------------1
.!i }------,C-'",
0.1
1
O.o.l-_
..
_ _...J:-._ _ _- - - I ~ 0.10.0 t======::;;;;;;;;;;;;;;;;;;;;;;;J
1
o 2 4 6 8 10 o 2 4 6 8 10
Time (sees) Time (secs)
(C) (D)
Figure 7.29: Dam break in Contractions and expansions test case (R.1): CanalFlow
predicted water depth time series using no limiter at (A) 81, (B) 82, (C) 83 and (D)
84.
216
CHAPTER 8
8.1 Conclusions
The canal hydraulic simulation softwares can be effectively used to increase efficiencies
of the irrigation schemes. The importance of such softwares in managing canal irrigation
flows is evident from the papers published in a special issue of Journal of Irrigation and
Drainage Engineering (119(4), July/August, 1993). As discussed in the Section 2.3,
there are certain problems and challenges which are unique to modelling canal irrigation
networks. Yet the technology adoption in this field has not kept a pace with advances
in computing and software. It was decided to review existing hydraulic simulation
softwares before deciding on the requirements of the new software. We looked at 19
existing hydraulic simulation softwares and their underlying numerical algorithm and
technology platform. The summary of these softwares is presented in the Table 2.1 on
page 45. Out of 20 softwares listed (including CanaIFlow), 75% of them fall under the
broad category of implicit finite difference schemes and within that broad category 80%
of them fall under four-point implicit Preissmann scheme. Meselhe and Holly (1997)
have clearly pointed out the invalidity of the Preissmann scheme in simulating trans-
critical flow regime which commonly occurs in canal irrigation networks. In recent
years, the use of finite volume method has become common in solving shallow water
equations. It has been widely applied to solve various problems in the field of hydraulic
217
8.1 Conclusions
engineering (see Section 2.4.4) but no known reference has been found where a finite
volume method with high-resolution techniques has been applied to develop a canal
hydraulic simulation software.
On the other hand, the work done by Kutija (1998); Kutija and Murray (2002);
Murray (2003) has stressed the need of implementing object-oriented paradigm in de-
veloping numerical algorithms. As such OOP is not new and it has been widely used
in developing GUI components around the existing numerical algorithms (Kutija and
Murray, 2002). In 2002, when this research work was started we looked at all available
technology options. In doing so we decided to set a new standard in the development of
hydraulic simulation softwares. One major criterion in selecting the technology options
was the possible migration of the new software to a web-based one with a multi-user
support. Thus, it was decided to use latest techniques in numerical algorithms and
latest programming tools available in software development.
The principal aim of this research was to develop a database-driven, object-oriented
hydraulic simulation model for canal irrigation networks using a finite volume method
with high-resolution shock capturing techniques. On the technology front, we decided
primarily to use Java, MySQL and CVS. The reasons behind selecting these technologies
are provided in the Section 3.3.2.1, for Java and in the Section 3.4.1, for MySQL. The
role of CVS is described in the Section 3.5. On the numerical side, it was decided to
use FVM with high-resolution shock capturing techniques (see Section 2.4.4) and an
exact Rlemann solver as described in the Section 4.4.2 was used.
It is evident from the Chapter 4 and 5 that we have achieved the above-mentioned
principal aim set for ourselves in developing a new canal hydraulic simulation software.
As the output of this research, a new category of hydraulic simulation software was
added (Category 6 in the Table 2.1, page 45). The CanalFlow model is integrated with
a relational database management system and in this thesis we present the data model
(Section 5.5) for canal hydraulic simulation software. The metadata is presented in the
Appendix A. This data model can be easily converted to XML schema which is fast
becoming a de facto standard in a data exchange.
218
8.1 Conclusions
The next obvious task was to test the model. As described in the Chapter 6, a
benchmark test suite developed by DEFRA/EA for 1-D river models was used (Crowder
et al., 2004a). As specified in Crowder et al. (2004i), the tests were classified under
the category of numerical accuracy, software capability and reproducibility (see Table
6.1, page 142). In all, 7 tests relevant to canal irrigation network were shortlisted from
theDEFRA/EA benchmark suite. A test case which models canal irrigation network
(Section 6.4.7) was added under the software capability category adapted from Misra
et al. (1992).
The CanalFlow results for the test suite are presented in the Chapter 7. We discuss
here conclusions of these tests:
• Ippen Wave: Discrepancies were found in the Ippen wave test case (Section
6.4.3). From the Section 7.6.2, it is seen that all calculated water levels
prior to any superimposition closely follow the analytical values. However,
as soon as the two waves meet, the calculated water levels show a varying
degree of divergence from the analytical values. Nevertheless, to a reasonable
219
8.1 Conclusions
degree of accuracy all the models predict the times at which the reflected
wave falls into phase with the incident wave. CanalFlow and other softwares
in the DEFRA/EA benchmark suite over-estimate the peak water level at
the closed boundary, which in part may be a result of the wave length being
slightly longer than the channel length, 800 m (Crowder et al., 2004f). Better
results could have been achieved by running this test case by reducing the
time step and running the model in a fixed time-step mode. In the current
work, this type of analysis in not carried out. In future, this test could be
run using different limiters to study the effect on water level and velocity
profiles and time series.
• Monoclinal Wave: CanaIFlow was able to model the monoclinal wave (Test
N.4 in Section 6.4.4) to a degree of accuracy that would be acceptable in
most practical situations. Some variations in predicting water levels were
seen as described in the Section 7.7.2. The CanalFlow results were compared
with the MIKE 11 results as ISIS and HEC-RAS were unable to model this
test case with Chezy's friction coefficient (the test specification stipulates
the use of Chezy's relation). It is evident from the Table 7.5 that error is
least when CanalFlow is run using van Leer limiter. In this case, CanalFlow
results are comparable to MIKE 11 results.
• Looped System: The Section 7.8.2 shows that CanalFlow is capable of han-
dling looped network. A minor adjustment of splitting Channel B in Bl
and B2 was done to configure this test case in CanalFlow. From the Table
7.6 and 7.6, it is evident that CanalFlow has shown some discrepancies in
mass balance at the junctions. The range of variation is between 0.02 % to
0.03 % for QSl and between -0.008 % to 0.02 % for QS2. This needs to be
addressed in the future versions.
• Weir: The Section 7.8.2 shows that CanalFlow is capable of modelling a weir
in the network. The results for steady state free flow condition match with
220
8.1 Conclusions
that of ISIS (Table 7.8) while the results for submerged flow condition are
different than other softwares but closer to ISIS .
• Irrigation Network: The results for canal irrigation network test (Section
6.4.7) are presented in the Section 7.10.2. It is evident from the results that
CanalFlow is capable of handling medium size canal irrigation networks with
a reservoir and control structures. It is also shown that the gates can be
operated as per the operation schedule. In the current version, CanalFlow
is capable of handling only rectangular gates. In future, other gates (for
example circular, radial) could be added. It would be interesting to integrate
the canal control algorithm and gate optimisation algorithm in the existing
code. Also a test case needs to be devised to simulate shock fronts in the
network.
221
8.2 Recommendations for the Future Work
using different limiters. It is noted here that CanalFlow did not require
the special setup of initial conditions as described in Crowder et al. (2004h)
to configure !SIS. HEC·RAS was unable to produce any results for this
test case and only MIKE 11 produced a stable solution for the dam-break
progressive wave but the results calculated were not truly representative of
the problem (Crowder et al., 2004h, p.23-25).
These tests show that CanalFlow is capable of handling variety of fiow scenarios,
viz., sub-critical, super-critical and trans-critical flow, flow in triangular channels, tidal
wave propagation, uniformly progressive wave, flow in a looped and dendritic network,
flow through weirs and gated structures and modelling a surge wave. It has outper-
formed ISIS, HEC-RAS and MIKE 11 in three of the benchmark tests and equally well
for the other four. The developed database schema allows great flexibility in inserting
and retrieving the data.
- More hydraulic structures can be added with modelling support for bridges,
culverts, aqueducts and siphons. Handling of structures in the existing code
needs to be improved, mainly, some more work is needed to take Care of
smooth transitions between free and submerged flow conditions at structures;
222
8.2 Recommendations for the Future Work
- Further test cases need to be run to simulate shock fronts in a large size
canal network;
- Verification of Contractions and Expansions test case with 2-D shallow water
model will be required to find out weaknesses in the l-D model;
- The existing model building process, which retrieves the data from the
RDBMS needs to be replaced with a service oriented layer to make it data
source independent. The constructors of network component classes (like
Network, Channel and ChannelStructure) will have to be refactored to ac-
commodate the database independent parameters;
- The existing Message class is a poor attempt to write the logging system.
This could be replaced by a better third party library called as log4j (http:
//logging.apache.org/log4j/docs/);
- jUnit test cases need to be added for all classes to improve the test coverage;
- The existing version control system CVS could be replaced with better sys-
tem called Subversion (http://subversion.tigris.org) .
223
Glossary
Abstraction Abstraction is a mechanism and practice to reduce and factor out details
so that one can focus on few concepts at a time, 54
Canal automation Canal automation refers to the actions performed without any
human intervention (closed-loop) in which a gate or pump changes its posi-
tion/setting in response to a water level, flow rate, or pressure because that
level/rate/pressure is different from the intended target value. The automation
may be performed through hydraulic, electrical, electronic, or a combination of
these means, 32, 46
Canal irrigation system The physical components of the irrigation facility, i. e., main
irrigation canal, distribution network and control structures, 1, 4
Canal pool A segment of canal bound between two controlling structures such as
cross-regulators, 9
Canal reach A segment of a canal along which the parameters like cross-sectional
shape, longitudinal slope and roughness coefficient are constant, 15, 31
224
GLOSSARY
Channel Also called as a reach is the segment of a canal along which the parame-
ters like cross-sectional shape, longitudinal slope and roughness coefficient are
constant, 85, 98, 107
Characteristics Characteristics are lines or curves in the (x, t) plane along which
disturbances propagate, 19
Checkout A checkout operation copies a working copy of a file from the repository to
the local disk, 65
Commit A commit occurs when you copy the changes you made on the local files
to the repository (the version control software takes care of knowing which files
changed since the last time the two were synchronised), 65
CPU A microprocessor chip that does most of the data processing; the CPU and
the memory form the central part of a computer to which the peripherals are
attached, 23
225
GLOSSARY
Distributary Canal taking off from a main canal (or sometimes from a branch canal),
usually supplying water to tertiary Iminor canals or 'directly to field offtakes, 11,
12
Dynamically-typed language The dynamically typed languages treat all data loca-
tions interchangeably, so inappropriate operations, like adding names or sorting
numbers alphabetically, will not cause errors until run-time, 53
Field application system The network of a field-level open channels in the area
below the outlet from a canal irrigation system, 9
Hydraulic jump A hydraulic jump represents a wave that desires to move in the
upstream direction, but which it is unable to do, as the flow velocity of the water
is larger than the wave's celerity (the speed of individual waves), 15
Irrigation The artificial application of water to land to promote the growth of crops,
1
Main canal A primary canal conveying water from the main intake of an irrigation
system to branch canals, 9
226
GLOSSARY
Metadata Metadata means data about data, it is information that describes another
set of data, 166
Performance ratio Ratio of the observed discharge at a specific location to the sched-
uled discharge at that location, 12
Repository The repository is where the files are stored, often on a server, 65
Riemann solvers The methods in which a numerical solution is built up from a series
of elementary Riemann problems at cell interfaces are called Riemann solvers, 24,
25
227
GLOSSARY
Rotational irrigation system This type of system involves dividing the command
area into two or more sections and diverting all available water to each section in
turns. The duration of a turn can be fixed or variable depending on the rotation
schedule, 6
Scalable It is the ease with which a system or component can be modified to fit the
problem area, 3
Trans-critical The term trans-critical flow denotes the existence of both sub-critical
and super-critical flow regimes in an open-channel system. In irrigation networks
this flow situation is caused because of dynamic gate operation, 20
Version control It is the management of multiple revisions of the same unit of in-
formation. It is most commonly used in engineering and software development
to manage ongoing evolution of digital documents like application source code
or electronic models and other critical information that may be worked on by a
team of people, 5, 64
228
APPENDIX A
Database Design
229
APPENDIX A. DATABASE DESIGN
230
APPENDIX A. DATABASE DESIGN
231
APPENDIX A. DATABASE DESIGN
232
APPENDIX A. DATABASE DESIGN
233
APPENDIX A. DATABASE DESIGN
234
APPENDIX A. DATABASE DESIGN
235
APPENDIX A. DATABASE DESIGN
236
APPENDIX A. DATABASE DESIGN
237
APPENDIX B
Test Dataset
Table B.I: Initial data for the sub-critical (N.l.I) sub-part in the sub-critical, super-
critical and transitional flows test case. Chainage, bed, depth are given in m while
discharge is in cumec.
238
APPENDIX B. TEST DATASET
239
APPENDIX B. TEST DATASET
Table B.2: Initial data for the super-critical flow (N.1.2) sub-part in the sub-critical,
super-critical and transitional flows test case. Chainage, bed, depth are given in m
while discharge is in cumec.
240
APPENDIX B. TEST DATASET
241
APPENDIX B. TEST DATASET
Table B.3: Initial data for the super-critical to sub-critical flow (N.1.3) sub-part in the
sub-critical, super-critical and transitional flows test case. Chainage, bed, depth are
given in m while discharge is in curnec.
242
APPENDIX B. TEST DATASET
243
APPENDIX B. TEST DATASET
Table B.4: Initial data for the sub-critical to super-critical to sub-critical flow (N.l.4)
sub-part in the sub-critical, super-critical and transitional flows test case. Chainage,
bed, depth are given in rn while discharge is in curnec.
244
APPENDIX B. TEST DATASET
245
APPENDIX B. TEST DATASET
Table B.5: Initial data for the super-critical to sub-critical to super-critical (N.1.5) flow
sub-part in the sub-critical, super-critical and transitional flows test case. Chainage,
bed, depth are given in m while discharge is in cumec.
246
APPENDIX B. TEST DATASET
247
APPENDIX B. TEST DATASET
Table B.6: Initial data for the sub-critical and super-critical flow test in triangular
channels (N.2.l and N.2.2). Chainage, bed, depth are given in m while discharge is in
cumec.
Test N.2.1 Test N.2.2
Chainage Bed Depth Discharge Chainage Bed Depth Discharge
0 13 3 20 0 13 1.7 20
300 12.7 3 20 15 12.7 1.7 20
600 12.4 3 20 30 12.4 1.7 20
900 12.1 3 20 45 12.1 1.7 20
1200 11.8 3 20 60 11.8 1.7 20
1500 11.5 3 20 75 11.5 1.7 20
1800 11.2 3 20 90 11.2 1.7 20
2100 10.9 3 20 105 10.9 1.7 20
2400 10.6
, 3 20 120 10.6 1.7 20
2700 10.3 3 20 135 10.3 1.7 20
3000 10 3 20 150 10 1.7 20
248
APPENDIX B. TEST DATASET
Table B.7: Initial data for the Ippen wave test (N.3). Chainage, bed, depth are given
in m while discharge is in cumec.
249
APPENDIX B. TEST DATASET
250
APPENDIX B. TEST DATASET
251
APPENDIX B. TEST DATASET
Table B.8: Initial data for Monoclinal rising wave test (N.4). Chainage, bed, depth are
given in ID while discharge is in cumec.
252
APPENDIX B. TEST DATASET
253
APPENDIX B. TEST DATASET
254
APPENDIX B. TEST DATASET
Table B.g: Initial data for the Looped System test(C.l). dx (distance step), bed, depth,
bottom are in m, discharge is in cumec and side slope is expressed as H:V.
255
APPENDIX B. TEST DATASET
Table B.IO: Channel geometry details used in the irrigation network (C.3) test case.
The length, bottom and depth are in m, slope is expressed as 1 in column value, side
slope is specified as horizontal to vertical while the discharge is in cumec.
256
APPENDIX B. TEST DATASET
257
APPENDIX B. TEST DATASET
Table B.Il: Initial data for test R.1 - Dambreak in converging and diverging channel
section. The chainage, bed, depth and bottom are given in m and discharge is in cumec.
258
APPENDIX B. TEST DATASET
259
APPENDIX B. TEST DATASET
260
APPENDIX B. TEST DATASET
261
APPENDIX B. TEST DATASET
262
APPENDIX B. TEST DATASET
263
References
ABBOTT, M.B. (1975). Unsteady Flow in Open Channels, K. Mohmmod and V. Yev-
jevich, Eds, vol. 1, chap. Method of Characteristics, 63-88. Water Resources Publi-
cations, Fort Collins, Colorado, USA. 2.4.1
ABBOTT, M.B. AND BASCO, D.R. (1989). Computational Fluid Dynamics An Intro-
duction for Engineers. Longman Group, UK. 1.1, 2.4.2, 2.6.14
264
REFERENCES
ANDERSON, JR., J.D. (1995). Computational Fluid Dynamics The Basics with Appli-
cations. McGraw-Hill, Inc. 2.4
265
REFERENCES
BERMUDEZ, A., DERVIEux, A., DESIDERI, J.A. AND VAZQUEZ, M.E. (1998). Upwind
Schemes for the Two-Dimensional Shallow Water Equations with Variable Depth us-
ing Unstructured Meshes. Computer Methods in Applied Mechanics and Engineering,
155,49-72. 2.4.4
BOISVERT, R.F., DONGARRA, J.J., Pozo, R., REMINGTON, K.A. AND STEWART,
G.W. (1998). Developing Numerical Libraries in Java. Concurrency: Practice and
Experience, 10, 1117-1129. 3.3.2.2, 3.3.2.2
BOISVERT, R.F., MORElRA, J., PHILIPPSEN, M. AND Pozo, R. (2001). Java and Nu-
merical Computing. Computing in Science and Engineering, IEEE, 3, 18-24. 3.3.2.2
BOOCH, G. (1994). Object-Oriented Analysis and Design with Applications. The Ben-
jamin Cummings Publishing Company. (document), 3.2, 3.1, 3.3, 5.1, 5.2.1
BORTHWICK, A.G.L., LEON, S.C. AND JOZSA, J. (2001). Adaptive Quadtree Model
of Shallow-Flow Hydrodynamics. Journal of Hydraulic Research, IAHR, 39, 413-424.
2.7
BRODIE, R.S. (1998). Integrating GIS and RDBMS Technologies During Construction
of a Regional Groundwater Model. Environmental Modelling and Software, 14, 119-
128. 3.4
BRUNNER, G.W. (2001a). HEC-RAS, River Analysis System Hydraulic Reference Man-
ual. Hydrologic Engineering Center, U. S. Army Corps of Engineers, California, USA.
2.6.9
266
REFERENCES
BRUNNER, G.W. (200lb). HEC-RAS, River Analysis System Users Manual. Hydro-
logic Engineering Center, U. S. Army Corps of Engineers, California, USA. 2.6.9
BURT, C.M. AND STYLES, S.W. (1999). Modern Water Control and Management
Practices in Irrigation: Impact on Performance. Food and Agricultural Organisation
Water Report Number 19. 1
CALEFFI, V., VALIANI, A. AND ZANNI, A. (2003). Finite Volume Method for Simu-
lating Extreme Flood Events in Natural Channels. Journal of Hydraulic Research,
IAHR, 41, 167-177. 2.4.4
CAPART, H., ELDHO, T.!., HUANG, S.Y., YOUNG, D.L. AND ZECH, Y. (2003).
Treatment of Natural Geometry in Finite Volume River Flow Computations. Journal
of Hydraulic Engineering, ASCE, 129, 385-393. 2.4.4
CARLETON, C.J., DAHLGREN, R.A. AND TATE, K.W. (2005a). A Relational Database
for the Monitoring and Analysis of Watershed Hydrologic Functions: I. Database
Design and Pertinent Queries. Computers and Geosciences, 31, 393-402. 3.4
CARLETON, C.J., DAHLGREN, R.A. AND TATE, K.W. (2005b). A Relational Data-
base for the Monitoring and Analysis of Watershed Hydrologic Functions: n. Data
Manipulation and Retrieval Programs. Computers and Geosciences, 31, 403-413. 3.4
267
REFERENCES
CASANOVA, H., DONGARRA, J. AND DOOLIN, D.M. (1997). Java Access to Numerical
Libraries. Concurrency: Practice and Experience, 9, 1279-1291. 3.3.2.2
CAUSON, D.M., MINGHAM, C.G. AND INGRAM, D.M. (1999). Advances in Calcu-
lation Methods for Supercritical Flow in Spillway Channels. Journal of Hydraulic
Engineering, ASCE, 125, 1039-1050. 2.4.4,2.7
CHIPPADA, S., DAWSON, C.N., MARTINEZ, M.L. AND WHEELER, M.F. (1998). A
Godunov-Type Finite Volume Method for the System of Shallow Water Equations.
Computer Methods in Applied Mechanics and Engineering, 151, 105-129. 2.4.4, 2.7
CHL (1997). HIVEL 2D v2.0 Users Manual. Coastal and Hydraulics Laboratory
(CHL), Waterways Experiment Station, Vicksburg, MS, USA. 2.6.10
CLEMMENS, A.J. AND REPLOGLE, J .A. (1989). Control of Irrigation Canal Networks.
Journal of Irrigation and Drainage Engineering, ASCE, 115,96-110. 2.2
CLEMMENS, A.J., HOLLY, F.M. AND SCHUURMANS, W. (1993). Description and Eval-
uation of Program: Duflow. Journal of Irrigation and Drainage Engineering, ASCE,
119, 724-734. 2.4.2, 2.6.5, 2.7
CLEMMENS, A.J., BAUTISTA, E., WAHLIN, B.T. AND STRAND, R.J. (2005). Sim-
ulation of Automatic Canal Control Systems. Journal of Irrigation and Drainage
Engineering, ASCE, 131, 324-335. 2.6.2
268
REFERENCES
CODD, E.F. (1970). A Relational Model of Data for Large Shared Data Banks. Com-
munications of the ACM, 13,377-387.3.4
CONTRACTOR, D.N. AND SCHUURMANS, W. (1993). Informed Use and Potential Pit-
falls of Canal Models. Journal of Irrigation and Drainage Engineering, ASCE, 119,
663-672. 2.7
CROWDER, R.A., PEPPER, A.T., WHITLOW, C., SLEIGH, A., WRIGHT, N. AND
TOMLIN, C. (2004a). Benchmarking of Hydraulic River Modelling Software Packages:
Project Overview. DEFRA/Environment Agency Flood and Coastal Defence R & D
Programme. R & D Technical Report: W5-105/TRO. 2.7,6.3,8.1
CROWDER, R.A., PEPPER, A.T., WHITLOW, C., SLEIGH, A., WRIGHT, N. AND
TOMLIN, C. (2004b). Benchmarking of Hydraulic River Modelling Software Pack-
ages: Results Test A (Subcritical, Supercritical and Thansitional Flows). DE-
FRA/Environment Agency Flood and Coastal Defence R & D Programme. R &
D Technical Report: W5-105/TR2A. 6.4.1.2, 2, 7.4.1, 3, 4
CROWDER, R.A., PEPPER, A.T., WHITLOW, C., SLEIGH, A., WRlGHT, N. AND
TOMLIN, C. (2004c). Benchmarking of Hydraulic River Modelling Software Packages:
Results Test B (Looped Systems). DEFRA/Environment Agency Flood and Coastal
Defence R & D Programme. R & D Technical Report: W5-105/TR2B. (document),
6.4.5,7.8, 7.6, 7.7
CROWD ER, R.A., PEPPER, A.T., WHITLOW, C., SLEIGH, A., WRlGHT, N. AND
TOMLIN, C. (2004d). Benchmarking of Hydraulic River Modelling Software Pack-
ages: Results Test C (Thiangular Channels). DEFRA/Environment Agency Flood
and Coastal Defence R & D Programme. R & D Technical Report: W5-105/TR2C.
(document), 6.4.2, 7.2
269
REFERENCES
CROWDER, R.A., PEPPER, A.T., WHITLOW, C., SLEIGH, A., WRlGHT, N. AND
TOMLIN, C. (2004e). Benchmarking of Hydraulic River Modelling Software Packages:
Results Test D (Weirs). DEFRA/Environment Agency Flood and Coastal Defence
R & D Programme. R & D Technical Report: W5-105/TR2D. 6.4.6,7.9
CROWDER, R.A., PEPPER, A.T., WHITLOW, C., SLEIGH, A., WRlGHT, N. AND
TOMLIN, C. (2004f). Benchmarking of Hydraulic River Modelling Software Packages:
Results Test E (Ippen Wave). DEFRA/Environment Agency Flood and Coastal De-
fence R & D Programme. R & D Technical Report: W5-105/TR2E. 6.4.3,7.6.1,1,
2, 1
CROWDER, R.A., PEPPER, A.T., WHITLOW, C., SLEIGH, A., WRlGHT, N. AND
TOMLIN, C. (2004g). Benchmarking of Hydraulic River Modelling Software Packages:
Results Test F (Monoclinal Wave). DEFRA/Environment Agency Flood and Coastal
Defence R & D Programme. R & D Technical Report: W5-105/TR2F. 6.4.4, 7.7,
7.7.1
CROWDER, R.A., PEPPER, A.T., WHITLOW, C., SLEIGH, A., WRlGHT, N. AND
TOMLIN, C. (2004h). Benchmarking of Hydraulic River Modelling Software Packages:
Results Test L (Contraction and Expansion). DEFRA/Environment Agency Flood
and Coastal Defence R & D Programme. R & D Technical Report: W1>-105/TR2L.
6.4.8, 6.4.8.2, 7.11, 7.11.1, 3
CROWDER, R.A., PEPPER, A.T., WHITLOW, C., SLEIGH, A., WRlGHT, N. AND
TOMLIN, C. (2004i). Benchmarking of Hydraulic River Modelling Software Packages:
Test Specifications. DEFRA/Environment Agency Flood and Coastal Defence R &
D Programme. R & D Technical Report: W5-105/TR1. 2.7,5.4.3, 6.3, 6.4.1, 6.4.1.2,
7.4,8.1
270
REFERENCES
DELlS, A.!, AND SKEELS, C.P. (1998). TVD Schemes for Open Channel Flow. Inter-
national Journal for Numerical Methods in Fluids, 26, 791-809. 2.4.2
DELlS, A.!., SKEELS, C.P. AND RYRIE, S.C. (2000). Implicit High-Resolution Meth-
ods for Modelling One-Dimensional Open Channel Flow. Journal of Hydraulic Re-
search, IAHR, 38, 369-382. 2.4.2, 2.4.2, 2.4.5, 2.7
DELONG, L.L., THOMPSON, D.B. AND LEE, J.K. (1997). The Computer Program
FourPT (Version 95.01) - A Model for Simulating One-Dimensional, Unsteady,
Open-Channel Flow. Water Resources Investigation Report 97-4016, U. S. Geological
Survey, Mississippi, USA. 2.4.2, 2.6.8
FIEDLER, F.R. AND RAMIREZ, J.A. (2000). A Numerical Method for Simulating Dis-
continuous Shallow Flow Over an Infiltrating Surface. International Journal for Nu-
merical Methods in Fluids, 32, 219-240. 2.4.2
FOGEL, K. AND BAR, M. (2001). Open Source Development with CVS Second Edition.
Coriolis Technology Press, USA. 3.5
Fox, G., LI, X., QIANG, Z. AND ZHIGANG, W. (1997). A Prototype Fortran-to-Java
Converter. Concurrency: Practice and Experience, 9, 1047-1061. 3.3.2.2
271
REFERENCES
FRANZ, D.D. AND MELCHING, C.S. (1997). Full Equations (FEQ) Model for the Solu-
tion of the Full, Dynamic Equations of Motion for One-Dimensional Unsteady Flow
in Open Channels and Through Control Structures. Water Resources Investigations
Report 96-4240, U. S. Geological Survey, California, USA. 2.4.2, 2.6.6
FROEHLICH, D.C. (2002). Finite Element Surface- Water Modeling System: Two-
Dimensional Flow in a Horizontal Plane Version 2, Draft Users Manual. U. S. De-
partment of Transportation, Federal Highway Administration, USA. 2.6.7
GARCIA-NAVARRO, M.P. AND SAVIRON, J.M. (1992a). McCormack's Method for the
Numerical Simulation of One-Dimensional Discontinuous Unsteady Open Channel
Flow. Journal of Hydraulic Research, IAHR, 30, 95-105.4.6.2.1,4.6.3.5,4.6.4,7.4.1
GICHUKI, F.N., WALKER, W.R. AND MERKLEY, G.P. (1990). Transient Hydraulic
Model for Simulating Canal-Network Operation. Journal of Irrigation and Drainage
Engineering, ASCE, 116, 67-81. 2.2
GOSLING, J., JOY, B., STEELE, G. AND BRACHA, G. (2005). The Java Language
Specification. Addison-Wesley. 3.3.2.2
272
REFERENCES
GRELLE, N., BERTRAND, G., FRAZAO, S., HIVER, J. AND ZECH, Y. (2003). Flows
Due to a Structure Failure on the Canal Du Centre. In Proceedings of xxx IAHR
Congress, Thessaloniki, Greece, Vol. D, 543-550. 2.2
HALCROW (1998). Irrigation Canal Control Vol. 1: User Manual. DFID Research
Project R6259. Halcrow Group Ltd., Swindon, UK. 2.2
HAUKE, G. (2002). A Stabilized Finite Element Method for the Saint-Venant Equa-
tions with Application to Irrigation. International Journal for Numerical Methods in
Fluids, 38, 963-984. 2.4.3
273
REFERENCES
HOLLY, F.M. AND MERKLEY, G.P. (1993). Unique Problems in Modelling Irrigation
Canals. Journal of Irrigation and Drainage Engineering, ASCE, 119,656-662. 2.3
HOLLY, F.M. AND PARRISH Ill, J.B. (1993). Description and Evaluation of Program:
CARIMA. Journal of Irrigation and Drainage Engineering, ASCE, 119, 703-713.
2.4.2, 2.6.2, 2.6.4, 2.7
HOLLY, F.M. AND PARRISH Ill, J.B. (n.d.). Numerical Simulation of Unsteady Hy-
drodynamics in Canals with (Optional) User-Supplied Gate Control Algorithms.
http://www . Ehr. uioya. edu/projects/canalcad/index.html. 2.6.2
Hsu, C.C., LEE, W.J. AND CHANG, C.H. (1998). Subcritical Open-Channel Junction
Flow. Journal of Hydraulic Engineering, ASCE, 124, 847-855. 4.6.4
Hu, K., MINGHAM, C.G. AND CAUSON, D.M. (1998). A Bore-Capturing Finite-
Volume Method for Open Channel Flows. International Journal for Numerical Meth-
ods in Fluids, 28, 1241-1261. 2.4.4,2.4.5,2.7
274
REFERENCES
Hu, Y.F., ALLAN, R.J. AND MAGUIRE, K.C.F. (2000). Comparing the Performance
of Java with Fortran and C for Numerical Computing. http://www.dl.ac . uk/TCSC/
UKHEC/JASPA/bench.html. 3.3.2.2
JOHNSON, R.E. AND FOOTE, B. (1988). Designing Reusable Classes. Journal of Object-
Oriented Programming, 1, 22-35. 3.3.1
JOLMA, A. (1998). Sharing the Information and the Tools: Web Suppport for Water
Resources Management. In Proceedings of the Third International Conference on
Hydroinformatics, Copenhagen, Denmark, 24-26 Aug., A.A. Balkema, Rotterdam,
The Netherlands, 935-940. 1.1.1
KARELSKY, K.V., PAPKOV, V.V. AND PETROSYAN, A.S. (2000). The Initial Discon-
tinuity Decay Problem for Shallow Water Equations on Slopes. Physics Letters A,
271, 349-357. 2.4.5
275
REFERENCES
KONG, X.A. AND CHEN, D.P. (1995). An Object-Oriented Design ofFEM Programs.
Computers and Structures, 57, 157-166. 3.3.1
KUMAR, P., MISHRA, A., RAGHUWANSHI, N.S. AND SINGH, R. (2002). Application
of Unsteady Flow Hydraulic-Model to a Large and Complex Irrigation System. Agri-
cultural and Water Management, 54, 49-66. 2.2, 2.7
276
REFERENCES
LAL, A.M.W., ZEE, R.V. AND BELNAP, M. (2005). Case Study: Model to Simulate
Regional Flow in South Florida. Journal of Hydraulic Engineering, ASCE, 131, 247-
258. 3.3.1
LEVEQUE, R.J. (2002). Finite Volume Methods for Hyperbolic Problems. Cambridge
University Press, UK. 1.1, 2.4.4, 2.4.5, 4.3, 4.3, 4.4, 4.4.2, 2, 6, 4.5, 4.5.1, 4..5.1
LIGGETT, J.A. AND CUNGE, J.A. (1975). Unsteady Flow in Open Channels, K. Mohm-
mod and V. Yevjevich, Eds, vo!. 1, chap. Numerical Methods of Solution of the
Unsteady Flow Equations, 89-182. Water Resources Publications, Fort Collins, Col-
orado, USA. 1.1,2.4.1, 4.0.2.1
LIN, G.F., LAI, J.S. AND Guo, W.D. (2003). Finite-Volume Component-Wise TVD
Schemes for 2D Shallow Water Equations. Advances in Water Resources, 26, 861-
873. 2.4.4
LIN, G.F., LAI, J.S. AND Guo, W.D. (2005). High-Resolution TVD Schemes in Finite
Volume Method for Hydraulic Shock Wave Modelling. Journal of Hydraulic Research,
IAHR, 43, 376-389. 2.4.4
LIU, J.L., LIN, I.J., SHIH, M.Z., CHEN, R.C. AND HSIEH, M.C. (1996). Object-
Oriented Programming of Adaptive Finite Element and Finite Volume Methods.
Applied Numerical Mathematics, 21, 439-467. 3.3.1
277
REFERENCES
MAcDoNALD, I. (1994). Test Problems with Analytic Solution for Steady Open Chan-
nel Flow. University of Reading, Numerical Analysis Report 6/94. 6.4.1.1, 7.4
MAcDoNALD, I., BAINES, M.J., NICHOLS, N.K. AND SAMUELS, P.G. (1997). Ana-
lytic Benchmark Solutions for Open-Channel Flows. Journal of Hydraulic Engineer-
ing, ASCE, 123, 1041-1045. 2.7
MALATERRE, P.O. AND BAUME, J.P. (1997). SIC 3.0, A Simulation Model for
Canal Automation Design. In International Workshop on the Regulation of Irri-
gation Canals: State of the Art of Research and Applications, RIG 97, Marrakech
(Morocco), April 22-24, 1997. 2.4.2,2.6.15,2.7
278
REFERENCES
MCKINNEY, D.C. AND CAI, X. (2002). Linking GIS and Water Resources Management
Models: An Object-Oriented Method. Environmental Modelling and Software, 17,
413-425. 3.3.1
MERKLEY, G.P. (1997). CanalMan: A Hydraulic Simulation Model for Unsteady Flow
in Branching Canal Networks. Department of Biological and Irrigation Engineering,
Utah State University, Utah, USA. 2.4.2, 2.6.3, 2.7
MERKLEY, G.P. AND ROGERS, D.C. (1993). Description and Evaluation of Program
Canal. Journal of Irrigation and Drainage Engineering, ASCE, 119, 714-723. 2.7
MESELHE, E.A. AND HOLLY, F.M. (1997). Invalidity ofPreissmann Scheme for Tran-
scritical Flow. Journal of Hydraulic Engineering, ASCE, 123,652-655. 2.4.2,8.1
MINGHAM, C.G. AND CAUSON, D.M. (2000). Calculation of Unsteady Bore Diffrac-
tion using a High-Resolution Finite-Volume Method. Journal of Hydraulic Research,
IAHR, 38, 49-56. 2.4.4, 2.4.5
MISHRA, A., ANAND, A., SINGH, R. AND RAGHUWANSHI, N.S. (2001). Hydraulic
Modeling of Kangsabati Main Canal for Performance Assessment. Journal of Irriga-
tion and Drainage Engineering, ASCE, 127, 27-34. 2.2, 2.2
MISRA, R. (1998). Recursive Algorithm for Steady Flow in a Canal Network. Advances
in Engineering Software, 29, 77-86. 4.6.4
MISRA, R., SRIDHARAN, K. AND KUMAR, M.S.M. (1992). Transients in Canal Net-
work. Journal of Irrigation and Drainage Engineering, ASCE, 118, 690-707. 2.2,
4.6.4,6.3,6.4.7,7.10,8.1
279
REFERENCES
MOREIRA, J.E., MIDKIFF, S.P., GUPTA, M., ARTIGAS, P.V., SNIR, M. AND
LAWRENCE, R.D. (2000). Java Programming for High-Performance Numerical Com-
puting. IBM Systems Journal, 39, 21-56. 3.3.2.2
MUTUA, B.M. AND MALANO, H.M. (2001). Analysis of Manual and Centralised Su-
pervisory Control Operations to Improve Level of Service: A Case Study of Pyramid
Hill No. 1 Channel, Victoria, Australia. Irrigation and Drainage Systems, 15, 1-19.
2.2
NAUGHTON, P. (1996). The Java Handbook. Osborne McGraw-Hill. McGraw Hill Inc.
3.3,3.3.2.1, 5.3.1
NOTo, L. AND TUCCIARELLI, T. (2001). DORA Algorithm for Network Flow Models
with Improved Stability and Convergence Properties. Journal of Hydraulic Engineer-
ing, ASCE, 127, 380-391. 4.6.4
280
REFERENCES
OMG (2003). Object Management Group. UML 2.0 Infrastructure Specification. http:
/ /www.omg.org. 5.2.1
PARRISH In, J.B. AND BURT, C.M. (1993). Cal Poly Model Canal. Journal of Irri-
gation and Dminage Engineering, ASCE, 119, 673-678. 2.7, 6.1
PERSONAL COMMUNICATION - CROWD ER (2005). Excel Data Sheet sent Over E-mail
by Richard Crowder on 13 April, 2005. 6.4.4.3, 7.5
PRICE, R.K. (1974). Comparison of Four Numerical Methods for Flood Routing. Jour-
nal of the Hydmulics Division, ASCE, 100, 879-899. 2.4.1,2.7
R (n.d.). The R Project for Statistical Computing. http://www .r-project. ~rg. 3.1
281
REFERENCES
REDDY, H.P. AND BHALLAMUDI, S.M. (2004). Gradually Varied Flow Computation in
Cyclic Looped Channel Networks. Journal of Irrigation and Drainage Engineering,
ASCE, 130,424-431. 4.6.4
RIBEIRO, F.L.B., GALEAO, A.C. AND LANDAU, L. (2001). Edge-Based Finite Element
Method for Shallow Water Equations. International Journal for Numerical Methods
in Fluids, 36, 659-685. 2.4.3
ROE, P .L. (1981). Approximate lliemann Solvers, Parameter Vectors and Difference
Schemes. Journal of Computational Physics, 43, 357-372. 2.4.2, 2.4.5
ROGERS, D.C. AND MERKLEY, G.P. (1993). Description and Evaluation of Program
USM. Journal of Irrigation and Drainage Engineering, ASCE, 119,693-702. 2.4.1,
2.6.19,2.7
ROGERS, D.C., KACEREK, T.F. AND GOOCH, R.S. (1993). Field Data for Verify-
ing Canal Unsteady Flow Models. Journal of Irrigation and Drainage Engineering,
ASCE, 119,679-692. 2.7, 6.1
Ross, T.J., WAGNER, L.R. AND LUGER, G.F. (1992a). Object-Oriented Program-
ming for Scientific Codes. I: Thoughts and Concepts. Journal of Computing in Civil
Engineering, 6, 480-495. 3.3.1
Ross, T.J., WAGNER, L.R. AND LUGER, G.F. (1992b). Object-Oriented Program-
ming for Scientific Codes. II: Examples in C++. Journal of Computing in Civil
Engineering, 6, 497-514. 3.3.1
282
REFERENCES
SANDERS, B.F. (2001). High-Resolution and Non-Oscillatory Solution of the St. Venant
Equations in Non-Rectangular and Non-Prismatic Channels. Journal of Hydraulic
Research, IAHR, 39, 321-330. 2.4.4, 2.4.5
SANDERS, B.F. (2002). Non-Reflecting Boundary Flux Function for Finite Volume
Shallow-Water Models. Advances in Water Resources, 25, 195-202. 2.4.5
SANDERS, B.F., GREEN, C.L., CHU, A.K. AND GRANT, S.B. (2001). Case Study:
Modeling Tidal Transport of Urban Runoff in Channels Using the Finite-Volume
Method. Journal of Hydraulic Engineering, ASCE, 127, 795-804. 2.4.4
SCHACH, S.R. (2004). Introduction to Object-Oriented Analysis and Design with UML
and the Unified Process. McGraw-Hill International Edition. 3.3, 5.1, 5.2.1
SEN, D.J. AND GARG, N.K. (2002). Efficient Algorithm for Gradually Varied Flows
in Channel Networks. Journal of Irrigation and Drainage Engineering, ASCE, 128,
351-357. 4.6.4
SHEN, J., PARKER, A. AND RIVERSON, J. (2005). A New Approach for a Windows-
based Watershed Modeling System based on a Database-supporting Architecture.
Environmental Modelling and Software, 20, 11271138. 3.4
283
REFERENCES
SHENDE, S., SMOUT, I. AND SCOTT, C. (2003). Hydraulic Simulation Model for Canal
Irrigation Systems Using Open Platform Technologies. First Year Progress Report
(Unpublished) Submitted to Loughborough University, UK. 5.3
SHEU, T.W.H. AND FANG, C.C. (2001). High Resolution Finite-Element Analysis of
Shallow Water Equations in Two Dimensions. Computer Methods in Applied Me-
chanics and Engineering, 190,2581-2601. 2.4.3,2.7
SLEIGH, P.A., GASKELL, P.H., BERZINS, M. AND WRIGHT, N.G. (1998). An Un-
structured Finite-Volume Algorithm for Predicting Flow in Rivers and Estuaries.
Computers and Fluids, 27, 479-508. 4.6.3.2
SMITH, T.H., GOWER, A.E. AND BONING, D.S. (1997). A Matrix Math Library for
Java. Concurrency: Practice and Experience, 9, 1127-1137. 3.3.2.2
SMOUT, I.K. AND GORANTIWAR, S.D. (2005). A Multilevel Approach for Optimizing
Land and Water Resources and Irrigation Deliveries for Tertiary Units in Large
Irrigation Schemes: 1. Method. Journal of Irrigation and Drainage Engineering,
ASCE, 131. 1.1.1
SOBEY, RJ. (2001). Evaluation of Numerical Models of Flood and Tide Propagation
in Channels. Journal of Hydraulic Engineering, ASCE, 127, 805-824. 4.6.4, 6.2
STOKER, J.J. (1957). Water Waves. Interscience Publishers, New York. 7.11.2
284
REFERENCES
STRELKOFF, T.S. AND FALVEY, H.T. (1993). Numerical Methods Used to Model
Unsteady Canal Flow. Journal of Irrigation and Drainage Engineering, ASCE, 119,
637-655. 2.5
SWAIN, E.D. AND CHIN, D.A. (1990). Model of Flow in Regulated Open-Channel
Networks. Journal of Irrigation and Drainage Engineering, ASCE, 116, 537-556.
4.6.4
TORO, E. (1997). Riemann Solvers and Numerical Methods for Fluid Dynamics. A
Practical Introduction. Springer-Verlag, Berlin, Germany. 4.3, 4.3, 4.3, 4.4.1
TORO, E. (2001). Shock-Capturing Methods for Free-Surface Shallow Flows. John Wi-
ley, USA. 2.4.5, 4.~1, 4.4
TSENG, M.H. (1999). Verification of 1-D 'Transcritical Flow Model in Channels. Pro-
ceedings of National Science Council, 23, 654-664. 2.4.2, 2.4.4
TSENG, M.H. AND CHU, C.R. (2000a). 2-Dimensional Shallow Water Flow Simulation
Using TVD-MacCormack Scheme. Journal of Hydraulic Research, 38, 123-131. 2.4.5,
2.7
TSENG, M.H. AND CHU, C.R. (2000b). The Simulation of Dam-Break Flows by an
Improved Predictor-Corrector TVD Scheme. Advances in Water Resources, 23, 637-
643. 2.4.2, 2.4.5
TSENG, M.H., Hsu, C.A. AND CHU, C.R. (2001). Channel Routing in Open-Channel
Flows with Surges. Journal of Hydraulic Engineering, ASCE, 127, 115-122. 2.4.4
285
REFERENCES
VALIANI, A., CALEFFI, V. AND ZANNI, A. (2002). Case Study: Malpasset Dam-Break
Simulation using a Two-Dimensional Finite Volume Method. Journal of Hydraulic
Engineering, ASCE, 128,460-472. 2.4.4, 2.4.5
VELICKOV, S., PRICE, R.K. AND SOLOMATINE, D.P. (1998). Internet for Hydroinfor-
matics - Practical Examples of Client /Server Modelling. In Proceedings of the Third
International Conference on Hydroinformatics, Copenhagen, Denmark, 24-26 Aug.,
A.A. Balkema, Rotterdam, The Netherlands. 1.1.1
WANG, J., NI, H. AND HE, Y. (2000). Finite Difference TVD Scheme for Computation
of Dam Break Problems. Journal of Hydraulic Engineering, ASCE, 126, 253-262.
2.4.2,2.7
WANG, J.W., HASSETT, J.M. AND ENDRENY, T.A. (2005). An Object-Oriented Ap-
proach to the Description and Simulation of Watershed Scale Hydrologic Processes.
Computers and Geosciences, 31, 425-435. 3.3.1
WIDENIUS, M., AXMARK, D. AND AB, M. (2002). MySQL Reference Manual: Doc-
umentation from Source. O'Reilly Community Press. O'Reilly and Associates, Inc.
3.4.1
YANG, L., LIN, B., KASHEFIPOUR, S.M. AND FALCONER, R.A. (2002). Integration
of a 1-D River Model with Object-Oriented Methodology. Environmental Modelling
and Software, 17, 693-701. 3.3.1
286
REFERENCES
YOON, T .H. AND KANG, S.K. (2004). Finite Volume Model for Two-Dimensional Shal-
low Water Flows on Unstructured Grids. Journal of Hydmulic Engineering, ASCE,
130, 678-688. 2.4.4, 4.6.3.2
YOST, S.A. AND RAO, P. (2000). A Moving Boundary Approach for One-Dimensional
Free Surface Flows. Advances in Water Resources, 23, 373-382. 2.4.2
ZANUTT!GH, B. AND LAMBERT!, A. (2002b). Roll Waves Simulation using Shallow Wa-
ter Equations and Weighted Average Flux Method. Journal of Hydmulic Research,
IAHR, 40, 610-622. 2.4.4
ZAWODNY, J.D. AND BALLING, D.J. (2004). High Performance MySQL Optimisation,
Backups, Replication and Load Balancing. O'Reilly and Associates. 3.4.1
ZHAO, D.H., SHEN, H.W., LAI, J.S. AND TABIOS In, G.Q. (1996). Approximate
Riemann Solvers in FVM for 2D Hydraulic Shock Wave Modelling. Journal of Hy-
dmulic Engineering, ASCE, 122, 692-702. 1.1
ZHOU, J.G. AND STANSBY, P.K. (1999). 2D Shallow Water Flow Model for the Hy-
draulic Jump. International Journal for Numerical Methods in Fluids, 29, 375-387.
2.4.4
287
REFERENCES
288