Wrong mass flow rate calculation in Modelica.Fluid.Fittings.GenericResistances.VolumeFlowRate

[joergu]

The mass flow rate of volumeFlowRate is always zero if the pressure drop is a linear function of volume flow rate.

The volume flow of is calculated by solving the quadric equation dp=aV_flow^2+bV_flow. This fails if a is zero. The actual problem is in Modelica.Fluid.Dissipation.PressureLoss.General.dp_volumeFlowRate_MFLOW:

1. a is silently replaced by max(a,Modelica.Constants.eps) (the same for b, but this is not an issue here)
2. The volume flow rate is calculated by (-b/(2a) + ((b/(2a))^2 + dp/a)^0.5 (I simplified a bit)
   The point is, for small a the right term sqrt(X^2 + eps) --> X, and this leads to V_flow=0 if b>0.

Proposed solution:

1. The replacement a,b=max(Modelica.Constants.eps, abs(IN_con.a,b)) should be removed. Instead, for (a*a)<=0 the linear solution dp/b should be returned. This should be relevant only if a is given as a parameter.
2. One should calculate only the solution with the greater absolute value with the formula above. In the case b>0 (which is the common case) this is the negative solution
   x1= (-b/(2a) - ((b/(2a))^2 + dp/a)^0.5
   The second root should be calculated by Vieta's formula
   x2 = c/a/x1
   In the common case b>0 the second solution x2 is the physically meaningful.
   I removed the regularized power function SmoothPower here only for simplicity of presentation. It should be kept in the actual code.

[casella]

@wischhusen can you take care of this?

Thanks!

[wischhusen]

Thanks for reporting this issue! I resolved the problem by allowing positive (incl. zero) numbers for a and b. One factor has to be greater zero to allow a proper inversion of the function.

[joergu]

Hi,
in Line 4606 I think it should be for the linear case M_FLOW = rho * dp/b
in Line 4602ff we have still the issue of -b/(2*a) has opposite sign as the square root which might lead to cancellation. Using Vietas formula would avoid that:

V_FLOW = c/a * 1/(b/(2*a) + sqrt(...))

[wischhusen]

I would propose the following solution. We do not need regularization since the argument of sqrt contains the constant term (b/(2*a))^2:

   if b>0 then
     M_FLOW := IN_var.rho*(if a>0 then (-b/(2*a) +
            sqrt((b/(2*a))^2 + (1/a)*dp)) else b*dp);
   else
     M_FLOW := IN_var.rho*sqrt(1/a)*
            Modelica.Fluid.Dissipation.Utilities.Functions.General.SmoothPower(
            dp,
            dp_min,
            0.5);
   end if;

[beutlich]

I would propose the following solution.

Hm, we already have a PR #2782 addressing this issue. 😕

[wischhusen]

Hm, but there is a flaw in that pull request. Shall I put the code there as well?

[beutlich]

Hm, but there is a flaw in that pull request.

It really is hard to tell or reproduce for outsiders without a SSCCE (Short, Self Contained, Correct (Compilable), Example).

Shall I put the code there as well?

You still can update this PR by pushing to the issue2770 branch.