R1233zd does not work in the saturation region close to the bubble point

[fedetftpolimi]

Try the following minimal working example

#include "CoolProp.h"
#include <stdio.h>

using namespace CoolProp;

int main()
{
    printf("%f\n",PropsSI("T", "P", 136730, "H", 260000.00, "R1233zd(E)"));
    printf("%f\n",PropsSI("T", "P", 136730, "H", 263277.27, "R1233zd(E)"));
    printf("%f\n",PropsSI("T", "P", 136730, "H", 263277.35, "R1233zd(E)"));
    printf("%f\n",PropsSI("T", "P", 136730, "H", 268000.00, "R1233zd(E)"));
    printf("%f\n",PropsSI("T", "P", 136730, "H", 275000.00, "R1233zd(E)"));
    printf("%f\n",PropsSI("T", "P", 136730, "H", 293000.00, "R1233zd(E)"));
    printf("%f\n",PropsSI("T", "P", 136730, "H", 300000.00, "R1233zd(E)"));
    return 0;
}

On my machine ti produces the following results

$ g++ -o t -DNO_FMTLIB -Iinclude test.cpp -L. -lCoolProp -ldl
$ ./t
296.914485
299.609517
inf
inf
inf
299.609875
299.609875

The expected behavior is to see 299.609875 in place of inf

@jowr could you please have a look, I'm keeping @casella and @albertoleva in the loop.

[fedetftpolimi]

If it could be of any help, the same issue occurs when using CoolProp through ExternalMedia,
attached you can find the Modelica code

model Mnwe
extends Modelica.Icons.Example;

package Medium
extends ExternalMedia.Media.CoolPropMedium(
    mediumName = "R1233",
    substanceNames = {"R1233zd(E)|calc_transport=0"},
    ThermoStates = Modelica.Media.Interfaces.Choices.IndependentVariables.ph,
    SpecificEnthalpy(start=2e5));
end Medium;

Medium.AbsolutePressure p1;
Medium.SpecificEnthalpy h1;
Medium.ThermodynamicState state1;
ExternalMedia.Test.GenericModels.CompleteThermodynamicState completeState1(
    redeclare package Medium = Medium, state=state1);
parameter Medium.SpecificEnthalpy hstart=260000;
parameter Medium.SpecificEnthalpy hend  =268000;
equation
state1 = Medium.setState_ph(p1, h1);
p1 = 136730;
h1 = hstart+(hend-hstart)*time;
end Mnwe;

It fails with the following error that we believe is the same underlying issue that causes inf to be printed when calling CoolProp directly:

ExternalMedia error: unable to solve 1phase PY flash with Tmin=195.149, Tmax=299.61 due to error: HSU_P_flash_singlephase_Brent could not find a solution because Hmolar [34356.6 J/mol] is above the maximum value of 34356.2103383 J/mol

and this error is cused by an exception being thrown from
/externals/CoolProp.git/src/Backends/Helmholtz/FlashRoutines.cpp:1548

[fedetftpolimi]

Hi,
tested on the newly released version 6.4.3, the problem is still present:

299.609517
inf
inf
inf
299.609875
299.609875

any updates on this issue?

[ibell]

@msaitta-mpr do you think you could dig into this a little bit?

[fedetftpolimi]

Thanks for looking into this!
Just to clarify, our main objective is to have the R1233zd fluid working in Modelica models through ExternalMedia, however the underlying issue was traced back to CoolProp.
Currently we are still stuck with the following error:
unable to solve 1phase PY flash with Tmin=195.149, Tmax=299.61 due to error: HSU_P_flash_singlephase_Brent could not find a solution because Hmolar [34356.6 J/mol] is above the maximum value of 34356.2103383 J/mol

[msaitta-mpr]

Looks another issue of the ancillary solver giving a bad limit when near the bubble point. I will need to spend some time to look into it further, but a few ideas off the top of my head:

• Adjust the ancillary solver guess by 0.05K to try to get a better bounding value.
• Try to force it to call the saturation solver to get a bound if the initial ancillary fails.
• Give it some sort of improvement in the initial singlephase solver to see if it can get closer.

[msaitta-mpr]

I looked into this a little bit more and there seem to be two issues at play:

1. When very close to the saturation line, the ancillary solver can give a guess that is off by a few thousandths of a Kelvin. Simply increasing the bounds on the Brent solver by 0.01 K fixes this issue.
2. The ancillary equation provides a bad solution for the phase. At the guessed temperature (299.609523 K), the ancillary equation provides a saturated liquid enthalpy of 37624.6 J/mol, but the actual saturated liquid enthalpy is 34356.2 J/mol. This is difference of more than 3000 J/mol. This is greater than the provided maximum error of 271 J/mol.

With regards to item 2, I am unsure what the provenance of the ancillary equation is here. @ibell do you know? I do not have access to the cited EOS paper: https://link.springer.com/article/10.1007/s10765-016-2040-6.

[msaitta-mpr]

See the following image to illustrate the issue with the current ancillary curve. The blue line indicates where the current ancillary function estimates the saturation line, and the orange line indicates where a TQ update indicates the saturation curve to be. The green curve illustrates the ancillary curve less the provided error limit. Any point that falls into the space between the curves will erroneously use the subcooled liquid solver, whereas it should use the two-phase solver.

[casella]

Is this behaviour specific of the R1233zd refrigerant, or is it a commonplace problem?

[msaitta-mpr]

I think it is just R1233zd. Looking at another arbitrary fluid (ammonia), the ancillary works well.

[casella]

Does that mean that the ancillary curve coefficients are not that good and should be improved?

[ibell]

You might want to just replace the ancillaries with the ones that are in the published EOS: https://aip.scitation.org/doi/suppl/10.1063/5.0083026 . That would be at least a worthwhile test.

[msaitta-mpr]

It looks like the linked ancillaries do not contain the function of interest. The function of current interest is the molar enthalpy of a saturated liquid as a function of temperature (hL). Do you know how these curves were originally generated? If not, I will go ahead and regress a polynomial of my own to see how well it performs and if it solves this issue.

[ibell]

First try deleting the ancillary for enthalpy. Then we can talk about refitting.

[msaitta-mpr]

Deleting the enthalpy ancillary fixes the issue and the expected values are returned.

[casella]

@ibell how can we get this improvement into ExternalMedia, eventually?

[ibell]

@msaitta-mpr knows how to make a pull request (PR) with the change, I will accept it, and then you can build ExternalMedia with the new CoolProp code

[casella]

Do I need some special git commands to pull the updates to CoolProp into ExternalMedia? Sorry for the dumb question 😅

[msaitta-mpr]

Once the PR is merged into master, you should be able to run the following commands to update your copy of the repository.

git checkout master
git pull

The above assumes that you have set the main copy of the repository as your default and you are not working off of a fork. Once you have pulled the changes, you should be able to build the library. See CoolProp's fantastic documentation on how to build. The library can then be compiled into your program.

I am not familiar with ExternalMedia, so I am not sure if you have a more robust build system that will automatically update and build CoolProp. If so, you should use the appropriate tools that are available.

[casella]

Thanks a lot for the fix, very much appreciated! We'll follow up on ExternalMedia ASAP.

[fedetftpolimi]

Thanks!
we already made a local build on ExternalMedia by replacing the externals/CoolProp.git directory with the lastest CoolProp and it works.

I think that for ExtenalMedia to import the change first a new CoolProp release will have to be made through SourceForge, and then this line needs updating:
https://github.com/modelica-3rdparty/ExternalMedia/blob/master/Projects/CMakeLists.txt#L53

[casella]

I'm not an expert in CMAKE, but it's not clear to me from https://github.com/modelica-3rdparty/ExternalMedia/blob/master/Projects/CMakeLists.txt if the idea was to use the sources from zipped released versions of CoolProp, downloaded from sourceforge (line #53), or rather to use git external to get the source code from git (line #52). And, in that case, which version of the GIT repository is used.

I think this is decided in some other script of the GitHub CI, but I'm only guessing, I'm not familiar with that either.

@jowr you originally wrote the CMAKE script, could you please comment on that? How would you recommend to proceed? Should we stick to released versions of CoolProp, or can we dare using development versions?

[fedetftpolimi]

I did not write the CMake script, I only reverse engineered it as much as I needed to get a working build locally.

However, the CoolProp.git at line 52 seems to have nothing to do with downloading it from git. It's just the name of the zipped directory downloaded from sourceforge once extracted and renamed, see line 64 of the same file.
The logic in cmake seems to be "if there is no directory externals/CoolProp.git then download it from sourceforge, extract the zip and rename the directory CoolProp.git"

Also note that I tried replacing the sourceforge download path with the github path, as also github has the feature to download zip files, but it cannot work as:

• CoolProp requires its git submodules
• the zip downloaded from github does not have them, and also does not include the .git directory
  so basically the zip of the source files downloaded through github is useless as it cannot be compiled because it does not have the submodule sources nor a way to download them.

The zip downloaded through sourceforge also does not have a .git directory, but at least it contains the submodules, I have no idea what magic happens in sourceforge.