Modèle de tuyau dynamique en MSL, volume fini Méthode

Question

J'essaie d'utiliser Modelica pour la modélisation d'un système composé de tuyaux élastiques. Pour l'instant, j'essaye d'implémenter mon propre modèle de pipe dynamique (rigide, pas encore élastique) en utilisant la même approche (volume fini, décalé) comme dans la bibliothèque Modelica.Fluid, mais bien sûr sans inclure toutes les options.

Ce modèle devrait être plus simple à comprendre, car il s'agit d'un modèle plat, ne s'étendant pas d'autres classes. Ceci est important parce que mes collègues peuvent donc comprendre le modèle même sans Modelica Knowhow et je peux les convaincre que Modelica est l'outil adéquat pour nos objectifs!

En tant que test, j'utilise une source de débit massique avec un signal de pas (waterhammer). Mon modèle ne donne pas les mêmes résultats que le composant Modelica.Fluid. J'apprécierais vraiment, si quelqu'un peut m'aider, de comprendre ce qui se passe!

Le système de test ressemble à ceci:

Les résultats pour 11 cellules sont ceci:

Comme vous pouvez le voir, le pic de pression est plus élevé pour le composant MSL et la fréquence/période n'est pas la même. Lorsque je choisis plus de cellules, l'erreur diminue.

Je suis certain que j'utilise exactement les mêmes équations. Pourrait-il être la cause de raisons numériques (j'ai essayé en utilisant des valeurs nominales)? J'ai également inclus mon propre modèle de flux "fixed zeta" pour le composant Modelica.Fluid afin que je puisse le comparer en cas de coefficient de perte de pression fixe zeta.

Le code de mon modèle de conduite est assez courte et ce serait vraiment bien si je reçois cela fonctionne comme ceci:

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; 

Je suis aux prises avec ce problème depuis un temps assez long, donc tous les commentaires ou astuces sont vraiment appréciés !! Si vous avez besoin de plus d'informations ou de résultats de test, s'il vous plaît dites-moi!

Ce code pour l'exemple de test:

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; 

J'ai couru le test avec 301 cellules:

Zoom pic 1 et 2:

Solution: Modifications comme suggéré par 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; 

Résultat droit:

Answer 1

Alrighty .. après quelques recherches, je compris. Ci-dessous, j'ai montré le code "tel que reçu" et ensuite mon édition ci-dessous. Espérons que cela corrige tout. Contexte, comme vous le savez, il existe une structure de modèle très importante. Celui que vous avez modélisé était av_vb.

1. Corriger la longueur du modèle de flux

La variable dL (longueur des segments de flux) est différent pour le premier et le dernier volume d'une structure de modèle av_vb. Cette correction est la plus importante pour le cas qui a été exécuté.

Ajouter la modification suivante:

// 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. Changement de dpf à MFlow calcul

J'ai couru un cas simple avec un calcul de débit constant et vérifié les résultats et a trouvé qu'ils étaient différents, même avec la première correction. Il semble que le calcul dpf = f (mflow) ait été utilisé lorsque, dans les paramètres spécifiés, la comparaison «one-to-one» serait d'utiliser mflow = f (dpf). C'est parce que vous avez choisi momentumDynamics=SteadyStateInitial ce qui fait from_dp = true dans PartialGenericPipeFlow. Si vous le modifiez, les résultats seront les mêmes pour l'exemple de flux constant (les différences entre les deux apparaîtront plus facilement puisqu'ils ne sont pas masqués par la dynamique dépendant du temps des changements de débit).

En outre, les densités moyennes qui étaient utilisées différaient de la conduite MSL je pense. Cela n'a pas d'impact sur le calcul de cet exemple, alors n'hésitez pas à vérifier ma conclusion.

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. Corriger référence port_b.m_flow

Ceci est une autre édition qui n'a pas d'incidence sur les résultats de ce calcul, mais pourrait dans d'autres.

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

Ci-dessous est le même graphique que vous avez généré. Les parcelles se chevauchent.