feat: BaseValve and Turbine update to aviod p_a=0 issues and non-linear closing law.

.gitignore

@@ -48,3 +48,12 @@ OpenHPL/Resources/Documents/UsersGuide_src/auto/
 OpenHPL/Resources/Documents/Nomenclature.html
 
 *.html
+
+
+# LaTeX intermdeiate documnts
+*.aux
+*.bbl
+*.blg
+*.log
+*.synctex.gz
+*.toc

OpenHPL/ElectroMech/BaseClasses/BaseValve.mo

@@ -1,76 +1,63 @@
 within OpenHPL.ElectroMech.BaseClasses;
 partial model BaseValve "Simple hydraulic valve (base class)"
-  outer Data data "Using standard class with global parameters";
-  extends OpenHPL.Interfaces.TwoContacts;
+   extends OpenHPL.Interfaces.TwoContacts;
+  outer OpenHPL.Data data "Using standard class with global parameters";
 
-  parameter Boolean ValveCapacity = true "If checked the valve capacity C_v should be specified,
-    otherwise specify the nominal values (net head and flow rate)"
+  parameter Boolean ValveCapacity = true
+    "If checked, the valve capacity C_v should be specified,
+    otherwise specify the nominal values (net head and flow rate at nominal opening)"
     annotation (Dialog(group = "Nominal values"), choices(checkBox = true));
   parameter Real C_v = 1 "Valve capacity"
     annotation (Dialog(group = "Nominal values", enable = ValveCapacity));
   parameter SI.Height H_n = 100 "Nominal net head"
     annotation (Dialog(group = "Nominal values", enable = not ValveCapacity));
-  parameter SI.VolumeFlowRate Vdot_n = 3 "Nominal flow rate"
+  parameter SI.VolumeFlowRate Vdot_n = 3 "Nominal volume flow rate"
     annotation (Dialog(group = "Nominal values", enable = not ValveCapacity));
-  parameter SI.PerUnit u_n = 0.95 "Nominal opening"
+  parameter SI.PerUnit u_n = 1 "Nominal opening"
     annotation (Dialog(group = "Nominal values", enable = not ValveCapacity));
-  parameter Boolean ConstEfficiency = true "If checked the constant efficiency eta_h is used,
-    otherwise specify lookup table for efficiency"
-    annotation (Dialog(group = "Efficiency data"), choices(checkBox = true));
-  parameter SI.Efficiency eta_h = 0.9 "Hydraulic efficiency"
-    annotation (Dialog(group = "Efficiency data", enable = ConstEfficiency));
-  parameter Real lookup_table[:, :] = [0, 0.4; 0.2, 0.7; 0.5, 0.9; 0.95, 0.95; 1.0, 0.93]
-    "Look-up table for the turbine/valve efficiency, described by a table matrix,
-     where the first column is a pu value of the guide vane opening,
-     and the second column is a pu value of the turbine efficiency."
-    annotation (Dialog(group = "Efficiency data", enable = not ConstEfficiency));
-  parameter Boolean WaterCompress = false "If checked, the water is compressible"
-    annotation (Dialog(tab = "Advanced"), choices(checkBox = true));
+  parameter Real alpha=1 "Exponent of closing curve (if different from 1, the closing law will be non-linear)" annotation(Dialog(tab="Advanced", group = "Closing law"));
 
   SI.Pressure dp "Pressure drop";
-  SI.EnergyFlowRate Kdot_i_tr "Kinetic energy flow";
   SI.MassFlowRate mdot "Mass flow rate";
-  SI.VolumeFlowRate Vdot "Flow rate";
-  Real C_v_ "Valve capacity";
+  SI.VolumeFlowRate Vdot "Volume flow rate";
 
-  output SI.EnergyFlowRate Wdot_s "Valve power";
-  Modelica.Blocks.Tables.CombiTable1Dv look_up_table(table=lookup_table);
 protected
-  Modelica.Blocks.Interfaces.RealInput u(min=0, max=1)
-    "=1: completely open, =0: completely closed"
-  annotation (Placement(transformation(
-        origin={0,70},
-        extent={{-10,-10},{10,10}},
-        rotation=270), iconTransformation(
-        extent={{-20,-20},{20,20}},
-        rotation=270,
-        origin={0,80})));
+  parameter Real C_v_ = if ValveCapacity then C_v else Vdot_n/sqrt(H_n*data.g*data.rho)/u_n
+    "Define 'valve capacity' based on the nominal values";
+  Modelica.Blocks.Interfaces.RealInput u(min = 0, max = 1) "=1: completely open, =0: completely closed"
+    annotation (Placement(transformation(origin = {0, 70}, extent = {{-10, -10}, {10, 10}}, rotation = 270), iconTransformation(extent = {{-20, -20}, {20, 20}}, rotation = 270, origin = {0, 80})));
+  constant Real epsilon = 5.0e-5 "Constant to ensure robust expression for dp vs flow. Trial and error to find suitable value.";
+
 equation
-  i.mdot+o.mdot=0;
+  i.mdot + o.mdot = 0;
   mdot = i.mdot;
-  Vdot = if WaterCompress then mdot / (data.rho * (1 + data.beta * (i.p - data.p_a))) else mdot / data.rho
-    "Checking for water compressibility";
-  look_up_table.u[1] = u "Link the valve opening";
-  C_v_ = if ValveCapacity then C_v else Vdot_n/sqrt(H_n*data.g*data.rho/data.p_a)/u_n
-    "Define 'valve capacity' base on the nominal values";
-  dp = Vdot ^ 2 * data.p_a / (C_v_ * max(1e-6, u)) ^ 2 "Valve equation for pressure drop";
+  Vdot = mdot/data.rho;
+  dp*(C_v_*max(epsilon, u^alpha))^2 = Vdot*abs(Vdot) "Valve equation for pressure drop";
   dp = i.p - o.p "Link the pressure drop to the ports";
-  Kdot_i_tr = dp * Vdot "Energy balance";
-  if ConstEfficiency then
-    Wdot_s = eta_h * Kdot_i_tr;
-  else
-    Wdot_s = look_up_table.y[1] * Kdot_i_tr;
-  end if;
-
-annotation (
-    Documentation(info="<html><p>
-This is a simple model of hydraulic valve.
-The model can use a constant efficiency or varying
-efficiency from a lookup-table.
+  annotation ( Documentation(info="<html>
+<p>
+This is a partial, simple model of hydraulic valve. &nbsp;</p><p>This model is based on the energy balance of a valve.
+</p>
+<ul>
+<li>Mass flow is equal at innflow and outflow</li>
+<li>The head loss and pressure difference is proportional to square of velocity</li>
+</ul>
+<p>
+Specifically:
+</p>
+<pre>
+dp*f(opening)=v|v|
+</pre>
+<p>
+The function f(opening) is expressed as:
 </p>
+<pre>
+(C_v_*max(epsilon, u^alpha))^2
+</pre>
 <p>
-This model is based on the energy balance of a valve.
-The valve capacity can either be specified
+When alpha is 1, this implies a linear relation between closing and head loss.
+</p>
+<p>The valve capacity can either be specified
 directly by the user by specifying <code>C_v</code> or it will be calculated from
 the nominal head <code>H_n</code> and nominal flow rate <code>Vdot_n</code>.
 </p>

OpenHPL/ElectroMech/Turbines/Turbine.mo

@@ -5,19 +5,42 @@ model Turbine "Simple turbine model with mechanical connectors"
   extends Interfaces.TurbineContacts;
   extends Icons.Turbine;
 
-   Modelica.Blocks.Math.Feedback lossCorrection annotation (Placement(transformation(extent={{-50,70},{-30,90}})));
+  parameter Boolean ConstEfficiency = true
+    "If checked the constant efficiency eta_h is used,
+     otherwise specify lookup table for efficiency"
+    annotation (Dialog(group = "Efficiency data"), choices(checkBox = true));
+  parameter SI.Efficiency eta_h = 0.9 "Hydraulic efficiency"
+    annotation (Dialog(group = "Efficiency data", enable = ConstEfficiency));
+  parameter Real lookup_table[:, :] = [0, 0.4; 0.2, 0.7; 0.5, 0.9; 0.95, 0.95; 1.0, 0.93]
+    "Look-up table for the turbine/valve efficiency, described by a table matrix, where the first column is a pu value of the guide vane opening, and the second column is a pu value of the turbine efficiency."
+    annotation (Dialog(group = "Efficiency data", enable = not ConstEfficiency));
+  Modelica.Blocks.Math.Feedback lossCorrection
+  annotation (Placement(transformation(extent={{-50,70},{-30,90}})));
+  Modelica.Blocks.Tables.CombiTable1Dv look_up_table(table = lookup_table, smoothness = Modelica.Blocks.Types.Smoothness.ContinuousDerivative, extrapolation = Modelica.Blocks.Types.Extrapolation.LastTwoPoints) annotation(Placement(transformation(origin = {-76, -74}, extent={{-10,-10},{10,10}})));
+  output Modelica.Units.SI.EnergyFlowRate Wdot_s "Turbine power";
+
+protected
+  SI.EnergyFlowRate Kdot_i_tr "Gross hydraulic power";
+
 equation
+  look_up_table.u[1] = u "Link the valve opening";
+  if ConstEfficiency then
+    Wdot_s = eta_h * Kdot_i_tr;
+  else
+    Wdot_s = look_up_table.y[1] * Kdot_i_tr;
+  end if;
+  Kdot_i_tr = dp * Vdot "Energy balance";
 
   connect(P_out, lossCorrection.y) annotation (Line(
       points={{40,110},{40,80},{-31,80}},
       color={0,0,127},
       pattern=LinePattern.Dash));
   connect(lossCorrection.u1, power.y) annotation (Line(points={{-48,80},{-88,80},{-88,30},{-81,30}},color={0,0,127}));
-  connect(frictionLoss.power, lossCorrection.u2) annotation (Line(points={{-1,12},{-40,12},{-40,72}},
-                                                                                                  color={0,0,127}));
-  connect(u_t, u) annotation (Line(points={{-80,120},{-80,90},{0,90},{0,70}},    color={0,0,127}));
+  connect(frictionLoss.power, lossCorrection.u2) annotation (Line(points={{-1,12},{-40,12},{-40,72}},color={0,0,127}));
+  connect(u_t, u) annotation (Line(points={{-80,120},{-80,90},{0,90},{0,70}},color={0,0,127}));
   annotation (
-    Documentation(info="<html><p>
+    Documentation(info="<html>
+<p>
 This is a simple model of the turbine that give possibilities for simplified
 modelling of the turbine unit. The model can use a constant efficiency or varying
 efficiency from a lookup-table.
@@ -51,5 +74,5 @@ there are inputs as the control signal for the valve opening and also output as
           textString="P"), Text(
           extent={{-96,100},{-60,80}},
           textColor={0,0,0},
-          textString="GVO")}));
+          textString="Opening")}));
 end Turbine;

OpenHPL/Examples/SimpleValve.mo

@@ -34,9 +34,7 @@ model SimpleValve "Model of a hydropower system with a simple turbine turbine"
         extent={{-10,-10},{10,10}})));
   Waterway.Valve valve(
     ValveCapacity=false,
-    C_v=1,
-    ConstEfficiency=false,
-    WaterCompress=false) annotation (Placement(transformation(extent={{20,20},{40,40}})));
+    C_v=1)               annotation (Placement(transformation(extent={{20,20},{40,40}})));
 equation
   connect(reservoir.o, intake.i) annotation (
     Line(points={{-80,30},{-70,30}}, color = {28, 108, 200}));