modelica - MSL 中的动态管道模型，有限体积法

Question

我正在尝试使用 Modelica 对由弹性管道组成的系统进行建模。
现在，我正在尝试使用与 Modelica.Fluid 库中相同的方法(有限体积，交错)来实现我自己的动态管道模型(刚性，尚未弹性)，但当然不包括所有选项。

这个模型应该更容易理解，因为它是一个平面模型，不是从其他类扩展而来的。这很重要，因为即使没有 Modelica Knowhow，我的同事也可以理解模型，我可以说服他们，Modelica 是满足我们目的的合适工具!

作为一个测试案例，我使用了一个带有阶跃信号(水锤)的质量流量源。
我的模型给出的结果与 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;

很长一段时间以来，我一直在努力解决这个问题，因此非常感谢任何评论或提示!!
如果您需要更多信息或测试结果，请告诉我!

这是 测试示例代码 :

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 个电池进行测试 :


变焦峰值 1 和 2:




解决方案:按照 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

好吧……经过一番挖掘，我想通了。下面我显示了“收到的”代码，然后在它下面进行了编辑。希望这能解决所有问题。

背景，如您所知，有一个非常非常重要的模型结构。您建模的那个是 av_vb .

1、修正流模型长度
av_vb 的第一个和最后一个体积的变量 dL(流段的长度)不同模型结构。对于运行的案例，此更正是最重要的。

添加以下修改:

// 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计算

我运行了一个带有恒流计算的简单案例，并检查了结果，发现即使进行了第一次修正，结果也有所不同。当在指定设置下“一对一”比较将使用 mflow=f(dpf) 时，似乎使用了 dpf = f(mflow) 计算。这是因为您选择了 momentumDynamics=SteadyStateInitial这使得 from_dp = true在 PartialGenericPipeFlow .如果你改变它，那么恒流示例的结果将是相同的(两者之间的差异将更容易显示出来，因为它们不会被变化的流速的时间相关动态所掩盖)。

此外，使用的平均密度与我认为的 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;

下面是您生成的同一个图。情节重叠。