From 506d32c468e065653acb469c5122fb9642dff166 Mon Sep 17 00:00:00 2001 From: AntoineGautier Date: Mon, 14 Oct 2024 17:19:03 +0200 Subject: [PATCH] Implement spline interpolation --- .../BasicFlowFunction_dp_DerivativeCheck.mo | 1 + .../FlowModels/basicFlowFunction_dp_der.mo | 2 +- .../FlowModels/basicFlowFunction_dp_der2.mo | 2 +- .../basicFlowFunction_m_flow_der.mo | 2 +- .../basicFlowFunction_m_flow_der2.mo | 2 +- IBPSA/Fluid/FixedResistances/CheckValve.mo | 99 ++++++++++++------- .../FixedResistances/Examples/CheckValve.mos | 1 - 7 files changed, 68 insertions(+), 41 deletions(-) diff --git a/IBPSA/Fluid/BaseClasses/FlowModels/Validation/BasicFlowFunction_dp_DerivativeCheck.mo b/IBPSA/Fluid/BaseClasses/FlowModels/Validation/BasicFlowFunction_dp_DerivativeCheck.mo index fdcc56b385..e3f1c0903c 100644 --- 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 index 335c0fa6cb..0777dbb0ab 100644 --- a/IBPSA/Fluid/BaseClasses/FlowModels/basicFlowFunction_dp_der.mo +++ b/IBPSA/Fluid/BaseClasses/FlowModels/basicFlowFunction_dp_der.mo @@ -31,7 +31,7 @@ Documentation(info=" Function that implements the first order derivative of IBPSA.Fluid.BaseClasses.FlowModels.basicFlowFunction_dp -with respect to the mass flow rate. +with respect to time.

", revisions=" diff --git a/IBPSA/Fluid/BaseClasses/FlowModels/basicFlowFunction_dp_der2.mo b/IBPSA/Fluid/BaseClasses/FlowModels/basicFlowFunction_dp_der2.mo index e64a9fbdb6..cd405fe2c2 100644 --- a/IBPSA/Fluid/BaseClasses/FlowModels/basicFlowFunction_dp_der2.mo +++ b/IBPSA/Fluid/BaseClasses/FlowModels/basicFlowFunction_dp_der2.mo @@ -36,7 +36,7 @@ Documentation(info=" Function that implements the second order derivative of IBPSA.Fluid.BaseClasses.FlowModels.basicFlowFunction_dp -with respect to the mass flow rate. +with respect to time.

", revisions=" diff --git a/IBPSA/Fluid/BaseClasses/FlowModels/basicFlowFunction_m_flow_der.mo b/IBPSA/Fluid/BaseClasses/FlowModels/basicFlowFunction_m_flow_der.mo index 05c27661a2..7169c2a5c4 100644 --- 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=" Function that implements the first order derivative of IBPSA.Fluid.BaseClasses.FlowModels.basicFlowFunction_m_flow -with respect to the mass flow rate. +with respect to time.

", revisions=" diff --git a/IBPSA/Fluid/BaseClasses/FlowModels/basicFlowFunction_m_flow_der2.mo b/IBPSA/Fluid/BaseClasses/FlowModels/basicFlowFunction_m_flow_der2.mo index 6c70c3ef0e..2b63cf1a76 100644 --- 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=" Function that implements the second order derivative of IBPSA.Fluid.BaseClasses.FlowModels.basicFlowFunction_m_flow -with respect to the mass flow rate. +with respect to time.

", revisions=" diff --git a/IBPSA/Fluid/FixedResistances/CheckValve.mo b/IBPSA/Fluid/FixedResistances/CheckValve.mo index 883f7acc25..6e5196f446 100644 --- 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 index 30e334f42f..20327d0c8d 100644 --- 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}});