Modeling Waves in Rivers and Canals

Note
This is a more advanced example that implements advanced channel hydraulics in RTC-Tools. It also capitalizes on the homotopy techniques available in RTC-Tools. If you are a first-time user of RTC-Tools, see Filling a Reservoir.
The RTC-Tools is capable of handling non-linear hydraulics. In this example, we model a river channel that receives a sudden pulse of higher-than-usual water volumes. We compare the results to those of an identical model built in HEC-RAS.
The Model
In this example, water is flowing through a single channel. There is an inflow at the upstream end and a water level bound at the downstream end.
In OpenModelica Connection Editor, the model looks like this in plan view:

In text mode, the Modelica model looks as follows (with annotation statements removed):
1model Example
2 // Elements
3 Deltares.ChannelFlow.Hydraulic.Branches.HomotopicTrapezoidal Channel(
4 Q_nominal = 100.0,
5 H_b_down = -5.0,
6 H_b_up = -5.0,
7 friction_coefficient = 0.045,
8 use_manning = true,
9 length = 10000,
10 theta = theta,
11 use_inertia = true,
12 use_convective_acceleration = false,
13 use_upwind = false,
14 n_level_nodes = 11,
15 uniform_nominal_depth = 5.0,
16 bottom_width_down = 30,
17 bottom_width_up = 30,
18 left_slope_angle_up = 45,
19 left_slope_angle_down = 45,
20 right_slope_angle_up = 45,
21 right_slope_angle_down = 45,
22 semi_implicit_step_size = step_size
23 ) ;
24 Deltares.ChannelFlow.Hydraulic.BoundaryConditions.Level Level;
25 Deltares.ChannelFlow.Hydraulic.BoundaryConditions.Discharge Discharge;
26 // Inputs
27 input Real Inflow_Q(fixed=true) = Discharge.Q;
28 input Real Level_H(fixed=true) = Level.H;
29 parameter Real theta;
30 parameter Real step_size;
31 // Output Channel states
32 output Real Channel_Q_up = Discharge.Q;
33 output Real Channel_Q_dn = Level.HQ.Q;
34 output Real Channel_H_up = Discharge.HQ.H;
35 output Real Channel_H_dn = Level.H;
36equation
37 connect(Channel.HQDown, Level.HQ);
38 connect(Discharge.HQ, Channel.HQUp);
39initial equation
40 Channel.Q = fill(Inflow_Q, Channel.n_level_nodes + 1);
41end Example;
The plan view of the model looks like this in HEC-RAS:

The channel cross-section is a simple trapezoidal shape. As rendered by HEC-RAS, here is a cross-section view of the channel being modeled:

The model was built with HEC-RAS version 5.0.6. In case you wish to verify the
HEC-RAS model yourself, a zip of the HEC-RAS model used in this comparison is
available: HEC-RAS.zip
The Python File
To keep this example simple and to allow for a 1:1 comparison with HEC-RAS, we will not have any decision variables in this model.
1from rtctools.optimization.collocated_integrated_optimization_problem import (
2 CollocatedIntegratedOptimizationProblem,
3)
4from rtctools.optimization.csv_mixin import CSVMixin
5from rtctools.optimization.homotopy_mixin import HomotopyMixin
6from rtctools.optimization.modelica_mixin import ModelicaMixin
7from rtctools.util import run_optimization_problem
8
9
10class Example(HomotopyMixin, CSVMixin, ModelicaMixin, CollocatedIntegratedOptimizationProblem):
11 def parameters(self, ensemble_member):
12 p = super().parameters(ensemble_member)
13 times = self.times()
14 if self.use_semi_implicit:
15 p["step_size"] = times[1] - times[0]
16 else:
17 p["step_size"] = 0.0
18 p["Channel.use_convective_acceleration"] = self.use_convective_acceleration
19 p["Channel.use_upwind"] = self.use_upwind
20 return p
21
22 def constraints(self, ensemble_member):
23 constraints = super().constraints(ensemble_member)
24 times = self.times()
25
26 # Extract the number of nodes in the channel
27 parameters = self.parameters(ensemble_member)
28 n_level_nodes = int(parameters["Channel.n_level_nodes"])
29
30 # To Mimic HEC-RAS behaviour, enforce steady state both at t0 and at t1.
31 for i in range(n_level_nodes):
32 state = "Channel.H[{}]".format(i + 1)
33 constraints.append(
34 (self.state_at(state, times[0]) - self.state_at(state, times[1]), 0, 0)
35 )
36 return constraints
37
38
39class ExampleInertialWave(Example):
40 """Inertial wave equation (no convective acceleration)"""
41
42 model_name = "Example"
43
44 use_semi_implicit = False
45 use_convective_acceleration = False
46 use_upwind = False
47
48 timeseries_export_basename = "timeseries_export_inertial_wave"
49
50
51class ExampleInertialWaveSemiImplicit(Example):
52 """Inertial wave equation (no convective acceleration)"""
53
54 model_name = "Example"
55
56 use_semi_implicit = True
57 use_convective_acceleration = False
58 use_upwind = False
59
60 timeseries_export_basename = "timeseries_export_inertial_wave_semi_implicit"
61
62
63class ExampleSaintVenant(Example):
64 """Saint Venant equation. Convective acceleration discretized with central differences"""
65
66 model_name = "Example"
67
68 use_semi_implicit = False
69 use_convective_acceleration = True
70 use_upwind = False
71
72 timeseries_export_basename = "timeseries_export_saint_venant"
73
74
75class ExampleSaintVenantUpwind(Example):
76 """Saint Venant equation. Convective acceleration discretized with upwind scheme"""
77
78 model_name = "Example"
79
80 use_semi_implicit = False
81 use_convective_acceleration = True
82 use_upwind = True
83
84 timeseries_export_basename = "timeseries_export_saint_venant_upwind"
85
86
87run_optimization_problem(ExampleInertialWave)
88run_optimization_problem(ExampleInertialWaveSemiImplicit)
89run_optimization_problem(ExampleSaintVenant)
90run_optimization_problem(ExampleSaintVenantUpwind)
As you can see, this model is as simple as it gets. We only add a constraint to keep the initialization states consistent with the HEC-RAS initialization.
Comparison of Discretizations and Numerical Schemes
HEC-RAS and RTC-Tools use different discretizations and numerical schemes, but also solve different equations. RTC-Tools solves the original nonlinear equations, whereas HEC-RAS solves a linearized momentum equation.
RTC-Tools 2 |
HEC-RAS |
|
---|---|---|
Momentum equation |
Saint-Venant / inertial wave (default) |
Linearized Saint-Venant |
Spatial discretization |
Staggered |
Collocated |
Numerical scheme (temporal) |
Semi-implicit / implicit (default) |
Centered Preissmann box scheme |
Numerical scheme (spatial) |
Central differences, upwind convective acceleration (optional) |
Centered Preissmann box scheme |
Note
For optimization, the recommended momentum equation and temporal scheme for RTC-Tools is semi-implicit inertial wave. Consult Baayen and Piovesan, A continuation approach to nonlinear model predictive control of open channel systems, 2018, for details. A preprint is available online as arXiv:1801.06507.
Comparison of Results
The results from the RTC-Tools run are found in the output directory with the
name timeseries_export.csv
, and the results generated by HEC-RAS have been
exported into the same directory under the name HEC-RAS_results.csv
. We can
compare the results using the Python library matplotlib
:
(Source code
, svg
, png
)
Both HEC-RAS and RTC-Tools were run with a spatial step size of 1000 m and a temporal step size of 15 min.