Modeling Waves in Rivers and Canals

../../_images/Lake-Hartwell-24012016.jpg

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:

../../_images/single_channel.png

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:

../../_images/single_channel_hec-ras.png

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:

../../_images/channel_impulse_xs.png

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)

../../_images/channel_pulse_results.svg

Both HEC-RAS and RTC-Tools were run with a spatial step size of 1000 m and a temporal step size of 15 min.