Implement spline interpolation

ibpsa · Oct 14, 2024 · 506d32c · 506d32c
1 parent 5136ecf
commit 506d32c
Show file tree

Hide file tree

Showing 7 changed files with 68 additions and 41 deletions.
diff --git a/IBPSA/Fluid/BaseClasses/FlowModels/Validation/BasicFlowFunction_dp_DerivativeCheck.mo b/IBPSA/Fluid/BaseClasses/FlowModels/Validation/BasicFlowFunction_dp_DerivativeCheck.mo
@@ -21,6 +21,7 @@ equation
     k = k,
     m_flow_turbulent=m_flow_turbulent);
   der(m_flow) = der(m_flow_comp);
+  der_m_flow_dp = IBPSA.Fluid.BaseClasses.FlowModels.basicFlowFunction_dp_der(dp);
   err = m_flow-m_flow_comp;
   assert(abs(err) < 1E-3, "Error in implementation.");
 annotation (

diff --git a/IBPSA/Fluid/BaseClasses/FlowModels/basicFlowFunction_dp_der.mo b/IBPSA/Fluid/BaseClasses/FlowModels/basicFlowFunction_dp_der.mo
@@ -31,7 +31,7 @@ Documentation(info="<html>
 Function that implements the first order derivative of
 <a href=\"modelica://IBPSA.Fluid.BaseClasses.FlowModels.basicFlowFunction_dp\">
 IBPSA.Fluid.BaseClasses.FlowModels.basicFlowFunction_dp</a>
-with respect to the mass flow rate.
+with respect to time.
 </p>
 </html>",
 revisions="<html>

diff --git a/IBPSA/Fluid/BaseClasses/FlowModels/basicFlowFunction_dp_der2.mo b/IBPSA/Fluid/BaseClasses/FlowModels/basicFlowFunction_dp_der2.mo
@@ -36,7 +36,7 @@ Documentation(info="<html>
 Function that implements the second order derivative of
 <a href=\"modelica://IBPSA.Fluid.BaseClasses.FlowModels.basicFlowFunction_dp\">
 IBPSA.Fluid.BaseClasses.FlowModels.basicFlowFunction_dp</a>
-with respect to the mass flow rate.
+with respect to time.
 </p>
 </html>",
 revisions="<html>

diff --git a/IBPSA/Fluid/BaseClasses/FlowModels/basicFlowFunction_m_flow_der.mo b/IBPSA/Fluid/BaseClasses/FlowModels/basicFlowFunction_m_flow_der.mo
@@ -32,7 +32,7 @@ Documentation(info="<html>
 Function that implements the first order derivative of
 <a href=\"modelica://IBPSA.Fluid.BaseClasses.FlowModels.basicFlowFunction_m_flow\">
 IBPSA.Fluid.BaseClasses.FlowModels.basicFlowFunction_m_flow</a>
-with respect to the mass flow rate.
+with respect to time.
 </p>
 </html>",
 revisions="<html>

diff --git a/IBPSA/Fluid/BaseClasses/FlowModels/basicFlowFunction_m_flow_der2.mo b/IBPSA/Fluid/BaseClasses/FlowModels/basicFlowFunction_m_flow_der2.mo
@@ -36,7 +36,7 @@ Documentation(info="<html>
 Function that implements the second order derivative of
 <a href=\"modelica://IBPSA.Fluid.BaseClasses.FlowModels.basicFlowFunction_m_flow\">
 IBPSA.Fluid.BaseClasses.FlowModels.basicFlowFunction_m_flow</a>
-with respect to the mass flow rate.
+with respect to time.
 </p>
 </html>",
 revisions="<html>

diff --git a/IBPSA/Fluid/FixedResistances/CheckValve.mo b/IBPSA/Fluid/FixedResistances/CheckValve.mo
@@ -27,48 +27,75 @@ model CheckValve "Check valve that avoids flow reversal"
     else 0
     "Flow coefficient of fixed resistance that may be in series with valve,
     k=m_flow/sqrt(dp), with unit=(kg.m)^(1/2).";
-
-  Real k(min=Modelica.Constants.small)
-    "Flow coefficient of valve and pipe in series in allowed/forward direction,
-    k=m_flow/sqrt(dp), with unit=(kg.m)^(1/2).";
-
 protected
-  Real a
-    "Scaled pressure variable";
-  Real cv
-    "Twice differentiable Heaviside check valve characteristic";
-  Real kCv
-    "Smoothed restriction characteristic";
-
+  parameter Real k_min=if dpFixed_nominal > Modelica.Constants.eps then
+    sqrt(1 / (1 / kFixed ^ 2 + 1 /(l * Kv_SI) ^ 2)) else l * Kv_SI
+    "Minimum flow coefficient (valve closed)";
+  parameter Real k_max=if dpFixed_nominal > Modelica.Constants.eps then
+    sqrt(1 / (1 / kFixed ^ 2 + 1 / Kv_SI ^ 2)) else Kv_SI
+    "Maximum flow coefficient (valve fully open)";
+  Modelica.Units.SI.MassFlowRate m_flow_min
+    "Flow rate through closed valve";
+  Modelica.Units.SI.MassFlowRate m_flow_max
+    "Flow rate through fully open valve";
+  Modelica.Units.SI.MassFlowRate m_flow_smooth
+    "Smooth interpolation result between two flow regimes";
 initial equation
-  assert(dpFixed_nominal > -Modelica.Constants.eps,
-    "In " + getInstanceName() + ": We require dpFixed_nominal >= 0.
+  assert(dpFixed_nominal > - Modelica.Constants.eps, "In " + getInstanceName() +
+    ": We require dpFixed_nominal >= 0.
     Received dpFixed_nominal = " + String(dpFixed_nominal) + " Pa.");
-  assert(l > -Modelica.Constants.eps,
-    "In " + getInstanceName() + ": We require l >= 0. Received l = " + String(l));
+  assert(l > - Modelica.Constants.eps, "In " + getInstanceName() +
+    ": We require l >= 0. Received l = " + String(l));
 equation
-  a = dp/dpValve_closing;
-  cv = smooth(2, max(0, min(1, a^3*(10+a*(-15+6*a)))));
-  kCv = Kv_SI*(cv*(1-l) + l);
-
-  if (dpFixed_nominal > Modelica.Constants.eps) then
-    k = sqrt(1/(1/kFixed^2 + 1/kCv^2));
-  else
-    k = kCv;
-  end if;
-
+  // min function ensures that m_flow_y1 does not increase further for dp_x > dp_x1
+  m_flow_min=IBPSA.Fluid.BaseClasses.FlowModels.basicFlowFunction_dp(
+    dp=min(dp, 0),
+    k=k_min,
+    m_flow_turbulent=m_flow_turbulent);
+  // max function ensures that m_flow_y2 does not decrease further for dp_x < dp_x2
+  m_flow_max=IBPSA.Fluid.BaseClasses.FlowModels.basicFlowFunction_dp(
+    dp=max(dp, dpValve_closing),
+    k=k_max,
+    m_flow_turbulent=m_flow_turbulent);
+  // We compute d(m_flow)/dp by calling d(m_flow)/dt with dp/dt=1,
+  // and d2(m_flow)/dp2 by calling d2(m_flow)/dt2 with dp/dt=1 and d2p/dt2=0
+  m_flow_smooth=noEvent(smooth(2,
+    if dp <= 0 then m_flow_min
+    elseif dp >= dpValve_closing then m_flow_max
+    else IBPSA.Utilities.Math.Functions.quinticHermite(
+      x=dp,
+      x1=0,
+      x2=dpValve_closing,
+      y1=m_flow_min,
+      y2=m_flow_max,
+      y1d=IBPSA.Fluid.BaseClasses.FlowModels.basicFlowFunction_dp_der(
+        dp=0,
+        k=k_min,
+        m_flow_turbulent=m_flow_turbulent,
+        dp_der=1),
+      y2d=IBPSA.Fluid.BaseClasses.FlowModels.basicFlowFunction_dp_der(
+        dp=dpValve_closing,
+        k=k_max,
+        m_flow_turbulent=m_flow_turbulent,
+        dp_der=1),
+      y1dd=IBPSA.Fluid.BaseClasses.FlowModels.basicFlowFunction_dp_der2(
+        dp=0,
+        k=k_min,
+        m_flow_turbulent=m_flow_turbulent,
+        dp_der=1,
+        dp_der2=0),
+      y2dd=IBPSA.Fluid.BaseClasses.FlowModels.basicFlowFunction_dp_der2(
+        dp=dpValve_closing,
+        k=k_max,
+        m_flow_turbulent=m_flow_turbulent,
+        dp_der=1,
+        dp_der2=0))));
   if homotopyInitialization then
-    m_flow = homotopy(
-      actual=IBPSA.Fluid.BaseClasses.FlowModels.basicFlowFunction_dp(
-            dp = dp,
-            k = k,
-            m_flow_turbulent = m_flow_turbulent),
-      simplified = m_flow_nominal_pos*dp/dp_nominal_pos);
+    m_flow=homotopy(
+      actual=m_flow_smooth,
+      simplified=m_flow_nominal_pos * dp / dp_nominal_pos);
   else
-    m_flow = IBPSA.Fluid.BaseClasses.FlowModels.basicFlowFunction_dp(
-        dp = dp,
-        k = k,
-        m_flow_turbulent = m_flow_turbulent);
+    m_flow=m_flow_smooth;
   end if;
   annotation (Icon(coordinateSystem(preserveAspectRatio=true, extent={{-100,-100},
             {100,100}}), graphics={

diff --git a/IBPSA/Resources/Scripts/Dymola/Fluid/FixedResistances/Examples/CheckValve.mos b/IBPSA/Resources/Scripts/Dymola/Fluid/FixedResistances/Examples/CheckValve.mos
@@ -1,4 +1,3 @@
 simulateModel("IBPSA.Fluid.FixedResistances.Examples.CheckValve", method="Cvode", tolerance=1e-06, stopTime=1, resultFile="CheckValve");
 createPlot(id=1, position={0, 0, 988, 599}, y={"checkValve.dp", "checkValve_m_flow.dp"}, range={0.0, 1.0, -20000.0, 40000.0}, grid=true, colors={{28,108,200}, {238,46,47}});
 createPlot(id=1, position={0, 0, 988, 197}, y={"checkValve.m_flow", "checkValveDpFix.m_flow", "checkValve_m_flow.m_flow"}, range={0.0, 1.0, -2.0, 6.0}, grid=true, subPlot=2, colors={{28,108,200}, {238,46,47}, {0,140,72}});
-createPlot(id=1, position={0, 0, 988, 197}, y={"checkValve.k", "checkValveDpFix.k", "checkValve_m_flow.k"}, range={0.0, 1.0, -0.02, 0.06}, grid=true, subPlot=3, colors={{28,108,200}, {238,46,47}, {0,140,72}});