Computation of friction coefficient by replaceable model

.gitignore

@@ -0,0 +1 @@
+**/dsin.txt

ThermoPower/Choices.mo

@@ -25,6 +25,7 @@ package Choices "Choice enumerations for ThermoPower models"
         OpPoint "Friction factor defined by operating point",
         Cfnom "Cfnom friction factor",
         Colebrook "Colebrook's equation",
+        External "External equation using replaceable function",
         NoFriction "No friction")
       "Type, constants and menu choices to select the friction factor";
 

ThermoPower/Friction.mo

@@ -0,0 +1,86 @@
+within ThermoPower;
+package Friction "Friction models for 1D gas and water components"
+  extends Modelica.Icons.Package;
+
+  partial model FrictionModelBase "Base model for friction"
+
+    input Medium.ThermodynamicState[N] fluidState;
+    input Medium.MassFlowRate[N] w;
+    replaceable package Medium = Modelica.Media.Interfaces.PartialMedium
+      "Medium model"
+      annotation(Dialog(enable=false, tab = "Set by Flow1D model"));
+    parameter Integer N "Number of elements (volumes if FV, nodes if FEM)"
+      annotation(Dialog(enable=false, tab = "Set by Flow1D model"));
+    parameter SI.Distance L "Tube length"
+      annotation(Dialog(enable=false, tab = "Set by Flow1D model"));
+    parameter SI.Area A "Cross-sectional area (single tube)"
+      annotation(Dialog(enable=false, tab = "Set by Flow1D model"));
+    final parameter SI.Length omega_hyd = 4*A/Dhyd
+      "Wet perimeter (single tube)"
+      annotation(Dialog(enable=false, tab = "Set by Flow1D model"));
+    parameter SI.Length Dhyd "Hydraulic Diameter (single tube)"
+      annotation(Dialog(enable=false, tab = "Set by Flow1D model"));
+    parameter SI.MassFlowRate wnom "Nominal mass flow rate (single tube)"
+      annotation(Dialog(enable=false, tab = "Set by Flow1D model"));
+
+    output SI.CoefficientOfFriction[N] Cf "Fanning friction factor";
+
+    annotation (Icon(coordinateSystem(preserveAspectRatio=false)), Diagram(
+          coordinateSystem(preserveAspectRatio=false)));
+  end FrictionModelBase;
+
+  model NoFriction "Ideal flow model"
+    extends FrictionModelBase;
+
+  equation
+
+    Cf = zeros(N);
+
+  end NoFriction;
+
+  model NominalFriction "Constant friction factor flow model"
+    extends FrictionModelBase;
+
+    parameter SI.CoefficientOfFriction Cfnom "Nominal friction factor value";
+
+  equation
+
+    for i in 1:N loop
+      Cf[i] = Cfnom;
+    end for;
+
+  end NominalFriction;
+
+  model NominalOperatingPoint
+    "Flow model with friction factor defined by operating point"
+    extends FrictionModelBase;
+
+    parameter SI.PressureDifference dpnom "Nominal pressure drop";
+    parameter Medium.Density rhonom "Nominal density";
+
+  equation
+
+    for i in 1:N loop
+      Cf[i] = 2*dpnom*A^3/(rhonom*wnom^2*omega_hyd*L);
+    end for;
+
+  end NominalOperatingPoint;
+
+  model Colebrook "Colebrook flow model"
+    extends FrictionModelBase;
+
+    parameter SI.PerUnit e "Relative roughness (ratio roughness/diameter)";
+
+  equation
+
+    // use f_colebrook from the Water package to keep it backwards-compatible
+    // (the function in the Gas package is identical)
+    for i in 1:N loop
+      Cf[i] = Water.f_colebrook(w,
+        Dhyd/A,
+        e,
+        Medium.dynamicViscosity(fluidState[i]));
+    end for;
+
+  end Colebrook;
+end Friction;

ThermoPower/Gas.mo

@@ -964,13 +964,27 @@ package Gas "Models of components with ideal gases as working fluid"
       final w=w*ones(N),
       final fluidState=gas.state) "Instantiated heat transfer model";
 
+    replaceable model FrictionModel = Friction.NoFriction
+      constrainedby Friction.FrictionModelBase
+      annotation (choicesAllMatching=true, Dialog(enable = useFrictionModel));
+
+    FrictionModel frictionModel(
+      redeclare package Medium = Medium,
+      final N=N-1,
+      final L=L,
+      final A=A,
+      final Dhyd=Dhyd,
+      final wnom=wnom/Nt,
+      final w=fill(w,N-1),
+      final fluidState=gas[1:N-1].state)
+      "Instantiated flow model";
+
     parameter SI.PerUnit wnm = 1e-2 "Maximum fraction of the nominal flow rate allowed as reverse flow";
     parameter Boolean fixedMassFlowSimplified = false "Fix flow rate = wnom for simplified homotopy model"
         annotation (Dialog(tab="Initialisation"));
 
     Medium.BaseProperties gas[N] "Gas nodal properties";
     SI.Pressure Dpfric "Pressure drop due to friction";
-    SI.Length omega_hyd "Wet perimeter (single tube)";
     Real Kf "Friction factor";
     Real Kfl "Linear friction factor";
     Real dwdt "Time derivative of mass flow rate";
@@ -988,7 +1002,7 @@ package Gas "Models of components with ideal gases as working fluid"
     SI.Power Q_single[N-1] = heatTransfer.Qvol/Nt
       "Heat flows entering the volumes from the lateral boundary (single tube)";
     SI.Velocity u[N] "Fluid velocity";
-    Medium.AbsolutePressure p(start=pstart, stateSelect=StateSelect.prefer);
+    Medium.AbsolutePressure p(start=pstart, fixed = false, stateSelect=StateSelect.prefer);
     SI.Time Tr "Residence time";
     SI.Mass M "Gas Mass (single tube)";
     SI.Mass Mtot "Gas Mass (total)";
@@ -1019,40 +1033,39 @@ package Gas "Models of components with ideal gases as working fluid"
       "dpnom=0 not supported, it is also used in the homotopy trasformation during the inizialization");
     //All equations are referred to a single tube
     // Friction factor selection
-    omega_hyd = 4*A/Dhyd;
-    if FFtype == ThermoPower.Choices.Flow1D.FFtypes.Kfnom then
+    if useFrictionModel then
+      Cf = sum(frictionModel.Cf)/(N-1)*Kfc;
+    elseif FFtype == ThermoPower.Choices.Flow1D.FFtypes.Kfnom then
       Kf = Kfnom*Kfc;
-      Cf = 2*Kf*A^3/(omega_hyd*L);
     elseif FFtype == ThermoPower.Choices.Flow1D.FFtypes.OpPoint then
       Kf = dpnom*rhonom/(wnom/Nt)^2*Kfc;
-      Cf = 2*Kf*A^3/(omega_hyd*L);
     elseif FFtype == ThermoPower.Choices.Flow1D.FFtypes.Cfnom then
-      Kf = Cfnom*omega_hyd*L/(2*A^3)*Kfc;
       Cf = Cfnom*Kfc;
     elseif FFtype == ThermoPower.Choices.Flow1D.FFtypes.Colebrook then
       Cf = f_colebrook(
           w,
           Dhyd/A,
           e,
           Medium.dynamicViscosity(gas[integer(N/2)].state))*Kfc;
-      Kf = Cf*omega_hyd*L/(2*A^3);
     elseif FFtype == ThermoPower.Choices.Flow1D.FFtypes.NoFriction then
       Cf = 0;
-      Kf = 0;
+    elseif FFtype == ThermoPower.Choices.Flow1D.FFtypes.External then
+      Cf = Kfc*Cfcorr(w,Dhyd,A,sum(Medium.dynamicViscosity(gas[:].state))/N);
     else
       assert(false, "Unsupported FFtype");
       Cf = 0;
-      Kf = 0;
     end if;
     assert(Kf >= 0, "Negative friction coefficient");
+    Kf = Cf*omega_hyd*L/(2*A^3)
+      "Relationship between friction coefficient and Fanning friction factor";
     Kfl = wnom/Nt*wnf*Kf "Linear friction factor";
 
     // Dynamic momentum term
     dwdt = if DynamicMomentum and not QuasiStatic then der(w) else 0;
 
     sum(dMdt) = (infl.m_flow + outfl.m_flow)/Nt "Mass balance";
     L/A*dwdt + (outfl.p - infl.p) + Dpfric = 0 "Momentum balance";
-    Dpfric = (if FFtype == ThermoPower.Choices.Flow1D.FFtypes.NoFriction then 0
+    Dpfric = (if not useFrictionModel and (FFtype == ThermoPower.Choices.Flow1D.FFtypes.NoFriction) then 0
               else homotopy((smooth(1, Kf*squareReg(w, wnom/Nt*wnf))*sum(vbar)/(N - 1)),
                              dpnom/(wnom/Nt)*w))
       "Pressure drop due to friction";
@@ -2064,9 +2077,9 @@ package Gas "Models of components with ideal gases as working fluid"
     parameter Medium.MassFlowRate wnom "Nominal mass flowrate"
       annotation (Dialog(group="Nominal operating point"));
     parameter Medium.Density rhonom=1000 "Nominal density" annotation (Dialog(group=
-            "Nominal operating point", enable=(CvData == CvTypes.OpPoint)));
+            "Nominal operating point", enable=(CvData == ThermoPower.Choices.Valve.CvTypes.OpPoint)));
     parameter SI.PerUnit thetanom=1 "Nominal valve opening" annotation (Dialog(group=
-            "Nominal operating point", enable=(CvData == CvTypes.OpPoint)));
+            "Nominal operating point", enable=(CvData == ThermoPower.Choices.Valve.CvTypes.OpPoint)));
     parameter Boolean CheckValve=false "Reverse flow stopped";
     parameter SI.PerUnit b=0.01 "Regularisation factor";
 
@@ -2850,18 +2863,20 @@ The packages Medium are redeclared and a mass balance determines the composition
       parameter SI.Length omega
         "Perimeter of heat transfer surface (single tube)";
       parameter SI.Length Dhyd "Hydraulic Diameter (single tube)";
+      final parameter SI.Length omega_hyd = 4*A/Dhyd "Wet perimeter (single tube)";
       parameter Medium.MassFlowRate wnom "Nominal mass flowrate (total)";
       parameter ThermoPower.Choices.Flow1D.FFtypes FFtype=ThermoPower.Choices.Flow1D.FFtypes.NoFriction
         "Friction Factor Type"
-        annotation(Evaluate=true);
+        annotation(Evaluate=true, Dialog(enable = not useFrictionModel));
       parameter SI.PressureDifference dpnom = 0 "Nominal pressure drop";
       parameter Real Kfnom=0 "Nominal hydraulic resistance coefficient"
-        annotation(Dialog(enable = (FFtype == ThermoPower.Choices.Flow1D.FFtypes.Kfnom)));
+        annotation(Dialog(enable = (FFtype == ThermoPower.Choices.Flow1D.FFtypes.Kfnom) and not useFrictionModel));
       parameter Medium.Density rhonom=0 "Nominal inlet density"
-        annotation(Dialog(enable = (FFtype == ThermoPower.Choices.Flow1D.FFtypes.OpPoint)));
+        annotation(Dialog(enable = (FFtype == ThermoPower.Choices.Flow1D.FFtypes.OpPoint) and not useFrictionModel));
       parameter SI.PerUnit Cfnom=0 "Nominal Fanning friction factor"
-        annotation(Dialog(enable = (FFtype == ThermoPower.Choices.Flow1D.FFtypes.Cfnom)));
-      parameter SI.PerUnit e=0 "Relative roughness (ratio roughness/diameter)";
+        annotation(Dialog(enable = (FFtype == ThermoPower.Choices.Flow1D.FFtypes.Cfnom) and not useFrictionModel));
+      parameter SI.PerUnit e=0 "Relative roughness (ratio roughness/diameter)"
+        annotation(Dialog(enable = (FFtype == ThermoPower.Choices.Flow1D.FFtypes.Colebrook) or useFrictionModel));
       parameter Real Kfc=1 "Friction factor correction coefficient";
       parameter Boolean DynamicMomentum=false
         "Inertial phenomena accounted for"
@@ -2905,6 +2920,15 @@ The packages Medium are redeclared and a mass balance determines the composition
       parameter Boolean noInitialPressure=false
         "Remove initial equation on pressure"
         annotation (Dialog(tab="Initialisation"),choices(checkBox=true));
+      parameter Boolean useFrictionModel=false
+        "Use replaceable model for friction factor"
+        annotation (choices(checkBox=true));
+
+      replaceable function Cfcorr =
+          ThermoPower.Functions.FrictionFactors.NoFriction constrainedby
+        ThermoPower.Functions.FrictionFactors.BaseFFcorr
+        "External Fanning friction factor correlation"
+        annotation(choicesAllMatching = true, Dialog(enable = (FFtype == ThermoPower.Choices.Flow1D.FFtypes.External)));
 
       function squareReg = ThermoPower.Functions.squareReg;
     protected
@@ -3793,6 +3817,7 @@ Several functions are provided in the package <tt>Functions.FanCharacteristics</
 </ul>
 </html>"));
     end FanBase;
+
   end BaseClasses;
 
   model SourceP "Pressure source for gas flows"
@@ -4553,6 +4578,17 @@ Several functions are provided in the package <tt>Functions.FanCharacteristics</
 </html>"),
       DymolaStoredErrors);
   end Flow1D;
+
+  model Fan "Model for a gas fan, with constant speed"
+    extends ThermoPower.Gas.BaseClasses.FanBase(
+      redeclare replaceable package Medium = DEMO.Media.Helium.DefaultHelium
+        constrainedby DEMO.Internals.HeliumMedium);
+
+  equation
+    n=n0;
+
+  end Fan;
+
   annotation (Documentation(info="<HTML>
 This package contains models of physical processes and components using ideal gases as working fluid.
 <p>All models with dynamic equations provide initialisation support. Set the <tt>initOpt</tt> parameter to the appropriate value: