Fix density in StaticHead models

[fcasella]

Currently StaticHead uses average density. This can be substantially wrong when flowing from a hot to a cold volume, or vice-versa.
A better solution could be using functions such as regStep to weight both densities according to the direction of flow, taking into account a nominal flow direction.
Test in two cases:
a) mass flow through component is non-zero at steady state
b) mass flow through component is zero at steady state (accumulator)

[msielemann]

This ticket was fixed in r1630. Will post more information soon.

[msielemann]

Comparison: Detailed pressure loss correlation

[msielemann]

What has been done?
The pressure loss correlations were rewritten because the existing ones did not address static head and a proper treatment was not possible without including static head exactly at this point. At the same time, several improvements were made to the correlations, most notably the Detailed correlation (which had a knee at flow reversal in the general case).

Why was is necessary?
In the fluid group, we thought about keeping the static head and friction contributions to the pressure drop separate. In the end, this was not possible due to the following reasons.

The sum of two separately regularized functions (e.g. regStep for static head, regRoot2 for friction) is not necessarily properly regularized (e.g. monotonicity).

Consider from_dp = true. Then

dp_grav = rho(sign(m_flow)) * g * h;
m_flow = f(dp_fric) = f(dp - dp_grav);

Essentially, this means that

dp_grav = f(m_flow);
m_flow = f(dp);

Note that it is not possible to use relations such as sign(m_flow)=sign(dp_fric) in a useful way because establishing dp_fric requires dp_grav, which would lead to an implicit equation if substituted above.
As a consequence, an inefficient fixed point iteration is the result of using from_dp = true with static head.

Details: Regularization algorithm
It was not possible to continue using the Fritsch Carlson algorithm for the regularization. The reason is that for regRoot2 and regSquare2 the monotonicity condition on the derivatives of the connection points between the root/square function and the cubic polynomials were fulfilled identically. This is not the case if considering static head. As a consequence, the regularization domain had to be enlarged in certain cases to fulfill the F-C conditions not only on the origin but also on the other two points. While it was possible to find a closed form solution for this lower limit of the regularization domain width for the laminar and quadratic turbulent correlations, no such solution was found for the detailed pressure loss correlation. Using a heuristic was judged as too crude.

Instead, a different algorithm was used (Gasparo Morandi). Here, it is first attempted to use a single cubic polynomial to regularize as demanded. If the cubic is not monotonic, two cubics with a linear section in between are implemented instead. This is very nice because, as a consequence, the exact non-regularized solution can be obtained in the limit of x_small -> 0.
The algorithm had to be slightly extended however to avoid zero/infinite slopes, which lead to integrator step size reduction (e.g. a saddle point does not conflict with the co-monotonicity requirement but has to be avoided). For this purpose, two additional conditions had to be introduced. One is that the slope of the linear section shall always be at least a tenth of the slope of the secant through the points of interest, the other that if the slope is increased using said condition that the start of the right cubic is to the right of the end of the left cubic.

Illustration
Consider the detailed correlation for example. The red line shows the old correlation where density and dynamic viscosity were averaged in several places. Most notably, the static head uses average density and thus presents an approximation only. The blue curve shows the newly implemented correlation.

Also note that the red curve is not C1 smooth. It has a knee at flow reversal, the most critical point for pressure loss correlations. The discontinuity in the first derivative of the old correlation may not seem to be that critical in the graph, but this is only due to the model parameters (IF97 water at atmospheric pressure and 10°C and 90°C temperature respectively).

[msielemann]

Another comment: The new functions end with _staticHead. My suggestion is to replace the old functions with them but I did not delete them yet to allow for comparisons (in WallFractionAndGravitiy just comment/uncomment the respective parts).

Also, I would like to ask everybody for feedback on how well the correlations work. I changed many things and did quite some testing but need the feedback to make everything bulletproof.

[msielemann]

There was a question why the following code is still included in WallFrictionAnGravity (note that d is not in the argument list of the _staticHead functions)

  if allowFlowReversal then
     d = (d_a + d_b)/2 
      "temporary solution that must be improved (BUT NOT WITH if port_a.m_flow>0 then d_a else d_b)";
  else
     d = d_a;
  end if;

The reason is to allow the other members of the fluid group to easily compare the new implementation to the old one. Just comment/uncomment any one of

  if from_dp and not WallFriction.dp_is_zero then
    m_flow = WallFriction.massFlowRate_dp_staticHead(dp, d_a, d_b, eta_a, eta_b, length, diameter,
      g_times_height_ab, roughness, if use_x_small_staticHead then dp_small_staticHead else dp_small);
  else
    dp = WallFriction.pressureLoss_m_flow_staticHead(m_flow, d_a, d_b, eta_a, eta_b, length, diameter,
      g_times_height_ab, roughness, if use_x_small_staticHead then m_flow_small_staticHead else m_flow_small);
  end if;

or

  if from_dp and not WallFriction.dp_is_zero then
     m_flow = WallFriction.massFlowRate_dp(dp-height_ab*ambient.g*d,
                                           d_a, d_b, eta_a, eta_b, length, diameter, roughness, dp_small);
  else
     dp = WallFriction.pressureLoss_m_flow(m_flow, d_a, d_b, eta_a, eta_b, length, diameter, roughness, m_flow_small)
          + height_ab*ambient.g*d;
  end if;

to get the results of the old or the new formulation.

As soon as the testing is complete, the *_staticHead function shall replace the old ones.