Wrong enthalpies calculated by DynamicPipe

[meyerg1m]

In ideal gases the enthalpies at in and outlet should be the same. This is valid only if nNodes in an instantiation of DynmamicPipe is set to 1, or a staticPipe is employed. If it is set to 2 the enthalpies differ and a temperature increase at the pipe outlet can be observed.

Increasing the number of nodes (nNodes>20) reduces the deviation but increases the calculation time as a matter of fact.

For real gases the behavior is analogous, but not so easy to detect due to the Joule Thompson effect.

An example is provided below.

model issue
  inner Modelica.Fluid.System system
    annotation (Placement(transformation(extent={{72,-96},{92,-76}})));
  Modelica.Fluid.Sources.MassFlowSource_T boundary(
    redeclare package Medium = Medium,
    m_flow=2.7778E-3,
    T=293.15,
    nPorts=1)
    annotation (Placement(transformation(extent={{-90,-10},{-70,10}})));
  replaceable package Medium = Modelica.Media.Air.DryAirNasa constrainedby
    Modelica.Media.Interfaces.PartialMedium annotation (
      __Dymola_choicesAllMatching=true);
  Modelica.Fluid.Pipes.DynamicPipe pipe(
    redeclare package Medium = Medium,
    allowFlowReversal=true,
    length=30,
    diameter=0.01,
    p_a_start=150000,
    p_b_start=100000,
    T_start=293.15,
    nNodes=2,
    modelStructure=Modelica.Fluid.Types.ModelStructure.a_v_b)
    annotation (Placement(transformation(extent={{-10,-10},{10,10}})));
  Modelica.Fluid.Sources.FixedBoundary boundary1(
    redeclare package Medium = Medium,
    p=100000,
    T=293.15,
    nPorts=1) annotation (Placement(transformation(extent={{84,-10},{64,10}})));
  Modelica.Fluid.Sensors.Temperature temperature(redeclare package Medium =
        Medium) annotation (Placement(transformation(extent={{8,18},{28,38}})));
equation
  connect(pipe.port_b, boundary1.ports[1])
    annotation (Line(points={{10,0},{64,0}}, color={0,127,255}));
  connect(boundary.ports[1], pipe.port_a)
    annotation (Line(points={{-70,0},{-10,0}}, color={0,127,255}));
  connect(temperature.port, pipe.port_b) annotation (Line(points={{18,18},{14,
          18},{14,0},{10,0}}, color={0,127,255}));
  annotation (
    Icon(coordinateSystem(preserveAspectRatio=false)),
    Diagram(coordinateSystem(preserveAspectRatio=false)),
    experiment(StopTime=10, __Dymola_Algorithm="Cvode"));
end issue;

[mestinso]

In ideal gases the enthalpies at in and outlet should be the same.

Whether or not there is an issue here or not in the Modelica Fluid library, I would comment this premise is not correct.

In general, total enthalpy at the inlet and outlet should be the same for this scenario (adiabatic, steady state, steady flow, no shaft work, 1 inlet/1 outlet pipe, z1=z2). So we are left with: h_in + V_in^2/2 = h_out + V_out^2/2

[meyerg1m]

There is probably some confusion here between ideal and perfect gases. Admittedly, this contributes only indirectly to the solution of the problem, therefore only a short statement of the thermodynamic facts is given below.

An ideal gas is defined as such: h=cp(T) *T Since the temperatures are the same, and cp, as the heat capacity, is a function of the temperature only, the enthalpies at in- and outlet have to be the same. The gas ‘DryAirNasa’ behaves mostly (i.e. other classes from Modelica-Standard-Library) as an ideal gas and since it extends ‘IdealGases.Common.SingleGasNasa’, I assume this gas is implemented as an ideal gas. This assumption is supported by the code snippets below.

redeclare function extends specificHeatCapacityCp
  "Return specific heat capacity at constant pressure"
algorithm 
  cp := Modelica.Media.IdealGases.Common.Functions.cp_T(
             data, state.T);
  annotation(Inline=true,smoothOrder=2);
end specificHeatCapacityCp;

function cp_T
  "Compute specific heat capacity at constant pressure from temperature and gas data"
  extends Modelica.Icons.Function;
  input IdealGases.Common.DataRecord data "Ideal gas data";
  input SI.Temperature T "Temperature";
  output SI.SpecificHeatCapacity cp "Specific heat capacity at temperature T";
algorithm 
  cp := smooth(0,if T < data.Tlimit then data.R*(1/(T*T)*(data.alow[1] + T*(
    data.alow[2] + T*(1.*data.alow[3] + T*(data.alow[4] + T*(data.alow[5] + T
    *(data.alow[6] + data.alow[7]*T))))))) else data.R*(1/(T*T)*(data.ahigh[1]
     + T*(data.ahigh[2] + T*(1.*data.ahigh[3] + T*(data.ahigh[4] + T*(data.
    ahigh[5] + T*(data.ahigh[6] + data.ahigh[7]*T))))))));
  annotation (Inline=true,smoothOrder=2);
end cp_T;

For perfect gases the heat capacity cp is a function of the pressure and for real gases cp is a function of pressure and temperature.

[AlesVojacek]

@mestinso You pointed that out that we are left with h_in + V_in^2/2 = h_out + V_out^2/2, but I believe that Modelica.Fluid.Pipes.DynamicPipe neglects the kinetic term in the energy balance so I believe that the h_in should be equal to h_out in this case.

[mestinso]

@mestinso You pointed that out that we are left with h_in + V_in^2/2 = h_out + V_out^2/2, but I believe that Modelica.Fluid.Pipes.DynamicPipe neglects the kinetic term in the energy balance so I believe that the h_in should be equal to h_out in this case.

Based on the comments here: Modelica.Fluid.UsersGuide.ComponentDefinition.BalanceEquations, I think the dynamic pipe is the one spot in the library where the kinetic energy effect is taken into account. As a test of this, if you swap the dynamic pipe with the static pipe (where the documentation indicates the kinetic effect is ignored), sure enough h_in = h_out.

[AlesVojacek]

Thanks for pointing that out. This "alternative" form of energy balance equation is quite new to me. What still puzzles me is that when you take the above mention model "issue" and increase number of elements, the enthalpy difference fades out, e.g. for 30 elements, (i.e. 1 element/1 meter length), it is almost zero. That should not happen if it has correctly implemented kinetic terms, or?

[rfranke]

I resigned from Fluid maintenance mid of last year. The paper: Franke, Casella, Sielemann, Proelss, Otter, Wetter: Standardization of Thermo-Fluid Modeling in Modelica.Fluid, Proceedings 7th Modelica Conference, Como, Italy, Sep. 20-22, 2009 describes the assumptions of the DynamicPipe model.

[mestinso]

One lead I'm seeing regarding the kinetic energy term:

parameter Boolean use_Ib_flows = momentumDynamics <> Types.Dynamics.SteadyState
          "= true to consider differences in flow of momentum through boundaries"
               annotation(Dialog(group="Advanced"), Evaluate=true);

... so apparently, the kinetic term is only enabled when momentumDynamics is not set to steady-state (by default). Not sure why that is the case, and arguably that should be made more clear. Regardless, h_in should equal h_out in the problem setup by @meyerg1m since momentumDynamics is set to steady-state. Also, I did run a few tests and was able to confirm that when momentum dynamics == false and significant discretization, h_in = h_out. And when momentumDynamics = true, h_in + V_in^2/2 = h_out + V_out^2/2.

Unfortunately, as OP described, it does seem surprising how much discretization is required to get an accurate result. Also, the model structure choices seem to have a very significant effect here. For example, av_vb (instead of a_v_b used in the example) seems to behave much better numerically. It's unclear to me if there are bugs here or not or if the dynamic pipe component just requires very careful setup.

[aci31]

I was quite intrigued by this, and by digging a bit, I discovered that it comes from (at least) the calculation of mechanical power Wb_flows (the last array value is always high whereas the other ones are close to 0). If you look at the definition of Wb_flows in the model Modelica.Fluid.Pipes.DynamicPipe, it is of course dependant of the schema selected for the pipe, but the formulation is quite strange at the boundary for structures av_b, a_vb and a_v_b

if n == 1 or useLumpedPressure then
    Wb_flows = dxs * ((vs*dxs)*(crossAreas*dxs)*((port_b.p - port_a.p) + sum(flowModel.dps_fg) - system.g*(dheights*mediums.d)))*nParallel;
  else
    if modelStructure == ModelStructure.av_vb or modelStructure == ModelStructure.av_b then
      Wb_flows[2:n-1] = {vs[i]*crossAreas[i]*((mediums[i+1].p - mediums[i-1].p)/2 + (flowModel.dps_fg[i-1]+flowModel.dps_fg[i])/2 - system.g*dheights[i]*mediums[i].d) for i in 2:n-1}*nParallel;
    else
      Wb_flows[2:n-1] = {vs[i]*crossAreas[i]*((mediums[i+1].p - mediums[i-1].p)/2 + (flowModel.dps_fg[i]+flowModel.dps_fg[i+1])/2 - system.g*dheights[i]*mediums[i].d) for i in 2:n-1}*nParallel;
    end if;
    if modelStructure == ModelStructure.av_vb then
      Wb_flows[1] = vs[1]*crossAreas[1]*((mediums[2].p - mediums[1].p)/2 + flowModel.dps_fg[1]/2 - system.g*dheights[1]*mediums[1].d)*nParallel;
      Wb_flows[n] = vs[n]*crossAreas[n]*((mediums[n].p - mediums[n-1].p)/2 + flowModel.dps_fg[n-1]/2 - system.g*dheights[n]*mediums[n].d)*nParallel;
    elseif modelStructure == ModelStructure.av_b then
      Wb_flows[1] = vs[1]*crossAreas[1]*((mediums[2].p - mediums[1].p)/2 + flowModel.dps_fg[1]/2 - system.g*dheights[1]*mediums[1].d)*nParallel;
      Wb_flows[n] = vs[n]*crossAreas[n]*((port_b.p - mediums[n-1].p)/1.5 + flowModel.dps_fg[n-1]/2+flowModel.dps_fg[n] - system.g*dheights[n]*mediums[n].d)*nParallel;
    elseif modelStructure == ModelStructure.a_vb then
      Wb_flows[1] = vs[1]*crossAreas[1]*((mediums[2].p - port_a.p)/1.5 + flowModel.dps_fg[1]+flowModel.dps_fg[2]/2 - system.g*dheights[1]*mediums[1].d)*nParallel;
      Wb_flows[n] = vs[n]*crossAreas[n]*((mediums[n].p - mediums[n-1].p)/2 + flowModel.dps_fg[n]/2 - system.g*dheights[n]*mediums[n].d)*nParallel;
    elseif modelStructure == ModelStructure.a_v_b then
      Wb_flows[1] = vs[1]*crossAreas[1]*((mediums[2].p - port_a.p)/1.5 + flowModel.dps_fg[1]+flowModel.dps_fg[2]/2 - system.g*dheights[1]*mediums[1].d)*nParallel;
      Wb_flows[n] = vs[n]*crossAreas[n]*((port_b.p - mediums[n-1].p)/1.5 + flowModel.dps_fg[n]/2+flowModel.dps_fg[n+1] - system.g*dheights[n]*mediums[n].d)*nParallel;
    else
      assert(false, "Unknown model structure");
    end if;
  end if;

It seems to me they are not using the right pressure reference. I corrected it somehow (as shown below) and it kind of solve the enthalpy conservation problem... (I am not sure my correction is covering every possible use case if some heights are used for instance):

if n == 1 or useLumpedPressure then
    Wb_flows = dxs * ((vs*dxs)*(crossAreas*dxs)*((port_b.p - port_a.p) + sum(flowModel.dps_fg) - system.g*(dheights*mediums.d)))*nParallel;
  else
    if modelStructure == ModelStructure.av_vb or modelStructure == ModelStructure.av_b then
      Wb_flows[2:n-1] = {vs[i]*crossAreas[i]*((mediums[i+1].p - mediums[i-1].p)/2 + (flowModel.dps_fg[i-1]+flowModel.dps_fg[i])/2 - system.g*dheights[i]*mediums[i].d) for i in 2:n-1}*nParallel;
    else
      Wb_flows[2:n-1] = {vs[i]*crossAreas[i]*((mediums[i+1].p - mediums[i-1].p)/2 + (flowModel.dps_fg[i]+flowModel.dps_fg[i+1])/2 - system.g*dheights[i]*mediums[i].d) for i in 2:n-1}*nParallel;
    end if;
    if modelStructure == ModelStructure.av_vb then
      Wb_flows[1] = vs[1]*crossAreas[1]*((mediums[2].p - mediums[1].p)/2 + flowModel.dps_fg[1]/2 - system.g*dheights[1]*mediums[1].d)*nParallel;
      Wb_flows[n] = vs[n]*crossAreas[n]*((mediums[n].p - mediums[n-1].p)/2 + flowModel.dps_fg[n-1]/2 - system.g*dheights[n]*mediums[n].d)*nParallel;
    elseif modelStructure == ModelStructure.av_b then
      Wb_flows[1] = vs[1]*crossAreas[1]*((mediums[2].p - mediums[1].p)/2 + flowModel.dps_fg[1]/2 - system.g*dheights[1]*mediums[1].d)*nParallel;
      //Wb_flows[n] = vs[n]*crossAreas[n]*((port_b.p - mediums[n-1].p)/1.5 + flowModel.dps_fg[n-1]/2+flowModel.dps_fg[n] - system.g*dheights[n]*mediums[n].d)*nParallel;
      Wb_flows[n] = vs[n]*crossAreas[n]*(
        (mediums[n].p - mediums[n-1].p)/2 + flowModel.dps_fg[n-1]/2+(port_b.p - mediums[n].p) + flowModel.dps_fg[n] - system.g*dheights[n]*mediums[n].d)*nParallel;
    elseif modelStructure == ModelStructure.a_vb then
      //Wb_flows[1] = vs[1]*crossAreas[1]*((mediums[2].p - port_a.p)/1.5 + flowModel.dps_fg[1]+flowModel.dps_fg[2]/2 - system.g*dheights[1]*mediums[1].d)*nParallel;
      Wb_flows[1] = vs[1]*crossAreas[1]*(
        (mediums[1].p - port_a.p) + flowModel.dps_fg[1]+(mediums[2].p - mediums[1].p)/2 + flowModel.dps_fg[2]/2 - system.g*dheights[1]*mediums[1].d)*nParallel;
      Wb_flows[n] = vs[n]*crossAreas[n]*((mediums[n].p - mediums[n-1].p)/2 + flowModel.dps_fg[n]/2 - system.g*dheights[n]*mediums[n].d)*nParallel;
    elseif modelStructure == ModelStructure.a_v_b then
      //Wb_flows[1] = vs[1]*crossAreas[1]*((mediums[2].p - port_a.p)/1.5 + flowModel.dps_fg[1]+flowModel.dps_fg[2]/2 - system.g*dheights[1]*mediums[1].d)*nParallel;
      Wb_flows[1] = vs[1]*crossAreas[1]*(
        (mediums[1].p - port_a.p) + flowModel.dps_fg[1]+(mediums[2].p - mediums[1].p)/2 + flowModel.dps_fg[2]/2 - system.g*dheights[1]*mediums[1].d)*nParallel;
      //Wb_flows[n] = vs[n]*crossAreas[n]*((port_b.p - mediums[n-1].p)/1.5 + flowModel.dps_fg[n]/2+flowModel.dps_fg[n+1] - system.g*dheights[n]*mediums[n].d)*nParallel;
      Wb_flows[n] = vs[n]*crossAreas[n]*(
        (mediums[n].p - mediums[n-1].p)/2 + flowModel.dps_fg[n]/2+(port_b.p - mediums[n].p) + flowModel.dps_fg[n+1] - system.g*dheights[n]*mediums[n].d)*nParallel;
    else
      assert(false, "Unknown model structure");
    end if;
  end if;

That could be a lead...

[hubertus65]

Arnaud, you might be right there. This error is, in my mind, a symptom of what I would call a design flaw of the dynamic pipe model. The problem is that there are, if I remember my counting correctly, 2^14 = 16384 different configurations of the model possible, depending on model structure and assumptions. Such a model is untestable for all combinations of model choices and assumptions, so there are bound to be errors in it. You probably found one of them, and there are probably more. Thanks for looking into this. I think our challenge is here to set up tests that cover a sufficient subset of the model option combinations (all are clearly unfeasible), and show that there are no errors at least in those. From visual inspection, your fix looks good, but unfortunately, I don't have the time to properly test and fix this right now. Many thanks for digging!

[mestinso]

@aci31 Yes, I think you nailed it. The idea of the 1.5 in the term you corrected is right in the sense that the length is 1.5, but the issue is the derivative is centered about the wrong location with that method. Using the idea of ghost nodes (ie a node outside the domain), I derived the same thing you did.

For example, for the a_vb Wb_flows[1] term, we want: v[1]*(medium.p[2] - medium.p[0])/2

medium.p[0] is the ghost node here. medium.p[0] can be estimated as port_a.p - (medium.p[1] - port_a.p) = 2*port_a.p - medium.p[1]

...plugging in above we get: v[1]*((medium.p[1] + medium.p[2])/2 - port_a.p), which effectively works out the same as what you have and you can tell by inspection that is centered around location 1.

[mestinso]

I investigated pretty thoroughly and issued a PR with the fix discussed. It seems to fully address the enthalpy balance when the kinetic energy term is turned off.

However, after investigating more, there are more discretization issues that affect a) the case when kinetic energy is turned on, and b) the pressure drop predictions in general. The most obvious hint of remaining issues is that the dp/dx[1] term has a 2 (clearly too large) in the denominator for av_vb and av_b cases. The first node [1] is off-center located at the port, and the next node [2] is centered, so the denominator should be 1.5. The same issue is present for the dps_fg term (should have a denominator of 1.5). Also, these issues are present for the last dp/dx and dps_fg terms for av_vb and a_vb. The issues don't stop there either, since for dp/dx[2], for the terms between between nodes [1] and [2] should have a 3 in the denominator instead of 2, etc. Sort of a pain to explain, but I have a fix lined up to get these remaining issues resolved. Will plan on adding to the existing PR.

[mestinso]

Update: the PR has been updated to resolve the remaining issues.

Additionally, I would comment that some of the previous commentary and proposed solutions for the a_vb and av_b cases were incorrect since they don't have "half momentum balances" like the a_v_b case does.

To summarize the cases:

av_vb: end nodes are on the far edge of the end volumes internal nodes are at the center of the internal volumes ...therefore there is some special handling via if statements to handle the different flow lengths that result

a_vb and av_b: the nodes are on the edges of the volumes ...therefore all the flow lengths are the same and these ones end up pretty simple looking

a_v_b: nodes are all at the center of the volumes Plus there are "half momentum balances" between the end nodes and the ports

Given the above facts plus the fact that the volume lengths are equal to LengthTotal/nVolumes, the fixed Wb_flows in the PR should be a natural conclusion. Ultimately, the fixes lead to fixed energy balance (resolving the original issue here), and more accurate results in general.