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

Question

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

这个模型应该更容易理解，因为它是一个平面模型，而不是从其他类扩展的。这一点很重要，因为因此，即使没有Modelica的知识，我的同事也可以理解模型，我可以让他们相信，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.更正流模型的长度

变量dL (流段的长度)对于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计算

我运行了一个简单的情况，用恒定的流量计算，并检查了结果，发现它们是不同的，即使是第一次修正。当在指定的设置下，“一对一”比较将使用mflow=f( dpf )时，似乎使用dpf= f(mflow)计算。这是因为您选择了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;

下面是您生成的相同的地块。情节重叠。

​

​