usnistgov / REFPROP-issues

A repository solely used for reporting issues with NIST REFPROP
26 stars 13 forks source link

MELT-PT does not work, nor SUBL-XX #362

Open ianhbell opened 3 years ago

ianhbell commented 3 years ago

The output of MELT-PT does not give sensible values. I passed in an invalid input for the state variables, and even when I supply a valid pressure, I never get a reasonable melting temperature. How is this supposed to work in REFPROP function? Should this be removed from the docs?

import os
from ctREFPROP.ctREFPROP import REFPROPFunctionLibrary
os.environ['RPPREFIX'] = 'D:/Program Files (x86)/REFPROPalpha'
root = os.environ['RPPREFIX']
RP = REFPROPFunctionLibrary(root)
print(RP.RPVersion())

r2 = RP.REFPROPdll('NITROGEN','XYZ','MELT-PT', RP.MASS_BASE_SI, 0, 0, 0.0, 0.0, [1.0])
print(r2.Output[0], r2.herr)

yields

10.0086
63.151 [MELTP error 143] Input value equal to or less than zero.
ianhbell commented 2 years ago

Let's just remove this from the docs. There are already other ways of getting at this information that are clear enough.

ianhbell commented 2 years ago

Also, if we are going to fix this, SUBL-TP doesn't work either

EricLemmon commented 2 years ago

I've been looking into this one, again..., and am still confused as to how it ever worked in the first place, if it did at all. The logic doesn't make sense. I abandoned the problem several months back because I did not see a good solution. But over the last few days as we were fixing a different problem dealing with the melting line, I wondered how those inputs to the REFPROP subroutine specified the calculation of a melting line state. Upon looking through the code, I realized it used a much more brilliant method! With MELT-TP, the idea (that doesn't work) is to use something like this:

      T=REFPROPF ('xyz','MELT-PT',1,0,23.45d0,0d0,z,q,ierr,herr)      

(where the xyz mimics what Ian did above and the 23.45 is an arbitrary number to aid the eye to see the pressure input.)

A different option is available to make the same calculation, but in a much better way:

      T=REFPROPF ('Pmelt','T',1,0,23.45d0,0d0,z,q,ierr,herr) 

This would be identical to using TQ, Tliq, or Tvap to obtain a saturation property. With this method you are not limited to just T as the output as is the case with MELT-PT, but you can obtain anything along the melting line (actually the freezing line, the properties returned are all in the liquid phase, not the solid phase). For example, the following command will return T, D, H, and S:

      call REFPROP (' ','Pmelt','T,D,H,S',1,0,0,23.45d0,0.d0,z,Output,hU,iU,x,y,x3,q,ierr,herr)

A little back and forth testing with inputs and outputs reveals that all is consistent:

      P=23.45d0
      T =REFPROPF ('Pmelt','T',2,0,P,0d0,z,q,ierr1,herr)      
      D =REFPROPF ('Pmelt','D',2,0,P,0d0,z,q,ierr2,herr)      
      P1=REFPROPF ('TD',   'P',2,0,T,D,  z,q,ierr3,herr)      
      D1=REFPROPF ('Tmelt','D',2,0,T,0d0,z,q,ierr4,herr)      
      H1=REFPROPF ('Tmelt','H',2,0,T,0d0,z,q,ierr5,herr)      
      T1=REFPROPF ('PH'   ,'T',2,0,P,H1 ,z,q,ierr6,herr)      

The outputs show that both T and T1 are almost identical (within the tolerance of the iteration scheme) and D and D1 are almost identical.

These methods would allow us to get rid of the MELT-TP, MELT-PT, and so on options. If we were to try to fix the code to make them work, many doors leading to new bugs would be opened.

For sublimation properties, everything is identical to the above commands aside from changing melt to subl:

      T=REFPROPF ('Psubl','T',1,0,23.45d0,0d0,z,q,ierr,herr) 

For the ancillary equations for liquid-vapor saturation states, I suggest we get rid of ANC-TP, ANC-TDL, ANC-TDV, and so on and instead use this same method, but with melt changed to anc:

      call REFPROP (' ','Panc','T,D,H,S',1,0,1,500d0,10.d0,z,Output,hU,iU,x,y,x3,q,ierr,herr)

In most cases the user should never use anc as an input, rather TQ with Q=0 or 1 should be used. But there are a few situations where this need might occur. (The ancillary equations are simplified equations dedicated to just saturated vapor pressure, liquid density, and vapor density states. For more information, see: https://trc.nist.gov/refprop/FAQ/R125.PDF)

Ian, if you are in agreement, I'll go ahead and make the changes.

ianhbell commented 2 years ago

I think I'm ok with this. It's a breaking change and needs to be documented since it worked in 10.0. I don't love the Panc thing though, where does the Q go there to specify which phase you mean? For subl and melt it is straightforward because you are always clear about what phase you mean (vapor and liquid, respectively), but that's not true for ancillaries. Your example uses a quality of 10, right?

EricLemmon commented 2 years ago

It would appear that we have good justification to remove access to these functions anyway. The routines from Version 10.0 are not converting the input temperature or pressure to the proper value based on the input unit system. For example, requesting the melting pressure at 68 C (which would be a ridiculously high pressure) returns the melting pressure at 68 K, which is the base unit system in Refprop.

      call GETENUM (0,'SIwithC',iU,ierr,herr)
      T=68d0
      P=REFPROPF('xyz','Melt-TP',iU,0,T,0d0,z,q,ierr1,herr)
      C=REFPROPF('xyz','Melt-PT',iU,0,0d0,P,z,q,ierr2,herr)

The first call to the REFPROP routine returns an incorrect pressure of 22.695 MPa rather than the correct value of 3119 MPa (for 68 C). With the incorrect pressure as input in the last line, the routine returns a temperature of -209.9967 C (which is the correct melting temperature in celsius for P=22.695 MPa).

Requests for sublimation or ancillary values are plagued with the same problems as above.

I'll go ahead and remove access to these routines for Melt-TP, Melt-PT, etc., with an error message to use Tmelt, Pmelt, and so on instead.

EricLemmon commented 2 years ago

For the ancillaries, I think the best thing to do would be to put the Panc, Tanc, and Danc strings into the requested output string. For example, given a temperature of 80 K for nitrogen, you could request both the vapor pressure from the equation of state as normal and in addition the vapor pressure from the ancillary equation like this:

      call REFPROP (' ','TQ','P,Panc',1,0,0,80d0,1d0,z,Output,hU,iU,x,y,x3,q,ierr,herr)

For densities, we could simply use Danc and let the value of Q decide which ancillary to call (DLSATT vs. DVSATT). Or we could use DLanc and DVanc so that both could be returned (which would mimic that done currently with the additional characters "liq" and "vap"). For example, given T, the densities from the equation of state and the ancillaries could be obtained with:

      call REFPROP (' ','Tsat','Dliq,Dvap,DLanc,DVanc',1,0,0,80d0,0.d0,z,Output,hU,iU,x,y,x3,q,ierr,herr)

I prefer the second option, the use of Q in the first option could be confusing when Q>0 and Q<1 (but not equal to 0 or 1).