Skip to content

Commit 269c442

Browse files
dietmarwdietmarw
committed
First draft of open channel model
1 parent 4799b1a commit 269c442

1 file changed

Lines changed: 161 additions & 75 deletions

File tree

‎OpenHPL/Waterway/OpenChannel.mo‎

Lines changed: 161 additions & 75 deletions
Original file line numberDiff line numberDiff line change
@@ -1,92 +1,178 @@
11
within OpenHPL.Waterway;
2-
model OpenChannel "Open channel model (use KP scheme)"
3-
extends Modelica.Icons.UnderConstruction;
2+
model OpenChannel "Open channel model with optional spatial discretization"
43
outer Data data "Using standard data set";
54
extends OpenHPL.Icons.OpenChannel;
6-
// geometrical parameters of the open channel
7-
parameter Integer N = 100 "Number of segments" annotation (Dialog(group = "Geometry"));
8-
parameter SI.Length W=180 "Channel width" annotation (Dialog(group="Geometry"));
9-
parameter SI.Length L = 5000 "Channel length" annotation (Dialog(group = "Geometry"));
10-
parameter SI.Height H[2] = {17.5, 0} "Channel bed geometry, height from the left and right sides" annotation (Dialog(group = "Geometry"));
11-
parameter Real f_n = 0.04 "Manning's roughness coefficient [s/m^1/3]" annotation (Dialog(group = "Geometry"));
12-
parameter Boolean SteadyState=data.SteadyState "If true, starts in steady state" annotation (Dialog(group="Initialization"));
13-
parameter SI.Height h_0[N]=ones(N)*5 "Initial water level" annotation (Dialog(group="Initialization"));
14-
parameter SI.VolumeFlowRate Vdot_0=data.Vdot_0 "Initial flow rate" annotation (Dialog(group="Initialization"));
15-
parameter Boolean BoundaryCondition[2,2] = [false, true; false, true] "Boundary conditions. Choose options for the boundaries in a matrix table, i.e., if the matrix element = true, this element is used as boundary. The element represent the following quantities: [inlet depth, inlet flow; outlet depth, outlet flow]" annotation (Dialog(group = "Boundary condition"));
16-
// variables
17-
SI.VolumeFlowRate Vdot_o "Outlet flow";
18-
SI.VolumeFlowRate Vdot_i "Inlet flow rate";
19-
SI.Height h[N] "Water level in each unit of the channel";
20-
// connector
215
extends OpenHPL.Interfaces.TwoContacts;
22-
// using open channel example from KP method class
23-
Internal.KPOpenChannel openChannel(
24-
N=N,
25-
W=W,
26-
L=L,
27-
Vdot_0=Vdot_0,
28-
f_n=f_n,
29-
h_0=h_0,
30-
boundaryValues=[h_0[1] + H[1],Vdot_i/W; h_0[N] + H[2],Vdot_o/W],
31-
boundaryCondition=BoundaryCondition,
32-
SteadyState=SteadyState) annotation (Placement(transformation(extent={{-10,-8},{10,12}})));
6+
7+
// Geometry
8+
parameter SI.Length L = 5000 "Channel length" annotation (Dialog(group = "Geometry"));
9+
parameter SI.Length W = 10 "Channel width" annotation (Dialog(group = "Geometry"));
10+
parameter SI.Height H_i = 10 "Bed elevation at inlet" annotation (Dialog(group = "Geometry"));
11+
parameter SI.Height H_o = 0 "Bed elevation at outlet" annotation (Dialog(group = "Geometry"));
12+
13+
// Friction
14+
parameter Real n_manning(unit = "s/m(1/3)") = 0.03 "Manning's roughness coefficient n"
15+
annotation (Dialog(group = "Friction"));
16+
17+
// Discretization
18+
parameter Boolean useSections = false "If true, discretize the channel into N sections with varying water levels"
19+
annotation (choices(checkBox = true), Dialog(group = "Discretization"));
20+
parameter Integer N = 10 "Number of sections (only used when useSections = true)"
21+
annotation (Dialog(group = "Discretization", enable = useSections));
22+
23+
// Initialization
24+
parameter Boolean SteadyState = data.SteadyState "If true, starts in steady state"
25+
annotation (Dialog(group = "Initialization"));
26+
parameter SI.VolumeFlowRate Vdot_0 = data.Vdot_0 "Initial volume flow rate"
27+
annotation (Dialog(group = "Initialization"));
28+
parameter SI.Height h_0 = 2 "Initial water depth"
29+
annotation (Dialog(group = "Initialization"));
30+
31+
// Variables — lumped (always computed)
32+
SI.MassFlowRate mdot "Mass flow rate through the channel";
33+
SI.VolumeFlowRate Vdot "Volume flow rate";
34+
SI.Velocity v "Average water velocity";
35+
SI.Height h_avg "Average water depth in the channel";
36+
SI.Pressure p_i "Inlet pressure";
37+
SI.Pressure p_o "Outlet pressure";
38+
SI.Force F_f "Friction force";
39+
40+
// Variables — sectional (only meaningful when useSections = true)
41+
SI.Height h_sec[if useSections then N else 0] "Water depth in each section";
42+
SI.Velocity v_sec[if useSections then N else 0] "Velocity in each section";
43+
44+
protected
45+
parameter SI.Height dH = H_i - H_o "Total height difference (positive = downhill inlet to outlet)";
46+
parameter SI.Length dx = L / max(N, 1) "Section length";
47+
parameter Real slope(unit = "1") = dH / L "Bed slope";
48+
49+
function manningVelocity "Compute velocity from Manning's equation"
50+
input SI.Height h "Water depth";
51+
input Real S "Slope (bed slope + water surface gradient)";
52+
input SI.Length w "Channel width";
53+
input Real n "Manning's roughness coefficient";
54+
output SI.Velocity v "Flow velocity";
55+
protected
56+
SI.Length R_h "Hydraulic radius";
57+
algorithm
58+
R_h := w * h / (w + 2 * h);
59+
v := sign(S) * R_h ^ (2.0 / 3) * abs(S) ^ 0.5 / n;
60+
end manningVelocity;
61+
62+
initial equation
63+
if SteadyState then
64+
der(mdot) = 0;
65+
if useSections then
66+
for j in 1:N loop
67+
der(h_sec[j]) = 0;
68+
end for;
69+
end if;
70+
else
71+
if useSections then
72+
for j in 1:N loop
73+
h_sec[j] = h_0;
74+
end for;
75+
end if;
76+
end if;
77+
3378
equation
34-
// define a vector of the water depth in the channel
35-
h = openChannel.h;
36-
// flow rate boundaries
37-
i.mdot =Vdot_i*data.rho;
38-
o.mdot =-Vdot_o*data.rho;
39-
// presurre boundaries
40-
i.p = h[1] * data.g * data.rho + data.p_a;
41-
o.p = h[N] * data.g * data.rho + data.p_a;
42-
o.elevation.z = i.elevation.z - (H[1] - H[2]) "Elevation propagation: channel bed drops from H[1] to H[2]";
43-
annotation (preferredView="info",
44-
Documentation(info="<html>
45-
<p style=\"color: #ff0000;\"><em>Note: Currently under investigation for plausibility.</em></p>
79+
// ----- Connector pressures -----
80+
p_i = i.p;
81+
p_o = o.p;
4682

47-
<h4>Open Channel Model</h4>
48-
<p>Model for open channels (rivers) that can be used for modeling run-of-river hydropower plants.
49-
The channel inlet and outlet are assumed to be at the bottom of the left and right sides, respectively.</p>
83+
// ----- Mass balance: incompressible, no storage in bulk -----
84+
i.mdot + o.mdot = 0;
85+
mdot = i.mdot;
86+
Vdot = mdot / data.rho;
5087

51-
<h5>Governing Equations</h5>
52-
<p>The open channel model is based on the following partial differential equation:</p>
53-
<p>$$ \\frac{\\partial U}{\\partial t}+\\frac{\\partial F}{\\partial x} = S $$</p>
54-
<p>where:</p>
55-
<ul>
56-
<li>\\(U=\\left[\\begin{matrix}q & z\\end{matrix}\\right]^T\\)</li>
57-
<li>\\(F=\\left[\\begin{matrix}q & \\frac{q^2}{z-B}+\\frac{g}{2}\\left(z-B\\right)^2\\end{matrix}\\right]^T\\)</li>
58-
<li>\\(S=\\left[\\begin{matrix}0 & -g\\left(z-B\\right)\\frac{\\partial B}{\\partial x}-\\frac{gf_n^2q|q|\\left(w+2\\left(z-B\\right)\\right)^\\frac{4}{3}}{w^\\frac{4}{3}}\\frac{1}{\\left(z-B\\right)^\\frac{7}{3}}\\end{matrix}\\right]^T\\)</li>
59-
</ul>
60-
<p>with: \\(z=h+B\\), and \\(q=\\frac{\\dot{V}}{w}\\). Here, h is water depth in the channel, B is the channel bed elevation,
61-
q is the discharge per unit width w of the open channel. f<sub>n</sub> is the Manning's roughness coefficient.</p>
88+
if useSections then
89+
// ===== Sectional mode: N sections with individual water levels =====
6290

63-
<h5>Eigenvalues</h5>
64-
<p>The eigenvalues for this model are defined as:</p>
65-
<p>$$ \\lambda_{1,2}=u\\pm\\sqrt{gh} $$</p>
66-
<p>where u is the cross-section average water velocity.</p>
91+
// Average depth and velocity from sections
92+
h_avg = sum(h_sec) / N;
93+
v = Vdot / (W * h_avg);
6794

68-
<h5>Desingularization</h5>
69-
<p>In dry or nearly dry channel areas, velocity at cell centers is recomputed using the desingularization formula:</p>
70-
<p>$$ \\bar{u}_j=\\frac{2\\bar{h}_j\\bar{q}_j}{\\bar{h}_j^2+\\max\\left(\\bar{h}_j^2,\\epsilon^2\\right)} $$</p>
71-
<p>applied when \\(h_{i\\pm\\frac{1}{2}}^\\pm<\\epsilon\\) (typically ε = 1e⁻⁵).</p>
95+
// Friction for overall momentum (Manning formula over full length)
96+
F_f = data.rho * data.g * n_manning ^ 2 * v * abs(v) * L
97+
* (W + 2 * h_avg) ^ (4.0 / 3) / (W * h_avg) ^ (4.0 / 3);
7298

73-
<h5>Implementation</h5>
74-
<p>Similar to <a href=\"modelica://OpenHPL.Waterway.PenstockKP\">PenstockKP</a>, this model uses the KP method
75-
(<a href=\"modelica://OpenHPL.Functions.KP07.KPmethod\">KPmethod</a> function) to discretize the PDEs into ODEs.</p>
99+
// Overall momentum balance determines flow rate
100+
L * der(mdot) = (p_i - p_o) * W * h_avg + data.rho * data.g * dH * W * h_avg - F_f;
76101

77-
<p>Boundary conditions specify inlet and outlet flows per unit width q₁ and q₂.
78-
Connectors should be connected to <a href=\"modelica://OpenHPL.Waterway.Pipe\">Pipe</a> elements from both sides.
79-
Connectors provide inlet/outlet flow rates and pressures (sum of atmospheric pressure and water depth-dependent pressure).</p>
102+
// Section velocities
103+
for j in 1:N loop
104+
v_sec[j] = Vdot / (W * h_sec[j]);
105+
end for;
80106

81-
<h5>Parameters</h5>
107+
// Water level dynamics per section: continuity for free surface
108+
// W * dx * dh/dt = Qdot_in - Qdot_out
109+
// where Qdot_in is driven by upstream depth and Qdot_out by downstream depth
110+
// using Manning equation locally for inter-section flow
111+
for j in 1:N loop
112+
W * dx * der(h_sec[j]) =
113+
(if j == 1 then Vdot
114+
else W * h_sec[j - 1] * manningVelocity(h_sec[j - 1], slope
115+
+ (h_sec[j - 1] - h_sec[j]) / dx, W, n_manning))
116+
-
117+
(if j == N then Vdot
118+
else W * h_sec[j] * manningVelocity(h_sec[j], slope
119+
+ (h_sec[j] - h_sec[j + 1]) / dx, W, n_manning));
120+
end for;
121+
122+
else
123+
// ===== Lumped mode: single control volume =====
124+
125+
// Water depth from connector pressures (average of inlet and outlet)
126+
h_avg = max((p_i + p_o) / (2 * data.rho * data.g), 0.01);
127+
128+
v = Vdot / (W * h_avg);
129+
130+
// Friction force using Manning equation for the full channel
131+
F_f = data.rho * data.g * n_manning ^ 2 * v * abs(v) * L
132+
* (W + 2 * h_avg) ^ (4.0 / 3) / (W * h_avg) ^ (4.0 / 3);
133+
134+
// Momentum balance
135+
L * der(mdot) = (p_i - p_o) * W * h_avg + data.rho * data.g * dH * W * h_avg - F_f;
136+
137+
end if;
138+
139+
annotation (
140+
Documentation(info="<html>
141+
<h4>Open Channel Model</h4>
142+
143+
<p>Model for open channels (rivers, canals) suitable for run-of-river hydropower plants.
144+
The channel connects an upstream and downstream component via standard
145+
<a href=\"modelica://OpenHPL.Interfaces.Contact\">Contact</a> connectors (pressure and mass flow rate).</p>
146+
147+
<h5>Geometry</h5>
148+
<p>The channel is defined by its length <code>L</code>, width <code>W</code>, and
149+
the bed elevations at inlet (<code>H_i</code>) and outlet (<code>H_o</code>).
150+
The bed slope is computed as (H_i &minus; H_o)/L.</p>
151+
152+
<h5>Governing Equations</h5>
153+
<p>The model is based on the momentum balance for incompressible flow:</p>
154+
<p>$$ L\\,\\frac{\\mathrm{d}\\dot{m}}{\\mathrm{d}t} = (p_\\mathrm{i} - p_\\mathrm{o})\\,A + \\rho\\,g\\,\\Delta H\\,A - F_\\mathrm{f} $$</p>
155+
<p>where A = W &middot; h is the cross-sectional flow area and F<sub>f</sub> is the friction force.</p>
156+
157+
<h5>Friction</h5>
158+
<p>Friction is computed using Manning's equation adapted for a rectangular cross-section:</p>
159+
<p>$$ F_\\mathrm{f} = \\rho\\,g\\,n^2\\,v\\,|v|\\,L\\,\\frac{(W + 2h)^{4/3}}{(W\\,h)^{4/3}} $$</p>
160+
<p>where n is Manning's roughness coefficient.</p>
161+
162+
<h5>Modes of Operation</h5>
82163
<ul>
83-
<li>Geometry: channel length L and width w, bed height vector H at left/right sides</li>
84-
<li>Manning's roughness coefficient f<sub>n</sub></li>
85-
<li>Number of discretization cells N</li>
86-
<li>Initialization: initial flow rate \\(\\dot{V}_0\\) and water depth h₀ for each cell</li>
164+
<li><strong>Lumped mode</strong> (default): The channel is treated as a single control volume.
165+
The flow rate responds to the pressure difference between inlet and outlet connectors,
166+
gravity, and friction.</li>
167+
<li><strong>Sectional mode</strong> (<code>useSections = true</code>): The channel is divided into
168+
<code>N</code> sections. Each section maintains its own water depth via a continuity equation
169+
for the free surface. Inter-section flows are computed using Manning's equation with the
170+
local water surface slope. This captures water level variations along the channel.</li>
87171
</ul>
88172
89-
<p><em>Note: This model is still under discussion and has not been tested properly.</em></p>
90-
<p>More details in <a href=\"modelica://OpenHPL.UsersGuide.References\">[Vytvytskyi2015]</a>.</p>
173+
<h5>Connectors</h5>
174+
<p>Inlet and outlet connectors carry pressure and mass flow rate.
175+
Connect upstream to the inlet <code>i</code> and downstream to the outlet <code>o</code>.</p>
91176
</html>"));
177+
92178
end OpenChannel;

0 commit comments

Comments
 (0)