MSL中的動態管道模型，有限體積法

Question

我正在嘗試使用Modelica對彈性管道組成的系統進行建模。 現在，我試圖使用與Modelica.Fluid庫中相同的方法（有限體積，交錯）來實現我自己的動態管道模型（剛性，還沒有彈性），但當然不包括所有選項。

該模型應該更易於理解，因爲它是一個平面模型，不能從其他類擴展。這很重要，因爲即使沒有Modelica Knowhow，我的同事也可以理解這個模型，我可以說服他們Modelica是適合我們用途的適當工具！

作爲一個測試案例，我使用帶有階躍信號（waterhammer）的質量流量源。 我的模型給出的結果與Modelica.Fluid組件不同。 我真的很感激，如果有人能幫助我，理解發生了什麼！

的測試系統看起來是這樣的：

的結果11個細胞是這樣的：

正如你所看到的，壓力峯值對於MSL組件而言更高，並且頻率/週期不相同。當我選擇更多的單元格時，錯誤會變小。

我很確定我使用的是完全相同的方程。 它可能是數字原因的原因（我嘗試使用名義值）？ 我還爲Modelica.Fluid組件包含了我自己的「固定zeta」流動模型，以便在固定壓力損失係數zeta的情況下進行比較。

我管模式的代碼很短，如果我得到它的這樣的工作，它會是非常好的：

model Pipe_FVM_staggered 

    // Import 
    import SI = Modelica.SIunits; 
    import Modelica.Constants.pi; 

    // Medium 
    replaceable package Medium = Modelica.Media.Interfaces.PartialMedium "Medium in the component" 
    annotation (choicesAllMatching = true); 

    // Interfaces, Ports 
    Modelica.Fluid.Interfaces.FluidPort_a port_a(redeclare package Medium = Medium) annotation (Placement(transformation(extent={{-110,-10},{-90,10}}))); 
    Modelica.Fluid.Interfaces.FluidPort_b port_b(redeclare package Medium = Medium) annotation (Placement(transformation(extent={{90,-10},{110,10}}))); 

    // Parameters 
    parameter Integer n(min=2) = 3 "Number of cells"; // No effect yet, only for icon 
    parameter SI.Length L = 1 "Length"; 
    parameter SI.Diameter D = 0.010 "Diameter"; 
    parameter SI.Height R = 2.5e-5 "Roughness"; 
    parameter Boolean use_fixed_zeta = false "Use fixed zeta value instead of Moody chart"; 
    parameter SI.CoefficientOfFriction zeta = 1; 

    // Initialization 
    parameter Medium.Temperature T_start = 293.15 "Start temperature" annotation(Dialog(tab="Initialization")); 
    parameter Medium.MassFlowRate mflow_start = 1 "Start mass flow rate in design direction" annotation(Dialog(tab="Initialization")); 
    parameter Medium.AbsolutePressure p_a_start = 2e5 "Start pressure p[1] at design inflow" annotation(Dialog(tab="Initialization")); 
    parameter Medium.AbsolutePressure p_b_start = 1e5 "Start pressure for p[n+1] at design outflow" annotation(Dialog(tab="Initialization")); 
    // parameter Medium.AbsolutePressure p_start = (p_a_start + p_b_start)/2 annotation(Dialog(tab="Initialization")); 
    parameter Medium.AbsolutePressure p_start[:] = linspace(p_a_start, p_b_start, n) annotation(Dialog(tab="Initialization")); 
    // parameter Medium.SpecificEnthalpy h_start[:] = Medium.specificEnthalpy_pTX(p_start, T_start, Medium.X_default); 
    parameter Medium.SpecificEnthalpy h_start = Medium.specificEnthalpy_pTX((p_a_start + p_b_start)/2, T_start, Medium.X_default) annotation(Dialog(tab="Initialization")); 
    parameter SI.AbsolutePressure dp_nominal = 1e5; 
    parameter SI.MassFlowRate m_flow_nominal = 1; 

    // Variables general 
    SI.Length dL = L/n; 
    SI.Area A(nominal=0.001) = D^2*pi/4; 
    SI.Volume V = A * dL; 

    // Variables cell centers: positiv in direction a -> b 
    Medium.AbsolutePressure p[n](start = p_start, each stateSelect=StateSelect.prefer) annotation(Dialog(tab="Initialization", showStartAttribute=true, enable=false)); 
    Medium.SpecificEnthalpy h[n](each start = h_start, each stateSelect=StateSelect.prefer) annotation(Dialog(tab="Initialization", showStartAttribute=true, enable=false)); 
    Medium.ThermodynamicState state[n] = Medium.setState_phX(p,h); 
    SI.Mass m[n] = rho .* V; 
    Medium.Density rho[n] = Medium.density(state); 
    SI.InternalEnergy U[n] = m .* u; 
    Medium.SpecificInternalEnergy u[n] = Medium.specificInternalEnergy(state); 
    Medium.Temperature T[n] = Medium.temperature(state); 
    Medium.DynamicViscosity mu[n] = Medium.dynamicViscosity(state); 
    SI.Velocity v[n](nominal=0.2) = 0.5 * (mflow[1:n] + mflow[2:n+1]) ./ rho ./ A; 
    SI.Power Wflow[n]; 
    SI.MomentumFlux Iflow[n] = v .* v .* rho * A; 

    // Variables faces: positiv in direction a -> b 
    Medium.MassFlowRate mflow[n+1](each start = mflow_start, each stateSelect=StateSelect.prefer, nominal=0.25) annotation(Dialog(tab="Initialization", showStartAttribute=true, enable=false)); 
    Medium.EnthalpyFlowRate Hflow[n+1]; 
    SI.Momentum I[n-1] = mflow[2:n] * dL; 
    SI.Force Fp[n-1]; 
    SI.Force Ff[n-1]; 
    SI.PressureDifference dpf[n-1](each start = (p_a_start - p_b_start)/(n-1), nominal=0.01e5) annotation(Dialog(tab="Initialization", showStartAttribute=true, enable=false)); 

equation 

    der(m) = mflow[1:n] - mflow[2:n+1];     // Mass balance 
    der(U) = Hflow[1:n] - Hflow[2:n+1] + Wflow;   // Energy balance 
    der(I) = Iflow[1:n-1] - Iflow[2:n] + Fp - Ff;   // Momentum balance, staggered 

    Hflow[1] = semiLinear(mflow[1], inStream(port_a.h_outflow), h[1]); 
    Hflow[2:n] = semiLinear(mflow[2:n], h[1:n-1], h[2:n]); 
    Hflow[n+1] = semiLinear(mflow[n+1], h[n], inStream(port_b.h_outflow)); 

    Wflow[1] = v[1] * A .* ((p[2] - p[1])/2 + dpf[1]/2); 
    Wflow[2:n-1] = v[2:n-1] * A .* ((p[3:n]-p[1:n-2])/2 + (dpf[1:n-2]+dpf[2:n-1])/2); 
    Wflow[n] = v[n] * A .* ((p[n] - p[n-1])/2 + dpf[n-1]/2); 

    Fp = A * (p[1:n-1] - p[2:n]); 
    Ff = A * dpf; // dpf = Ff ./ A; 

    if use_fixed_zeta then 
    dpf = 1/2 * zeta/(n-1) * (mflow[2:n]).^2 ./ (0.5*(rho[1:n-1] + rho[2:n]) * A * A); 
    else 
    dpf = homotopy(
     actual = Modelica.Fluid.Pipes.BaseClasses.WallFriction.Detailed.pressureLoss_m_flow(
     m_flow = mflow[2:n], 
     rho_a = rho[1:n-1], 
     rho_b = rho[2:n], 
     mu_a = mu[1:n-1], 
     mu_b = mu[2:n], 
     length = dL, 
     diameter = D, 
     roughness = R, 
     m_flow_small = 0.001), 
     simplified = dp_nominal/(n-1)/m_flow_nominal*mflow[2:n]); 
    end if; 

    // Boundary conditions 
    mflow[1] = port_a.m_flow; 
    mflow[n] = -port_b.m_flow; 
    p[1] = port_a.p; 
    p[n] = port_b.p; 
    port_a.h_outflow = h[1]; 
    port_b.h_outflow = h[n]; 

initial equation 
    der(mflow[2:n]) = zeros(n-1); 
    der(p) = zeros(n); 
    der(h) = zeros(n); 

    annotation (Icon(coordinateSystem(preserveAspectRatio=false), graphics={Rectangle(
      extent={{-100,60},{100,-60}}, 
      fillColor={255,255,255}, 
      fillPattern=FillPattern.HorizontalCylinder, 
      lineColor={0,0,0}), 
     Line(
      points={{-100,60},{-100,-60}}, 
      color={0,0,0}, 
      thickness=0.5), 
     Line(
      points={{-60,60},{-60,-60}}, 
      color={0,0,0}, 
      thickness=0.5), 
     Line(
      points={{-20,60},{-20,-60}}, 
      color={0,0,0}, 
      thickness=0.5), 
     Line(
      points={{20,60},{20,-60}}, 
      color={0,0,0}, 
      thickness=0.5), 
     Line(
      points={{60,60},{60,-60}}, 
      color={0,0,0}, 
      thickness=0.5), 
     Line(
      points={{100,60},{100,-60}}, 
      color={0,0,0}, 
      thickness=0.5), 
     Line(
      points={{60,-80},{-60,-80}}, 
      color={0,128,255}, 
      visible=showDesignFlowDirection), 
     Polygon(
      points={{20,-65},{60,-80},{20,-95},{20,-65}}, 
      lineColor={0,128,255}, 
      fillColor={0,128,255}, 
      fillPattern=FillPattern.Solid, 
      visible=showDesignFlowDirection), 
     Text(
      extent={{-150,100},{150,60}}, 
      lineColor={0,0,255}, 
      textString="%name"), 
     Text(
      extent={{-40,22},{40,-18}}, 
      lineColor={0,0,0}, 
      textString="n = %n")}),        Diagram(
     coordinateSystem(preserveAspectRatio=false))); 
end Pipe_FVM_staggered; 

我這個問題，因爲相當長一段時間掙扎，所以任何意見或提示都非常感謝！ 如果您需要更多信息或測試結果，請告訴我！

這是試驗例的代碼：

變焦峯1和2：：

model Test_Waterhammer 

    extends Modelica.Icons.Example; 
    import SI = Modelica.SIunits; 
    import g = Modelica.Constants.g_n; 

    replaceable package Medium = Modelica.Media.Water.StandardWater; 

    Modelica.Fluid.Sources.Boundary_pT outlet(
    redeclare package Medium = Medium, 
    nPorts=1, 
    p=2000000, 
    T=293.15) 
    annotation (Placement(transformation(extent={{90,-10},{70,10}}))); 

    inner Modelica.Fluid.System system(
    allowFlowReversal=true, 
    energyDynamics=Modelica.Fluid.Types.Dynamics.SteadyStateInitial, 
    massDynamics=Modelica.Fluid.Types.Dynamics.SteadyStateInitial, 
    momentumDynamics=Modelica.Fluid.Types.Dynamics.SteadyStateInitial, 
    m_flow_start=0.1, 
    m_flow_small=0.0001) 
    annotation (Placement(transformation(extent={{60,60},{80,80}}))); 

    Modelica.Fluid.Sources.MassFlowSource_T inlet(
    redeclare package Medium = Medium, 
    nPorts=1, 
    m_flow=0.1, 
    use_m_flow_in=true, 
    T=293.15) 
    annotation (Placement(transformation(extent={{-50,-10},{-30,10}}))); 

    Modelica.Blocks.Sources.TimeTable timeTable(table=[0,0.1; 1,0.1; 1,0.25; 
     40,0.25; 40,0.35; 60,0.35]) 
    annotation (Placement(transformation(extent={{-90,10},{-70,30}}))); 

    Pipe_FVM_staggered          pipe(
    redeclare package Medium = Medium, 
    R=0.035*0.005, 
    mflow_start=0.1, 
    L=1000, 
    m_flow_nominal=0.1, 
    D=0.035, 
    zeta=2000, 
    n=11, 
    use_fixed_zeta=false, 
    T_start=293.15, 
    p_a_start=2010000, 
    p_b_start=2000000, 
    dp_nominal=10000) 
    annotation (Placement(transformation(extent={{10,-10},{30,10}}))); 

    Modelica.Fluid.Pipes.DynamicPipe pipeMSL(
    redeclare package Medium = Medium, 
    allowFlowReversal=true, 
    length=1000, 
    roughness=0.035*0.005, 
    m_flow_start=0.1, 
    energyDynamics=Modelica.Fluid.Types.Dynamics.SteadyStateInitial, 
    massDynamics=Modelica.Fluid.Types.Dynamics.SteadyStateInitial, 
    momentumDynamics=Modelica.Fluid.Types.Dynamics.SteadyStateInitial, 
    diameter=0.035, 
    modelStructure=Modelica.Fluid.Types.ModelStructure.av_vb, 
    redeclare model FlowModel = 
     Modelica.Fluid.Pipes.BaseClasses.FlowModels.DetailedPipeFlow (
      useUpstreamScheme=false, use_Ib_flows=true), 
    p_a_start=2010000, 
    p_b_start=2000000, 
    T_start=293.15, 
    nNodes=11) 
    annotation (Placement(transformation(extent={{10,-50},{30,-30}}))); 

    Modelica.Fluid.Sources.MassFlowSource_T inlet1(
    redeclare package Medium = Medium, 
    nPorts=1, 
    m_flow=0.1, 
    use_m_flow_in=true, 
    T=293.15) 
    annotation (Placement(transformation(extent={{-48,-50},{-28,-30}}))); 

    Modelica.Fluid.Sources.Boundary_pT outlet1(
    redeclare package Medium = Medium, 
    nPorts=1, 
    p=2000000, 
    T=293.15) 
    annotation (Placement(transformation(extent={{90,-50},{70,-30}}))); 


equation 
    connect(inlet.ports[1], pipe.port_a) 
    annotation (Line(points={{-30,0},{-10,0},{10,0}}, color={0,127,255})); 
    connect(pipe.port_b, outlet.ports[1]) 
    annotation (Line(points={{30,0},{50,0},{70,0}}, color={0,127,255})); 
    connect(inlet1.ports[1], pipeMSL.port_a) 
    annotation (Line(points={{-28,-40},{-10,-40},{10,-40}}, color={0,127,255})); 
    connect(pipeMSL.port_b, outlet1.ports[1]) 
    annotation (Line(points={{30,-40},{50,-40},{70,-40}}, color={0,127,255})); 
    connect(timeTable.y, inlet.m_flow_in) 
    annotation (Line(points={{-69,20},{-60,20},{-60,8},{-50,8}}, color={0,0,127})); 
    connect(inlet1.m_flow_in, inlet.m_flow_in) 
    annotation (Line(points={{-48,-32},{-60,-32},{-60,8},{-50,8}}, color={0,0,127})); 


    annotation (Icon(coordinateSystem(preserveAspectRatio=false)), Diagram(
     coordinateSystem(preserveAspectRatio=false)), 
    experiment(
     StopTime=15, 
     __Dymola_NumberOfIntervals=6000, 
     Tolerance=1e-005, 
     __Dymola_Algorithm="Dassl")); 

end Test_Waterhammer; 

我已經與301個細胞運行測試

解決方法：修改由scottG

model FVM_staggered_Ncells 

    // Import 
    import SI = Modelica.SIunits; 
    import Modelica.Constants.pi; 

    // Medium 
    replaceable package Medium = Modelica.Media.Interfaces.PartialMedium "Medium in the component" 
    annotation (choicesAllMatching = true); 

    // Interfaces, Ports 
    Modelica.Fluid.Interfaces.FluidPort_a port_a(redeclare package Medium = Medium) annotation (Placement(transformation(extent={{-110,-10},{-90,10}}))); 
    Modelica.Fluid.Interfaces.FluidPort_b port_b(redeclare package Medium = Medium) annotation (Placement(transformation(extent={{90,-10},{110,10}}))); 

    // Parameters 
    parameter Integer n(min=2) = 3 "Number of cells"; // No effect yet, only for icon 
    parameter SI.Length L = 1 "Length"; 
    parameter SI.Diameter D = 0.010 "Diameter"; 
    parameter SI.Height R = 2.5e-5 "Roughness"; 
    parameter Boolean use_fixed_zeta = false "Use fixed zeta value instead of Moody chart"; 
    parameter SI.CoefficientOfFriction zeta = 1; 

    // Initialization 
    parameter Medium.Temperature T_start = 293.15 "Start temperature" annotation(Dialog(tab="Initialization")); 
    parameter Medium.MassFlowRate mflow_start = 1 "Start mass flow rate in design direction" annotation(Dialog(tab="Initialization")); 
    parameter Medium.AbsolutePressure p_a_start = 2e5 "Start pressure p[1] at design inflow" annotation(Dialog(tab="Initialization")); 
    parameter Medium.AbsolutePressure p_b_start = 1e5 "Start pressure for p[n+1] at design outflow" annotation(Dialog(tab="Initialization")); 
    parameter Medium.AbsolutePressure p_start[:] = linspace(p_a_start, p_b_start, n) annotation(Dialog(tab="Initialization")); 
    // parameter Medium.SpecificEnthalpy h_start[:] = Medium.specificEnthalpy_pTX(p_start, T_start, Medium.X_default); 
    parameter Medium.SpecificEnthalpy h_start = Medium.specificEnthalpy_pTX((p_a_start + p_b_start)/2, T_start, Medium.X_default) annotation(Dialog(tab="Initialization")); 
    parameter SI.AbsolutePressure dp_nominal = 1e5; 
    parameter SI.MassFlowRate m_flow_nominal = 1; 

    // Variables general 
    SI.Length dL = L/n; 
    SI.Length dLs[n-1] = cat(1,{1.5*dL}, fill(dL,n-3), {1.5*dL}); 
    SI.Area A = D^2*pi/4; 
    SI.Volume V = A * dL; 

    // Variables cell centers: positiv in direction a -> b 
    Medium.AbsolutePressure p[n](start = p_start, each stateSelect=StateSelect.prefer) annotation(Dialog(tab="Initialization", showStartAttribute=true, enable=false)); 
    Medium.SpecificEnthalpy h[n](each start = h_start, each stateSelect=StateSelect.prefer) annotation(Dialog(tab="Initialization", showStartAttribute=true, enable=false)); 
    Medium.ThermodynamicState state[n] = Medium.setState_phX(p,h); 
    SI.Mass m[n] = rho .* V; 
    Medium.Density rho[n] = Medium.density(state); 
    SI.InternalEnergy U[n] = m .* u; 
    Medium.SpecificInternalEnergy u[n] = Medium.specificInternalEnergy(state); 
    Medium.Temperature T[n] = Medium.temperature(state); 
    Medium.DynamicViscosity mu[n] = Medium.dynamicViscosity(state); 
    SI.Velocity v[n] = 0.5 * (mflow[1:n] + mflow[2:n+1]) ./ rho ./ A; 
    SI.Power Wflow[n]; 
    SI.MomentumFlux Iflow[n] = v .* v .* rho * A; 

    // Variables faces: positiv in direction a -> b 
    Medium.MassFlowRate mflow[n+1](each start = mflow_start, each stateSelect=StateSelect.prefer) annotation(Dialog(tab="Initialization", showStartAttribute=true, enable=false)); 
    Medium.EnthalpyFlowRate Hflow[n+1]; 
    SI.Momentum I[n-1] = mflow[2:n] .* dLs; 
    SI.Force Fp[n-1]; 
    SI.Force Ff[n-1]; 
    SI.PressureDifference dpf[n-1](each start = (p_a_start - p_b_start)/(n-1)) annotation(Dialog(tab="Initialization", showStartAttribute=true, enable=false)); 


equation 

    der(m) = mflow[1:n] - mflow[2:n+1];     // Mass balance 
    der(U) = Hflow[1:n] - Hflow[2:n+1] + Wflow;   // Energy balance 
    der(I) = Iflow[1:n-1] - Iflow[2:n] + Fp - Ff;   // Momentum balance, staggered 

    Hflow[1] = semiLinear(mflow[1], inStream(port_a.h_outflow), h[1]); 
    Hflow[2:n] = semiLinear(mflow[2:n], h[1:n-1], h[2:n]); 
    Hflow[n+1] = semiLinear(mflow[n+1], h[n], inStream(port_b.h_outflow)); 

    Wflow[1] = v[1] * A .* ((p[2] - p[1])/2 + dpf[1]/2); 
    Wflow[2:n-1] = v[2:n-1] * A .* ((p[3:n]-p[1:n-2])/2 + (dpf[1:n-2]+dpf[2:n-1])/2); 
    Wflow[n] = v[n] * A .* ((p[n] - p[n-1])/2 + dpf[n-1]/2); 

    Fp = A * (p[1:n-1] - p[2:n]); 
    Ff = A * dpf; 

    if use_fixed_zeta then 
    dpf = 0.5 * zeta/(n-1) * abs(mflow[2:n]) .* mflow[2:n] ./ (0.5*(rho[1:n-1] + rho[2:n]) * A * A); 
    else 
    dpf = homotopy(
     actual = Modelica.Fluid.Pipes.BaseClasses.WallFriction.Detailed.pressureLoss_m_flow(
     m_flow = mflow[2:n], 
     rho_a = 0.5*(rho[1:n-1] + rho[2:n]), 
     rho_b = 0.5*(rho[1:n-1] + rho[2:n]), 
     mu_a = 0.5*(mu[1:n-1] + mu[2:n]), 
     mu_b = 0.5*(mu[1:n-1] + mu[2:n]), 
     length = dLs, 
     diameter = D, 
     roughness = R, 
     m_flow_small = 0.001), 
     simplified = dp_nominal/(n-1)/m_flow_nominal*mflow[2:n]); 
    end if; 

    // Boundary conditions 
    mflow[1] = port_a.m_flow; 
    mflow[n+1] = -port_b.m_flow; 
    p[1] = port_a.p; 
    p[n] = port_b.p; 
    port_a.h_outflow = h[1]; 
    port_b.h_outflow = h[n]; 

initial equation 
    der(mflow[2:n]) = zeros(n-1); 
    der(p) = zeros(n); 
    der(h) = zeros(n); 

    annotation (Icon(coordinateSystem(preserveAspectRatio=false), graphics={Rectangle(
      extent={{-100,60},{100,-60}}, 
      fillColor={255,255,255}, 
      fillPattern=FillPattern.HorizontalCylinder, 
      lineColor={0,0,0}), 
     Line(
      points={{-100,60},{-100,-60}}, 
      color={0,0,0}, 
      thickness=0.5), 
     Line(
      points={{-60,60},{-60,-60}}, 
      color={0,0,0}, 
      thickness=0.5), 
     Line(
      points={{-20,60},{-20,-60}}, 
      color={0,0,0}, 
      thickness=0.5), 
     Line(
      points={{20,60},{20,-60}}, 
      color={0,0,0}, 
      thickness=0.5), 
     Line(
      points={{60,60},{60,-60}}, 
      color={0,0,0}, 
      thickness=0.5), 
     Line(
      points={{100,60},{100,-60}}, 
      color={0,0,0}, 
      thickness=0.5), 
     Line(
      points={{60,-80},{-60,-80}}, 
      color={0,128,255}, 
      visible=showDesignFlowDirection), 
     Polygon(
      points={{20,-65},{60,-80},{20,-95},{20,-65}}, 
      lineColor={0,128,255}, 
      fillColor={0,128,255}, 
      fillPattern=FillPattern.Solid, 
      visible=showDesignFlowDirection), 
     Text(
      extent={{-150,100},{150,60}}, 
      lineColor={0,0,255}, 
      textString="%name"), 
     Text(
      extent={{-40,22},{40,-18}}, 
      lineColor={0,0,0}, 
      textString="n = %n")}), 
     Diagram(coordinateSystem(preserveAspectRatio=false))); 

end FVM_staggered_Ncells; 

正確的結果表明：

Answer 1

Alrighty ..一些挖後，我想通了。下面我顯示了「收到」代碼，然後在我的編輯下面。希望這能解決所有問題。

背景，如你所知，有一個非常重要的模型結構。你建模的是av_vb。

1.糾正流動模型

的可變分升（流段的長度）的長度爲av_vb模型結構的第一和最後一個體積不同。這種修正對於運行的案例來說是最重要的。

添加了以下修改：

// Define the variable 
SI.Length dLs[n-1]; 
SI.Momentum I[n-1] = mflow[2:n] .* dLs; // Changed from *dL to .*dLs 

// Add to equation section 
dLs[1] = dL + 0.5*dL; 
dLs[2:n-2] = fill(dL,n-3); 
dLs[n-1] = dL + 0.5*dL; 

2.從DPF更改爲MFLOW計算

我跑了一個恆定的流量計算的簡單情況和檢查結果，發現他們是不同的，甚至第一次更正。看起來在指定設置下使用dpf = f（mflow）計算時，「一對一」比較將使用mflow = f（dpf）。這是因爲您選擇了momentumDynamics=SteadyStateInitial，它在PartialGenericPipeFlow中產生了from_dp = true。如果改變它，那麼結果對於恆定流量示例將是相同的（兩者之間的差異將更容易顯示，因爲它們不被流速變化的與時間相關的動力學所掩蓋）。

此外，使用的平均密度與我認爲的MSL管道不同。這並沒有影響這個例子的計算，所以請隨時仔細檢查我的結論。

if use_fixed_zeta then 
    dpf = 1/2*zeta/(n - 1)*(mflow[2:n]) .^ 2 ./ (0.5*(rho[1:n - 1] + rho[2:n])* 
     A*A); 
    else 

// This was the original 
    //  dpf = homotopy(
    //  actual = Modelica.Fluid.Pipes.BaseClasses.WallFriction.Detailed.pressureLoss_m_flow(
    //   m_flow = mflow[2:n], 
    //   rho_a = rho[1:n-1], 
    //   rho_b = rho[2:n], 
    //   mu_a = mu[1:n-1], 
    //   mu_b = mu[2:n], 
    //   length = dLs, //Notice changed dL to dLs 
    //   diameter = D, 
    //   roughness = R, 
    //   m_flow_small = 0.001), 
    //  simplified = dp_nominal/(n-1)/m_flow_nominal*mflow[2:n]); 

// This is the correct model for "one-to-one" comparison for the chosen conditions. Averaged rho and mu was used since useUpstreamScheme = false. 
    mflow[2:n] = homotopy(actual= 
     Modelica.Fluid.Pipes.BaseClasses.WallFriction.Detailed.massFlowRate_dp(
     dpf, 
     0.5*(rho[1:n - 1] + rho[2:n]), 
     0.5*(rho[1:n - 1] + rho[2:n]), 
     0.5*(mu[1:n - 1] + mu[2:n]), 
     0.5*(mu[1:n - 1] + mu[2:n]), 
     dLs, 
     D, 
     A, 
     R, 
     1e-5, 
     4000), simplified=m_flow_nominal/dp_nominal .* dpf); 
    end if; 

3.正確port_b.m_flow參考

這是另一個編輯並沒有影響這個計算的結果，但可能在其他國家。

// Original 
    mflow[n] = -port_b.m_flow; 
// Fixed to reference proper flow variable 
    mflow[n+1] = -port_b.m_flow; 

下面是您生成的圖。地塊重疊。