feat: New version of OpenChannel

OpenHPL/Functions/manningVelocity.mo

@@ -0,0 +1,14 @@
+within OpenHPL.Functions;
+function manningVelocity "Compute velocity from Manning's equation"
+  extends Modelica.Icons.Function;
+  input SI.Height h "Water depth";
+  input Real S "Slope (bed slope + water surface gradient)";
+  input SI.Length w "Channel width";
+  input Real n "Manning's roughness coefficient";
+  output SI.Velocity v "Flow velocity";
+protected
+  SI.Length R_h "Hydraulic radius";
+algorithm
+  R_h := w * h / (w + 2 * h);
+  v := sign(S) * R_h ^ (2.0 / 3) * abs(S) ^ 0.5 / n;
+end manningVelocity;

OpenHPL/Functions/package.order

@@ -1,3 +1,4 @@
 Fitting
 DarcyFriction
+manningVelocity
 KP07

OpenHPL/Waterway/OpenChannel.mo

@@ -1,91 +1,162 @@
 within OpenHPL.Waterway;
-model OpenChannel "Open channel model (use KP scheme)"
-  extends Modelica.Icons.UnderConstruction;
+model OpenChannel "Open channel model with optional spatial discretization"
   outer Data data "Using standard data set";
   extends OpenHPL.Icons.OpenChannel;
-  // geometrical parameters of the open channel
-  parameter Integer N = 100 "Number of segments" annotation (Dialog(group = "Geometry"));
-  parameter SI.Length W=180 "Channel width" annotation (Dialog(group="Geometry"));
-  parameter SI.Length L = 5000 "Channel length" annotation (Dialog(group = "Geometry"));
-  parameter SI.Height H[2] = {17.5, 0} "Channel bed geometry, height from the left and right sides" annotation (Dialog(group = "Geometry"));
-  parameter Real f_n = 0.04 "Manning's roughness coefficient [s/m^1/3]" annotation (Dialog(group = "Geometry"));
-  parameter Boolean SteadyState=data.SteadyState "If true, starts in steady state" annotation (Dialog(group="Initialization"));
-  parameter SI.Height h_0[N]=ones(N)*5 "Initial water level" annotation (Dialog(group="Initialization"));
-  parameter SI.VolumeFlowRate Vdot_0=data.Vdot_0 "Initial flow rate" annotation (Dialog(group="Initialization"));
-  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"));
-  // variables
-  SI.VolumeFlowRate Vdot_o "Outlet flow";
-  SI.VolumeFlowRate Vdot_i "Inlet flow rate";
-  SI.Height h[N] "Water level in each unit of the channel";
-  // connector
   extends OpenHPL.Interfaces.TwoContacts;
-  // using open channel example from KP method class
-  Internal.KPOpenChannel openChannel(
-    N=N,
-    W=W,
-    L=L,
-    Vdot_0=Vdot_0,
-    f_n=f_n,
-    h_0=h_0,
-    boundaryValues=[h_0[1] + H[1],Vdot_i/W; h_0[N] + H[2],Vdot_o/W],
-    boundaryCondition=BoundaryCondition,
-    SteadyState=SteadyState) annotation (Placement(transformation(extent={{-10,-8},{10,12}})));
+
+  // Geometry
+  parameter SI.Length L = 5000 "Channel length" annotation (Dialog(group = "Geometry"));
+  parameter SI.Length W = 10 "Channel width" annotation (Dialog(group = "Geometry"));
+  parameter SI.Length H = 10 "Bed elevation difference between inlet and outlet (positive = downhill inlet to outlet)" annotation (Dialog(group = "Geometry"));
+
+  // Friction
+  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 (
+    Dialog(group = "Friction", enable = not use_n));
+  parameter Boolean use_n = true "If true, use Mannings coefficient n (=1/M) instead of Manning's M (Strickler)" annotation (
+    Dialog(group = "Friction"), choices(checkBox=true));
+  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 (
+    Dialog(group = "Friction", enable = use_n));
+
+  // Discretization
+  parameter Boolean useSections = false "If true, discretize the channel into N sections with varying water levels"
+    annotation (choices(checkBox = true), Dialog(group = "Discretization"));
+  parameter Integer N = 10 "Number of sections (only used when useSections = true)"
+    annotation (Dialog(group = "Discretization", enable = useSections));
+
+  // Initialization
+  parameter Boolean SteadyState = data.SteadyState "If true, starts in steady state"
+    annotation (Dialog(group = "Initialization"));
+  parameter SI.VolumeFlowRate Vdot_0 = data.Vdot_0 "Initial volume flow rate"
+    annotation (Dialog(group = "Initialization"));
+  parameter SI.Height h_0 = 2 "Initial water depth"
+    annotation (Dialog(group = "Initialization"));
+
+  // Variables — lumped (always computed)
+  SI.MassFlowRate mdot "Mass flow rate through the channel";
+  SI.VolumeFlowRate Vdot "Volume flow rate";
+  SI.Velocity v "Average water velocity";
+  SI.Height h_avg "Average water depth in the channel";
+  SI.Pressure p_i "Inlet pressure";
+  SI.Pressure p_o "Outlet pressure";
+  SI.Force F_f "Friction force";
+
+  // Variables — sectional (only meaningful when useSections = true)
+  SI.Height h_sec[if useSections then N else 0] "Water depth in each section";
+  SI.Velocity v_sec[if useSections then N else 0] "Velocity in each section";
+
+protected
+  parameter Real n_eff(unit="s/m(1/3)") = if use_n then n_manning else 1/m_manning "Effective Manning's n coefficient";
+  parameter SI.Length dx = L / max(N, 1) "Section length";
+  parameter SI.PerUnit slope = H / L "Bed slope";
+
+initial equation
+  if SteadyState then
+    der(mdot) = 0;
+    if useSections then
+      for j in 1:N loop
+        der(h_sec[j]) = 0;
+      end for;
+    end if;
+  else
+    if useSections then
+      for j in 1:N loop
+        h_sec[j] = h_0;
+      end for;
+    end if;
+  end if;
+
 equation
-// define a vector of the water depth in the channel
-  h = openChannel.h;
-// flow rate boundaries
-  i.mdot =Vdot_i*data.rho;
-  o.mdot =-Vdot_o*data.rho;
-// presurre boundaries
-  i.p = h[1] * data.g * data.rho + data.p_a;
-  o.p = h[N] * data.g * data.rho + data.p_a;
-  annotation (preferredView="info",
-    Documentation(info="<html>
-<p style=\"color: #ff0000;\"><em>Note: Currently under investigation for plausibility.</em></p>
+  p_i = i.p "Inlet connector pressure";
+  p_o = o.p "Outlet connector pressure";
 
-<h4>Open Channel Model</h4>
-<p>Model for open channels (rivers) that can be used for modeling run-of-river hydropower plants.
-The channel inlet and outlet are assumed to be at the bottom of the left and right sides, respectively.</p>
+  // ----- Mass balance: incompressible, no storage in bulk -----
+  i.mdot + o.mdot = 0;
+  mdot = i.mdot;
+  Vdot = mdot / data.rho;
 
-<h5>Governing Equations</h5>
-<p>The open channel model is based on the following partial differential equation:</p>
-<p>$$ \\frac{\\partial U}{\\partial t}+\\frac{\\partial F}{\\partial x} = S $$</p>
-<p>where:</p>
-<ul>
-<li>\\(U=\\left[\\begin{matrix}q & z\\end{matrix}\\right]^T\\)</li>
-<li>\\(F=\\left[\\begin{matrix}q & \\frac{q^2}{z-B}+\\frac{g}{2}\\left(z-B\\right)^2\\end{matrix}\\right]^T\\)</li>
-<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>
-</ul>
-<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,
-q is the discharge per unit width w of the open channel. f<sub>n</sub> is the Manning's roughness coefficient.</p>
+  if useSections then
+    // ===== Sectional mode: N sections with individual water levels =====
+
+    h_avg = sum(h_sec) / N "Average depth from sections";
+    v = Vdot / (W * h_avg) "Average velocity from sections";
+
+    F_f = data.rho * data.g * n_eff ^ 2 * v * abs(v) * L
+          * (W + 2 * h_avg) ^ (4.0 / 3) / (W * h_avg) ^ (4.0 / 3)
+          "Friction for overall momentum (Manning formula over full length)";
+
+    L * der(mdot) = (p_i - p_o) * W * h_avg + data.rho * data.g * H * W * h_avg - F_f
+                    "Overall momentum balance determines flow rate";
+
+    for j in 1:N loop
+      v_sec[j] = Vdot / (W * h_sec[j]) "Section velocities";
+    end for;
+
+    // Water level dynamics per section: continuity for free surface
+    //   W * dx * H/dt = Vdot_in - Vdot_out
+    // where Vdot_in is driven by upstream depth and Vdot_out by downstream depth
+    // using Manning equation locally for inter-section flow
+    for j in 1:N loop
+      W * dx * der(h_sec[j]) =
+        (if j == 1 then Vdot
+         else W * h_sec[j - 1] * Functions.manningVelocity(h_sec[j - 1], slope
+              + (h_sec[j - 1] - h_sec[j]) / dx, W, n_eff))
+        -
+        (if j == N then Vdot
+         else W * h_sec[j] * Functions.manningVelocity(h_sec[j], slope
+              + (h_sec[j] - h_sec[j + 1]) / dx, W, n_eff));
+    end for;
+
+  else
+    // ===== Lumped mode: single control volume =====
 
-<h5>Eigenvalues</h5>
-<p>The eigenvalues for this model are defined as:</p>
-<p>$$ \\lambda_{1,2}=u\\pm\\sqrt{gh} $$</p>
-<p>where u is the cross-section average water velocity.</p>
+    h_avg = max((p_i + p_o) / (2 * data.rho * data.g), 0.01)
+            "Water depth from connector pressures (average of inlet and outlet)";
 
-<h5>Desingularization</h5>
-<p>In dry or nearly dry channel areas, velocity at cell centers is recomputed using the desingularization formula:</p>
-<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>
-<p>applied when \\(h_{i\\pm\\frac{1}{2}}^\\pm<\\epsilon\\) (typically ε = 1e⁻⁵).</p>
+    v = Vdot / (W * h_avg);
 
-<h5>Implementation</h5>
-<p>Similar to <a href=\"modelica://OpenHPL.Waterway.PenstockKP\">PenstockKP</a>, this model uses the KP method
-(<a href=\"modelica://OpenHPL.Functions.KP07.KPmethod\">KPmethod</a> function) to discretize the PDEs into ODEs.</p>
+    F_f = data.rho * data.g * n_eff ^ 2 * v * abs(v) * L
+          * (W + 2 * h_avg) ^ (4.0 / 3) / (W * h_avg) ^ (4.0 / 3)
+          "Friction force using Manning equation for the full channel";
+    L * der(mdot) = (p_i - p_o) * W * h_avg + data.rho * data.g * H * W * h_avg - F_f
+                    "Momentum balance";
+  end if;
+
+  annotation (
+    Documentation(info="<html>
+<h4>Open Channel Model</h4>
+
+<p>Model for open channels (rivers, canals) suitable for run-of-river hydropower plants.
+The channel connects an upstream and downstream component via standard
+<a href=\"modelica://OpenHPL.Interfaces.Contact\">Contact</a> connectors (pressure and mass flow rate).</p>
+
+<h5>Geometry</h5>
+<p>The channel is defined by its length <code>L</code>, width <code>W</code>, and
+the difference in bed elevations between inlet and outlet <code>H</code>.
+The bed slope is computed as H/L.</p>
+
+<h5>Governing Equations</h5>
+<p>The model is based on the momentum balance for incompressible flow:</p>
+<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>
+<p>where A = W &middot; h is the cross-sectional flow area and F<sub>f</sub> is the friction force.</p>
 
-<p>Boundary conditions specify inlet and outlet flows per unit width q₁ and q₂.
-Connectors should be connected to <a href=\"modelica://OpenHPL.Waterway.Pipe\">Pipe</a> elements from both sides.
-Connectors provide inlet/outlet flow rates and pressures (sum of atmospheric pressure and water depth-dependent pressure).</p>
+<h5>Friction</h5>
+<p>Friction is computed using Manning's equation adapted for a rectangular cross-section:</p>
+<p>$$ F_\\mathrm{f} = \\rho\\,g\\,n^2\\,v\\,|v|\\,L\\,\\frac{(W + 2h)^{4/3}}{(W\\,h)^{4/3}} $$</p>
+<p>where n is Manning's roughness coefficient.</p>
 
-<h5>Parameters</h5>
+<h5>Modes of Operation</h5>
 <ul>
-<li>Geometry: channel length L and width w, bed height vector H at left/right sides</li>
-<li>Manning's roughness coefficient f<sub>n</sub></li>
-<li>Number of discretization cells N</li>
-<li>Initialization: initial flow rate \\(\\dot{V}_0\\) and water depth h₀ for each cell</li>
+<li><strong>Lumped mode</strong> (default): The channel is treated as a single control volume.
+The flow rate responds to the pressure difference between inlet and outlet connectors,
+gravity, and friction.</li>
+<li><strong>Sectional mode</strong> (<code>useSections = true</code>): The channel is divided into
+<code>N</code> sections. Each section maintains its own water depth via a continuity equation
+for the free surface. Inter-section flows are computed using Manning's equation with the
+local water surface slope. This captures water level variations along the channel.</li>
 </ul>
 
-<p><em>Note: This model is still under discussion and has not been tested properly.</em></p>
-<p>More details in <a href=\"modelica://OpenHPL.UsersGuide.References\">[Vytvytskyi2015]</a>.</p>
+<h5>Connectors</h5>
+<p>Inlet and outlet connectors carry pressure and mass flow rate.
+Connect upstream to the inlet <code>i</code> and downstream to the outlet <code>o</code>.</p>
 </html>"));
 end OpenChannel;

OpenHPL/Waterway/package.order

@@ -12,6 +12,6 @@ Gate_HR
 RunOff_zones
 RunOff_zones_input
 VolumeFlowSource
-Penstock
 OpenChannel
+Penstock
 ReservoirChannel

OpenHPLTest/OpenChannel.mo

@@ -16,12 +16,14 @@ model OpenChannel "Model of a hydropower system with open channel model"
                           annotation (Placement(transformation(
         origin={-90,90},
         extent={{-10,-10},{10,10}})));
-  OpenHPL.Waterway.OpenChannel openChannel(N=10, H={0,0}) annotation (Placement(transformation(extent={{-10,-10},{10,10}})));
+  OpenHPL.Waterway.OpenChannel openChannel(
+    H=2,
+    useSections=true,                      N=10) annotation (Placement(transformation(extent={{-10,-10},{10,10}})));
   OpenHPL.Waterway.Pipe pipe(H=0, L=10) annotation (Placement(transformation(extent={{20,-10},{40,10}})));
 equation
   connect(discharge.o, openChannel.i) annotation (Line(points={{-20,0},{-10,0}}, color={0,128,255}));
   connect(openChannel.o, pipe.i) annotation (Line(points={{10,0},{20,0}}, color={0,128,255}));
   connect(reservoir.o, discharge.i) annotation (Line(points={{-60,0},{-40,0}}, color={0,128,255}));
   connect(pipe.o, tail.o) annotation (Line(points={{40,0},{60,0}}, color={0,128,255}));
-  annotation (experiment(StopTime=1000));
+  annotation (experiment(StopTime=10000));
 end OpenChannel;