Skip to content

Commit 7f62d43

Browse files
committed
Extracts the Manning velocity calculation into a dedicated, reusable function
to improve modularity. Introduces the option to specify the Manning friction coefficient 'n' or the Strickler coefficient 'M' (M = 1/n), providing more flexibility for model configuration. Simplifies the channel bed slope definition by replacing separate inlet and outlet bed elevations with a single total elevation difference 'H'.
1 parent 5c50497 commit 7f62d43

5 files changed

Lines changed: 51 additions & 52 deletions

File tree

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,14 @@
1+
within OpenHPL.Functions;
2+
function manningVelocity "Compute velocity from Manning's equation"
3+
extends Modelica.Icons.Function;
4+
input SI.Height h "Water depth";
5+
input Real S "Slope (bed slope + water surface gradient)";
6+
input SI.Length w "Channel width";
7+
input Real n "Manning's roughness coefficient";
8+
output SI.Velocity v "Flow velocity";
9+
protected
10+
SI.Length R_h "Hydraulic radius";
11+
algorithm
12+
R_h := w * h / (w + 2 * h);
13+
v := sign(S) * R_h ^ (2.0 / 3) * abs(S) ^ 0.5 / n;
14+
end manningVelocity;

OpenHPL/Functions/package.order

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,4 @@
11
Fitting
22
DarcyFriction
3+
manningVelocity
34
KP07

OpenHPL/Waterway/OpenChannel.mo

Lines changed: 34 additions & 50 deletions
Original file line numberDiff line numberDiff line change
@@ -7,12 +7,15 @@ model OpenChannel "Open channel model with optional spatial discretization"
77
// Geometry
88
parameter SI.Length L = 5000 "Channel length" annotation (Dialog(group = "Geometry"));
99
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"));
10+
parameter SI.Length H = 10 "Bed elevation difference between inlet and outlet (positive = downhill inlet to outlet)" annotation (Dialog(group = "Geometry"));
1211

1312
// Friction
14-
parameter Real n_manning(unit = "s/m(1/3)") = 0.03 "Manning's roughness coefficient n"
15-
annotation (Dialog(group = "Friction"));
13+
parameter Real m_manning(unit="m(1/3)/s", min=0) = 33 "Manning M (Strickler) coefficient M=1/n (typically 60-110 for steel, 30-60 for rock tunnels)" annotation (
14+
Dialog(group = "Friction", enable = not use_n));
15+
parameter Boolean use_n = true "If true, use Mannings coefficient n (=1/M) instead of Manning's M (Strickler)" annotation (
16+
Dialog(group = "Friction"), choices(checkBox=true));
17+
parameter Real n_manning(unit="s/m(1/3)", min=0) = 0.03 "Manning's n coefficient (typically 0.009-0.017 for steel/concrete, 0.017-0.030 for rock tunnels)" annotation (
18+
Dialog(group = "Friction", enable = use_n));
1619

1720
// Discretization
1821
parameter Boolean useSections = false "If true, discretize the channel into N sections with varying water levels"
@@ -42,22 +45,9 @@ model OpenChannel "Open channel model with optional spatial discretization"
4245
SI.Velocity v_sec[if useSections then N else 0] "Velocity in each section";
4346

4447
protected
45-
parameter SI.Height dH = H_i - H_o "Total height difference (positive = downhill inlet to outlet)";
48+
parameter Real n_eff(unit="s/m(1/3)") = if use_n then n_manning else 1/m_manning "Effective Manning's n coefficient";
4649
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;
50+
parameter SI.PerUnit slope = H / L "Bed slope";
6151

6252
initial equation
6353
if SteadyState then
@@ -76,9 +66,8 @@ initial equation
7666
end if;
7767

7868
equation
79-
// ----- Connector pressures -----
80-
p_i = i.p;
81-
p_o = o.p;
69+
p_i = i.p "Inlet connector pressure";
70+
p_o = o.p "Outlet connector pressure";
8271

8372
// ----- Mass balance: incompressible, no storage in bulk -----
8473
i.mdot + o.mdot = 0;
@@ -88,52 +77,48 @@ equation
8877
if useSections then
8978
// ===== Sectional mode: N sections with individual water levels =====
9079

91-
// Average depth and velocity from sections
92-
h_avg = sum(h_sec) / N;
93-
v = Vdot / (W * h_avg);
80+
h_avg = sum(h_sec) / N "Average depth from sections";
81+
v = Vdot / (W * h_avg) "Average velocity from sections";
9482

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);
83+
F_f = data.rho * data.g * n_eff ^ 2 * v * abs(v) * L
84+
* (W + 2 * h_avg) ^ (4.0 / 3) / (W * h_avg) ^ (4.0 / 3)
85+
"Friction for overall momentum (Manning formula over full length)";
9886

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;
87+
L * der(mdot) = (p_i - p_o) * W * h_avg + data.rho * data.g * H * W * h_avg - F_f
88+
"Overall momentum balance determines flow rate";
10189

102-
// Section velocities
10390
for j in 1:N loop
104-
v_sec[j] = Vdot / (W * h_sec[j]);
91+
v_sec[j] = Vdot / (W * h_sec[j]) "Section velocities";
10592
end for;
10693

10794
// 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
95+
// W * dx * H/dt = Vdot_in - Vdot_out
96+
// where Vdot_in is driven by upstream depth and Vdot_out by downstream depth
11097
// using Manning equation locally for inter-section flow
11198
for j in 1:N loop
11299
W * dx * der(h_sec[j]) =
113100
(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))
101+
else W * h_sec[j - 1] * Functions.manningVelocity(h_sec[j - 1], slope
102+
+ (h_sec[j - 1] - h_sec[j]) / dx, W, n_eff))
116103
-
117104
(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));
105+
else W * h_sec[j] * Functions.manningVelocity(h_sec[j], slope
106+
+ (h_sec[j] - h_sec[j + 1]) / dx, W, n_eff));
120107
end for;
121108

122109
else
123110
// ===== Lumped mode: single control volume =====
124111

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);
112+
h_avg = max((p_i + p_o) / (2 * data.rho * data.g), 0.01)
113+
"Water depth from connector pressures (average of inlet and outlet)";
127114

128115
v = Vdot / (W * h_avg);
129116

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-
117+
F_f = data.rho * data.g * n_eff ^ 2 * v * abs(v) * L
118+
* (W + 2 * h_avg) ^ (4.0 / 3) / (W * h_avg) ^ (4.0 / 3)
119+
"Friction force using Manning equation for the full channel";
120+
L * der(mdot) = (p_i - p_o) * W * h_avg + data.rho * data.g * H * W * h_avg - F_f
121+
"Momentum balance";
137122
end if;
138123

139124
annotation (
@@ -146,8 +131,8 @@ The channel connects an upstream and downstream component via standard
146131
147132
<h5>Geometry</h5>
148133
<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>
134+
the difference in bed elevations between inlet and outlet <code>H</code>.
135+
The bed slope is computed as H/L.</p>
151136
152137
<h5>Governing Equations</h5>
153138
<p>The model is based on the momentum balance for incompressible flow:</p>
@@ -174,5 +159,4 @@ local water surface slope. This captures water level variations along the channe
174159
<p>Inlet and outlet connectors carry pressure and mass flow rate.
175160
Connect upstream to the inlet <code>i</code> and downstream to the outlet <code>o</code>.</p>
176161
</html>"));
177-
178162
end OpenChannel;

OpenHPL/Waterway/package.order

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,6 @@ Gate_HR
1212
RunOff_zones
1313
RunOff_zones_input
1414
VolumeFlowSource
15-
Penstock
1615
OpenChannel
16+
Penstock
1717
ReservoirChannel

OpenHPLTest/OpenChannel.mo

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -17,7 +17,7 @@ model OpenChannel "Model of a hydropower system with open channel model"
1717
origin={-90,90},
1818
extent={{-10,-10},{10,10}})));
1919
OpenHPL.Waterway.OpenChannel openChannel(
20-
H_i=2,
20+
H=2,
2121
useSections=true, N=10) annotation (Placement(transformation(extent={{-10,-10},{10,10}})));
2222
OpenHPL.Waterway.Pipe pipe(H=0, L=10) annotation (Placement(transformation(extent={{20,-10},{40,10}})));
2323
equation

0 commit comments

Comments
 (0)