some models fail to translate for Glycol47

[mwetter]

This branch is to integrate https://github.com/lbl-srg/modelica-buildings/issues/1163

Todo:

• [x] add validation case
• [x] merge change in ConservationEquation.mo

[Mathadon]

@mwetter Some analysis on this:

    if simplify_mWat_flow then
      // If moisture is neglected in mass balance, assume for computation
      // of the mass of air that the air is at Medium.X_default.
      m = fluidVolume*Medium.density(Medium.setState_phX(
        p = medium.p,
        h = hOut,
        X = Medium.X_default));
    else
      // Use actual density
      m = fluidVolume*medium.d;
    end if;

The top code branch fails, the bottom one works.

The reason is that the bottom one implicitly defines the temperature (and from this the density) using

h = if enthalpyOfT then h_T(T) else  h_pT(p,T,densityOfT);

function h_pT is differentiable. Dymola then generates an algebraic loop to solve the set of implicit equations, which works.

The top code branch uses Medium.setState_phX containing

state :=ThermodynamicState(p=p,T=T_ph(p,h));

which calls T_ph which computes

T := Internal.solve(h, T_min, T_max, p, {1}, Internal.f_nonlinear_Data());

the solve() function is not differentiable. Therefore the solver fails using

Cannot find differentiation function:
Modelica.Media.Incompressible.Examples.Glycol47.T_ph_Unique5(vol.ports[1].p, vol1.ports[1].h_outflow)
with respect to time

so that makes sense.

So to solve this properly we should be able to define the function T_ph implicitly as h=h_pT(p,T) (such that an algebraic loop is generated instead of using the derivative of solve()), which is not allowed by the Modelica specification as far as I know. I tested it in Dymola and that returns an error Assignment to input variable: h. An alternative is to define the derivative function for solve(), which is something that has to be solved on the MSL side.

Based on this analysis, the pull request is a bit of a workaround of a problem in MSL or in the Modelica specification. This workaround can cause unexpected results in the future and clutters the code. So I want to ask if this is the approach that we want to take here? If yes, fine, but I would also make a Modelica trac issue for both the Modelica specification issue and the MSL issue?

edit: Setting massDynamics=SteadyState also solves the problem so that's another workaround that does not require changing the models. And I suppose Massimo's media also won't experience these problems.

[mwetter]

@Mathadon

The suggestion of defining a derivative function for solve() won't work: This function uses iterations, hence it only approximates its result. The same will be for the derivative, e.g., it is an approximate derivative. Now, some Newton solver implementations require to control the approximation error of the derivative as a function of the rate of convergence. But solve() does not allow this, hence this suggestion is not what one wants to implement.

I agree about the cluttering, and I am not thrilled about it, but don't have a better suggestion. I don't see how this could cause "unexpected results": My change essentially says for media with one component, always use the actual X rather than X_default. But in this case, we always have X = {1}, hence the formulations are mathematically equivalent.

I am not sure what you mean by assigning the input variable returns an error. There must be something else going on here. This model works:

model Unnamed1
  Real x;
equation
  1 = IBPSA.Utilities.Math.Functions.average(nin=2, u = {x, 1.5});
end Unnamed1;

[Mathadon]

In a model it works, but not in a function :)