alternating results using python wrapper in matlab

[jastema]

Hello,

I'm running into some trouble using the python wrapper in Matlab. I've been using the refpropm function in Matalb up to this point, but since there were some issues with this, I'm now switching to the python wrapper.

I'm trying to calculate different state points of an 50/50 (mass basis) ammonia/water mixture.

When using the REFPROP GUI it yields following state points:

1.) p = 1500 kPa; h = 170 kJ/kg --> T = 328.94 K

2.) p = 1600 kPa; h = 250 kJ/kg --> T = 344.51 K

which I'm now trying to reproduce using the REFPROPdll call in the Matlab environment:

RP = py.ctREFPROP.ctREFPROP.REFPROPFunctionLibrary('C:\Program Files (x86)\REFPROP\REFPRP64.DLL');
iUnit = RP.GETENUMdll(int8(0), 'MKS').iEnum;
iMass = int8(1);
iFlag = int8(1);

r = RP.REFPROPdll('AMMONIA;WATER','PH','T',iUnit,iMass,iFlag,1500,170,{0.5 0.5}); 
o = double(r.Output); 
T_1 = o(1)

q = RP.REFPROPdll('AMMONIA;WATER','PT','H',iUnit,iMass,iFlag,1500,T_1,{0.5 0.5}); 
p = double(q.Output); 
h_1 = p(1)

w = RP.REFPROPdll('AMMONIA;WATER','PH','T',iUnit,iMass,iFlag,1600,250,{0.5 0.5}); 
g = double(w.Output); 
T_2 = g(1)

%% deactivate these lines to alter results   %%

s = RP.REFPROPdll('Water','PQ','T',iUnit,iMass,iFlag,101.325,0.0,{1.0});
k = double(s.Output);
T_3 = k(1)

%% first state point calculated with refpropm does not converge %%
% T = refpropm('T','H',1.7e5,'P',1.5e3,'water','ammonia',[1-0.5 0.5])

%% second one yields same results  %%
% T = refpropm('T','H',2.5e5,'P',1.6e3,'water','ammonia',[1-0.5 0.5])

When I first run this code the results are as follow:

T_1 = 328.94 K h_1 = 170.00 kJ/kg T_2 = 344.50 K T_3 = 373.12 K

where T_1 and T_2 are the above mentioned state points 1.) and 2.) so this is matching with the results of the GUI.

But when I run the code a second time, the results change to:

T_1 = 330.05 K h_1 = 175.71 kJ/kg T_2 = 344.50 K T_3 = 373.12 K

so T_1 and h_1 are changed.

Strange thing is, when i comment/deactivate the lines

s = RP.REFPROPdll('Water','PQ','T',iUnit,iMass,iFlag,101.325,0.0,{1.0});
k = double(s.Output);
T_3 = k(1)

and run the code twice, T_1 and h_1 are matching the GUI results again. When the lines are activated again, the results will change again with the 2nd run of the code.

What could be the reason for this ? Is it possible that is has something to do with the change of fluid ?

As I've said, I'm quite new to the python wrapper. I've also run into troubles calculating the first state point with the refpropm function (does not converge) which was the reason for switching to the python wrapper in the first place. The second state point can be calculated correctly with the refpropm function as well. Maybe that information is of help.

Versions

REFPROP Version: 10.0.0.02 Operating System and Version: Win10 & Win7 Matlab Version: R2016a & R2018a
Python Version: 3.7

[ianhbell]

Have you read this part of the FAQ: https://pages.nist.gov/REFPROP-docs/#reference-states-enthalpy-and-entropy-differences ? The GUI and the DLL handle this differently. So you should never input enthalpy, you should always start with a thermodynamic state point (in T,p, for instance), use that to calculate enthalpy, and then try the backwards calculation with enthalpy and pressure. If that doesn't come back to the same T, that's a bug.

[jastema]

Hello Ian, thanks for your reply. Yes, I was aware of that. The above mentioned code is just an example, which only represents a small part of my whole code. The used enthalpy is calculated beforehand with given temperature, pressure and a change in energy due to heat transfer. But this point was were the problems occured, so I took this as an example. Anyways, if i rearrange the code of my example in order to first calculate enthalpy and the temperature afterwards, the same behaviour occurs.

RP = py.ctREFPROP.ctREFPROP.REFPROPFunctionLibrary('C:\Program Files (x86)\REFPROP\REFPRP64.DLL');
iUnit = RP.GETENUMdll(int8(0), 'MKS').iEnum;
iMass = int8(1);
iFlag = int8(1);

T_0 = 328.9438;

q = RP.REFPROPdll('AMMONIA;WATER','PT','H',iUnit,iMass,iFlag,1500,T_0,{0.5 0.5}); 
p = double(q.Output); 
h_1 = p(1)

r = RP.REFPROPdll('AMMONIA;WATER','PH','T',iUnit,iMass,iFlag,1500,h_1,{0.5 0.5}); 
o = double(r.Output); 
T_1 = o(1)

T_2 = 344.5068;

w = RP.REFPROPdll('AMMONIA;WATER','PT','H',iUnit,iMass,iFlag,1600,T_2,{0.5 0.5}); 
g = double(w.Output); 
h_2 = g(1)

d = RP.REFPROPdll('AMMONIA;WATER','PH','T',iUnit,iMass,iFlag,1600,h_2,{0.5 0.5}); 
z = double(d.Output); 
T_3 = z(1)

%% deactivate these lines to alter results   %%

s = RP.REFPROPdll('Water','PQ','T',iUnit,iMass,iFlag,101.325,0.0,{1.0});
k = double(s.Output);
T_4 = k(1)

I’m now setting the Temperature T_0 to be 328.9498 K. With this, and a pressure of 1500 kPa I’m calculating the corresponding enthalpy and afterwards again the temperature. So T_1 should equal T_0. The second state point of the ammonia/water mixture with given T_2 = 344.5068 K is just an example of a state point where the mentioned behaviour does not occur. With a first run of the code, the results are as follows:

h_1 = 170.02 kJ/kg
T_1 = 328.94 K
h_2 = 250.01 kJ/kg
T_3 = 344.50 K
T_4 = 373.12 K 

So this is matching the expactations, where T_1 equals T_0 and T_3 equals T_2, as well as all of these matching the results of the GUI mentioned in my first post to this thread. However, when I run the exact same code - without any changes - a second time, the results of h_1 and T_1 will change again to be:

h_1 = 164.29 kJ/kg
T_1 = 327.83 K
h_2 = 250.01 kJ/kg
T_3 = 344.50 K
T_4 = 373.12 K

So now T_1 and T_0 don’t match. Which would be a bug as of your definition. Yet T_2 and T_3 still match as the result of the calculation of h_2 has not changed. Keep in mind, up until this point I haven’t changed anything about the code. I’ve just run the same code a second time. And the results are changing. Which seems very odd, because the same code should always yield the same results no matter how often you run it?

I can recreate this whole behaviour again by deactivating the mentioned lines of code which calculates the properties of water. (As described in my first post)

Of course reference states for ethalpy may be different when comparing the GUI and the DLL, but should be consistent in either case itself? What I don’t understand is, that how can the exact same code yield two different results when run twice? When I don’t change anything about the code. The used reference state for the enthalpy shouldn’t change during the first run of the code, right? Otherwise calculated differences in enthalpy would become meaningless as well.

[ianhbell]

Smells like a bug to me. @EricLemmon will have to dig into this one. Especially fishy is the fact that just running the same calculation again results in this behavior, which seems like a caching issue.

[EricLemmon]

It is possible that you have another ammonia.fld, water.fld, or hmx.bnc file on your machine that is not from version 10. Can you do a full system search looking for each of these files (also look for Refprop.dll) and get rid of anything that is not in your Refprop directory? Be sure to search Windows, System, and hidden files to find them all.

[jastema]

Hello Eric, sorry it took me so long to answer.

I've tried your suggested fix and I'm pretty sure I deleted all of the mentioned files which haven't been in the Refprop folder. Yet this doesn't change anything concerning the behaviour of this code.

Could there be another solution to this problem ?

Thanks a lot

[BenMarkmann]

Hi Eric & Ian, I have also run into the problem that jastema described above. Actually, we're working on the same topic. Just tried the beta 10.0.0.61. But the results still differ. Have you released any new version that might fix the problem? All the best! Benjamin

[EricLemmon]

I just tried this directly in the Fortran code and did not find any problems. The fortran code is below, perhaps I implemented it differently than you did? I then did the same thing in VB6 to test the calls to the DLL and got the same thing as with the Fortran calls. I tried both 10.0 and 10.0.0.61. The code is also shown below. We may need Ian to test the calls from the Python wrapper to see if he finds a difference.

 1    hIn='MKS'
      call GETENUM(0,hIn,iUnit,ierr,herr)
      iMass = 1
      iFlag = 1

      T_0 = 328.9438
      z(1:2)=0.5d0

      call REFPROP ('ammonia;water','pt','h',iUnit,iMass,iFlag,
     &1500d0,T_0,z,Ou,hfle,iUnitx,x,y,x3,q,ierr,herr)
      h_1 = ou(1)

      call REFPROP ('ammonia;water','ph','t',iUnit,iMass,iFlag,
     &1500d0,h_1,z,Ou,hfle,iUnitx,x,y,x3,q,ierr,herr)

      T_1 = ou(1)
      T_2 = 344.5068d0
      write (*,*) T_1

      call REFPROP ('ammonia;water','pt','h',iUnit,iMass,iFlag,
     &1600d0,T_2,z,Ou,hfle,iUnitx,x,y,x3,q,ierr,herr)
      h_2 = ou(1)

      call REFPROP ('ammonia;water','ph','t',iUnit,iMass,iFlag,
     &1600d0,h_2,z,Ou,hfle,iUnitx,x,y,x3,q,ierr,herr)
      T_3 = ou(1)
      write (*,*) T_3
c %% deactivate these lines to alter results   %%

      z(1)=1
      z(2)=0
      call REFPROP ('water','pq','t',iUnit,iMass,iFlag,
     &.101325d0,0d0,z,Ou,hfle,iUnit,x,y,x3,q,ierr,herr)
      T_4 = ou(1)
      write (*,*) T_4
      goto 1

----------------VB6 code:

      Call REFPROP("", "", "DLL#", iUnit, iMass, iFlag, P, T, z, Output, hUnits, q, ierr, herr)

1     iMass = 1
      iFlag = 1
      iUnit = 7

      T = 328.9438
      P = 1500
      z(1) = 0.5
      z(2) = 0.5
      Call REFPROP("ammonia;water", "PT", "H", iUnit, iMass, iFlag, P, T, z, Output, hUnits, q, ierr, herr)
      h = Output(1)
      Call REFPROP("ammonia;water", "PH", "T", iUnit, iMass, iFlag, P, h, z, Output, hUnits, q, ierr, herr)
      T = Output(1)

      T = 344.5068
      P = 1600
      Call REFPROP("ammonia;water", "PT", "H", iUnit, iMass, iFlag, P, T, z, Output, hUnits, q, ierr, herr)
      h = Output(1)
      Call REFPROP("ammonia;water", "PH", "T", iUnit, iMass, iFlag, P, h, z, Output, hUnits, q, ierr, herr)
      T = Output(1)

      z(1) = 1
      z(2) = 0
      P = 0.101325
      qq = 0
      Call REFPROP("water", "PQ", "T", iUnit, iMass, iFlag, P, q, z, Output, hUnits, q, ierr, herr)
      T = Output(1)

      GoTo 1

[ianhbell]

We have identified a serious bug in the REFPROPdll function when iMass=1 and the mass fractions are passed as an array. It seems that this bug does not occur when the compositions are passed in as part of the mixture composition string. We want to make sure we fix this properly, so we are not rushing a fix out the door. In short, to workaround this bug, use mole fractions, or if that is not possible, specify the mass fractions in the fluid string. An alternative solution is to call REFPROP via CoolProp (http://www.coolprop.org/coolprop/REFPROP.html), as the mass->mole conversions are done correctly internally in CoolProp.

[ianhbell]

I should have been more clear. The high-level interface in CoolProp ALWAYS uses mole fractions, but the low-level interface gives you options:

import CoolProp.CoolProp as CP
root = 'c:/Program Files (x86)/REFPROP10'
CP.set_config_string(CP.ALTERNATIVE_REFPROP_PATH, root)
AS = CP.AbstractState('REFPROP', 'Ammonia&Water')
AS.set_mass_fractions([0.5,0.5])
AS.update(CP.HmassP_INPUTS, 170e3, 1500e3)
print(AS.T())

yielding

328.94376677674626

[EricLemmon]

This is definitely a bug that I did not see in my previous code that I pasted above because I did not reset the value of z to the 50/50 mass fractions as you did for every single call. In my case what happened is that the 50/50 mass fractions sent in were converted to mole fractions for use in the code, but were also returned on a mole basis. Because I did not change it back to 50/50 as you did between each call, the final calculation of T returned to the original value that I sent in. However, the enthalpy that came back was incorrect (it was for the 0.5/0.5 mole fraction state).

(Refprop tip: In the GUI, if you are using mass fractions, table titles and so on will show "50/50", but if you are using mole fractions, they will show "0.5/0.5")

The best way to avoid this issue completely is to first call the Refprop subroutine with just the names of the fluid, like this (in VB6):

      iMass = 1
      iFlag = 1
      iUnit = 7
      z(1) = 0.5
      z(2) = 0.5
      Call REFPROP("ammonia;water", "PT", "H", iUnit, iMass, iFlag, 0, 0, z, Output, hUnits, q, ierr, herr)

The values of z only need to be sent if you want the SATSPLN routine to be called, but realize that the output values for z will have been incorrectly changed. If iFlag is 0, you do not need to include z on this line.

Then on subsequent calls do not include the fluid names, just calculate the properties like this:

      z(1) = 0.5
      z(2) = 0.5
      T = 328.9438
      P = 1500
      Call REFPROP("", "PT", "H", iUnit, iMass, iFlag, P, T, z, Output, hUnits, q, ierr, herr)
      h = Output(1)
      Call REFPROP("", "PH", "T", iUnit, iMass, iFlag, P, h, z, Output, hUnits, q, ierr, herr)
      T = Output(1)

It doesn't matter if you set z again before the second call, either way will produce the same results.

When switching to water, it is not necessary to do as above, you can include the fluid name with the calculation requests, like this;

      P = 0.101325
      q = 0
      Call REFPROP("water", "PQ", "T", iUnit, iMass, iFlag, P, q, z, Output, hUnits, qq, ierr, herr)

But when switching back to the mixture, do the same procedure as I showed above.

We will get this bug fixed, but it will need some time for testing because such a change could impact other calculations and we need to be sure all is working before releasing an update. Having a "one routine does all" in Version 10 has been great, but our testing obviously did not hit all possible scenarios.

[Puya-Javidmand]

I have a similar issue when using python wrapper in python. When using mixtures, when I call the wrapper for the second time it gives different results or actually errors (as you can see in the second result below) REFPROPdlloutput(z=array('d', [0.697614699375863, 0.302385300624138, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]), Output=array('d', [104.10321447468603, -999

REFPROPdlloutput(z=array('d', [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]), Output=array('d', [-9999990.0, -9999990.0, -9999990.0, -9999990.0, -999999

[ianhbell]

@japuay It is usually recommended to open a new issue and link it to this one rather than adding to it unless you can be 100% sure you have the same problem. Also, please provide the code that generated these outputs so we have a possibility to debug what happened.

[Puya-Javidmand]

import statistics
import importlib
import sys
from openpyxl.chart import BarChart, Reference, LineChart
import glob
import pandas
import numpy as np
import os; os.environ['RPPREFIX'] = r'C:\Users\pjavidmand\Desktop\TRA\P test'
from ctREFPROP.ctREFPROP import REFPROPFunctionLibrary

RP = REFPROPFunctionLibrary(os.environ['RPPREFIX'])
RP.SETPATHdll(os.environ['RPPREFIX'])

p = 351.290828571429
Q = 1.0

print (RP.REFPROPdll("R410A","PQ","T",5,0,0,p,Q,[1.0]))
print (RP.REFPROPdll("R410A","PQ","T",5,0,0,p,Q,[1.0]))

[ianhbell]

Confirmed on my side with the code below. Some notes for you:

• As noted in the docs, you should not hard code the unit system, see code below for recommended solution
• You should always specify either R410A.MIX or R410A.PPF, otherwise it is ambiguous. FWIW, the pseudo-pure works fine; this is a composition caching bug. Personally I think REFPROP should not allow R410A without suffix, but we allow it for historical reasons.

import os; os.environ['RPPREFIX'] = r'D:/Program Files (x86)/REFPROP'
from ctREFPROP.ctREFPROP import REFPROPFunctionLibrary

RP = REFPROPFunctionLibrary(os.environ['RPPREFIX'])
RP.SETPATHdll(os.environ['RPPREFIX'])

p = 351.290828571429
Q = 1.0

print (RP.REFPROPdll("R410A","PQ","T",RP.ENGLISH,0,0,p,Q,[1.0]))
print (RP.REFPROPdll("R410A","PQ","T",RP.ENGLISH,0,0,p,Q,[1.0]))

[Puya-Javidmand]

I used your code, but it still has the same issue in my side

[Puya-Javidmand]

Ok, PPF is working. Mix does not work Thank you

[ianhbell]

That's correct, this is a bug in REFPROP. My comments were more best-practices recommendations than a workaround.

[Puya-Javidmand]

When I am using it for R454B I still have same problem even with .PPF

[ianhbell]

Correct, R454B.PPF should not work at all, since that we don't have a PPF for that predefined mixture.

[Puya-Javidmand]

Yes there is no PPF, but I have the same problem with .mix as well

[Puya-Javidmand]

Do you suggest any solution for R454B?

[ianhbell]

Yes, the workaround is to use the molar composition of the mixture directly, rather than specifying the predefined mixture file.

[Puya-Javidmand]

Can you do it on this command for me please: print (RP.REFPROPdll("R454B.PPF","PQ","T",RP.ENGLISH,0,0,p,Q,[1.0]).Output[0])

[ianhbell]

Open predefined mixture file:

R454B    [R-32/1234yf (68.9/31.1)]
 62.61363340943545 351.2537406768218 5266.938468522863 7.075660535893995
 2
R32.FLD
R1234YF.FLD
 .829247912869081
 .170752087130919
 0

Components are R32 and R1234yf, mole fractions are [0.829247912869081, 0.170752087130919]

[ianhbell]

So:

print (RP.REFPROPdll("R32*R1234yf","PQ","T",RP.ENGLISH,0,0,p,Q,[0.829247912869081, 0.170752087130919]))

[Puya-Javidmand]

Thank you. It is working.