Pipe: Implement a improved handling of tapering.

 ref: #63

   This  extends the simple pipe momentum equation to include
      a local loss term F_taper = K_c · (ρ·v_o^2/2)·A_o for the diameter contraction.
      Adjust K_c based on your taper geometry: 0.05–0.15 for gentle cones,
      up to 0.6 for sharp contractions.

Author: dietmarw

OpenHPL/Waterway/Pipe.mo

@@ -15,6 +15,8 @@ model Pipe "Model of a pipe"
     Dialog(group = "Geometry"));
   parameter SI.Height p_eps = data.p_eps "Pipe roughness height" annotation (
     Dialog(group = "Geometry"));
+  parameter Real K_c = 0.1 "Loss coefficient for contraction"
+    annotation (Dialog(group = "Geometry"));
   // Steady state:
   parameter Boolean SteadyState=data.SteadyState "If true, starts in steady state" annotation (Dialog(group="Initialization"));
   parameter SI.VolumeFlowRate Vdot_0=data.Vdot_0 "Initial flow rate of the pipe" annotation (Dialog(group="Initialization"));
@@ -26,7 +28,9 @@ model Pipe "Model of a pipe"
   SI.Area A_ = D_ ^ 2 * C.pi / 4 "Average cross-sectional area";
   Real cos_theta = H / L "Slope ratio";
   SI.Velocity v "Water velocity";
+  SI.Velocity v_o "Outlet water velocity";
   SI.Force F_f "Friction force";
+  SI.Force F_taper "Tapering (local contraction) loss";
   SI.Momentum M "Water momentum";
   SI.Pressure p_i "Inlet pressure";
   SI.Pressure p_o "Outlet pressure";
@@ -41,18 +45,32 @@ model Pipe "Model of a pipe"
   SI.Temperature T( start = T_0);
   */
 
+protected
+  parameter SI.Diameter D_eff=
+    if abs(D_i - D_o)< C.eps then
+      D_i
+    else
+      (D_i - D_o) / log(D_i/D_o) "Effective diameter for a linear taper";
+
 initial equation
   if SteadyState then
     der(M) = 0;
   end if;
 equation
   Vdot = mdot / data.rho "Volumetric flow rate through the pipe";
-  v = Vdot / A_ "Water velocity";
+  v = Vdot / A_ "Average water velocity";
+  v_o = Vdot / A_o "Outlet water velocity";
   M = data.rho * L * Vdot "Momentum of water";
   m = data.rho * A_ * L "Mass of water";
-  F_f = Functions.DarcyFriction.Friction(v, D_, L, data.rho, data.mu, p_eps) "Friction force";
-  der(M) = data.rho * Vdot ^ 2 * (1 / A_i - 1 / A_o) + p_i * A_i - p_o * A_o - F_f + m * data.g * cos_theta
-   "momentum balance";
+  F_f = Functions.DarcyFriction.Friction(v, D_eff, L, data.rho, data.mu, p_eps)
+    "Friction force";
+  F_taper = K_c * 0.5 * data.rho * A_o * v_o * abs(v_o)
+    "Tapering (local contraction) loss";
+  der(M) = data.rho * Vdot^2 * (1/A_i - 1/A_o)
+           + p_i * A_i - p_o * A_o
+           - F_f - F_taper
+           + m * data.g * cos_theta
+    "Momentum balance including tapering loss";
   p_i = i.p "Inlet pressure";
   p_o = o.p "Outlet pressure";
 
@@ -81,6 +99,9 @@ equation
     (inlet and outlet pressure from connectors).</p>
     <p>It should be noted that this pipe model provides possibilities for modelling
     of pipes with both a positive and a negative slopes (positive or negative height difference).</p>
+    <p>If the pipe is slightly tapered then this can be taken into account by adjusting
+    <code>K_c</code> based on your taper geometry: 0.05–0.15 for gentle cones,
+      up to 0.6 for sharp contractions.</p>
     <p>More info about the pipe model can be found in
         <a href=\"modelica://OpenHPL.UsersGuide.References\">[Vytvytskyi2017]</a>
     and <a href=\"modelica://OpenHPL.UsersGuide.References\">[Splavska2017a]</a>.</p>