usnistgov / REFPROP-issues

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

ABFLSH requires molar inputs despite iFlag/iMass settings #293

Open HCSpin opened 4 years ago

HCSpin commented 4 years ago

When using ABFLSHdll with inputs such as T and H, to get proper output values for D, for example, H must be entered in molar units. This contradicts what is written in issue 260: "Only the very top level routines (REFPROP, ALLPROPS, and ABFLSH) allow the user to set the input and output units." It also contradicts the guidance in the REFPROP Documentation, Release 10: image

I am using REFPROP version 10.002 on MacOSX 10.15.4, with python.

For example, with Nitrogen at 300K and 101.325 kPa, (MKS units) ABFLSHdll("TP",300,101.325,[1],22) yields Density = 1.138 kg/m^3 and Enthalpy = 433.2 J/g

and

ABFLSHdll("TP",300,101.325,[1],20) yields Density = 0.0406 mol/liter and Enthalpy = 12135.8 J/mol

Both of these above calculations are consistent with the iFlag/iMass guidance in the manual for this function. The correct output units are given.

However, running the reverse calculation with: ABFLSHdll("TH",300,433.2,[1],22) yields the following incorrect values: Density = 159.8 kg/m^3 Pressure = 14560.9 kPa and Enthalpy = 407.2 J/g

Changing the value of the enthalpy input to the molar form: ABFLSHdll("TH",300,12135.8,[1],20) yields proper molar output values: Density = 0.0406 mol/liter Pressure = 101.325 kPa and Enthalpy = 12135.8 J/mol

and, proper mass-based output values:

ABFLSH("TH",300,12135.8,[1],22) Density = 1.138 kg/m^3 Pressure = 101.325 kPa and Enthalpy = 433.2 J/g

It appears that the language in the documentation should refer only to output values of this function and state that all inputs must be in the molar form. Is that the desired behavior?

EricLemmon commented 4 years ago

Take a look at these: https://github.com/usnistgov/REFPROP-wrappers/issues/143 https://github.com/usnistgov/REFPROP-wrappers/issues/272 https://github.com/usnistgov/REFPROP-wrappers/issues/256

Please send one of us a private email and we'll send you the updated DLL:

eric.lemmon@nist.gov or ian.bell@nist.gov

ianhbell commented 4 years ago

Notes:

  1. TH are possibly multi-valued inputs, so you might get a different solution than you expect!
  2. Never hardcode unit system codes, always request them from GETENUMdll. The enumerated values might change in the future.
  3. Please provide minimal runnable examples for easier debugging on our side.
HCSpin commented 4 years ago

Thanks for the feedback.
Point 1: Could you explain what you mean by TH being possibly multi-valued inputs? Are you referring to 2-phase states? (I intentionally chose this simple vapor-phase example to familiarize myself with the code.) Point 2: You are referring to the units used (i.e., "MKS" should be replaced by GETENUMdll(0,"MKS").iEnum) and are not suggesting that the hardcoding of iFlag as 20 or 22 is improper, correct? Point 3: Would that code be helpful now, or is it too late?

ianhbell commented 4 years ago

1) Here's an example for propane at 400 K:

import numpy as np
import CoolProp.CoolProp as CP
import matplotlib.pyplot as plt

rhomolar = np.geomspace(1e-10, 10000, 1000)
hmolar = CP.PropsSI('Hmolar','T',400,'D',rhomolar,'REFPROP::Propane')

plt.plot(rhomolar, hmolar)
plt.xlabel(r'$\rho$ / mol/m$^3$')
plt.ylabel(r'$h$ / J/mol')
plt.savefig('TH.png')
plt.tight_layout(pad=0.2)
plt.show()

TH

2) Yes, that's right, as is clearly spelled out in the docs 3) Never too late, and useful for others that might land here in the future

HCSpin commented 4 years ago

Thank you for the example.

Here is the code used:

import  os
import  numpy as np
from    ctREFPROP.ctREFPROP import REFPROPFunctionLibrary

RP =                REFPROPFunctionLibrary(os.environ['DYLD_LIBRARY_PATH'])
RP.SETPATHdll(os.environ['DYLD_LIBRARY_PATH'])
print("REFPROP version: ",RP.RPVersion())

Units=              RP.GETENUMdll(0,"MKS").iEnum

#set reference
RP.REFPROPdll("Nitrogen","SETREF","NBP",Units,0,0,0,0,[1])

T=              300     #K
P=              101.325 #kPa

r =                 RP.REFPROPdll("Nitrogen","TP","D",Units, 0,0,T,P, [1])
D=              r.Output[0]
units_D=            r.hUnits
r =                 RP.REFPROPdll("Nitrogen","TP","H",Units, 0,0,T,P, [1])
H=              r.Output[0]
units_H=            r.hUnits

print()
print("------------------------------")
print("Original Density and Enthalpy using REFPROP")
print("    Density: {:.3f} {:s}".format(D,units_D))
print("    Enthalpy: {:.3f} {:s}".format(H,units_H))

#initial calculation of density,enthalpy using T and P and iFlag signifying mass units
iFlag=              22      #vapor phase (tens digit); mass units except for composition (ones digit)
D_mass=             RP.ABFLSHdll("TP",T,P,[1],iFlag).D
H_mass=             RP.ABFLSHdll("TP",T,P,[1],iFlag).h
print("------------------------------")
print("ABFLSHdll calculations with mass-based i/o (T,P inputs):")
print("    Density: {:.3f} {:s}".format(D_mass,units_D))
print("    Enthalpy: {:,.1f} {:s}".format(H_mass,units_H))

#initial calculation of density,enthalpy using T and P and iFlag signifying molar units
iFlag=              20      #vapor phase (tens digit); molar units  (ones digit)
D_molar=            RP.ABFLSHdll("TP",T,P,[1],iFlag).D
H_molar=            RP.ABFLSHdll("TP",T,P,[1],iFlag).h
print()
print("ABFLSHdll calculations with molar-based i/o (T,P inputs):")
print("    Density: {:.4f} {:s}".format(D_molar,'mol/liter'))
print("    Enthalpy: {:,.1f} {:s}".format(H_molar,'J/mol'))

print("------------------------------")
print()
print("Reverse calculations using T and H (MASS units) as inputs:")
print("    (iFlag indicates mass-based inputs expected)")
iFlag=              22      #vapor phase (tens digit); mass units except for composition (ones digit)
D_mass_reverse1=    RP.ABFLSHdll("TH",T,H_mass,[1],iFlag).D
P_mass_reverse1=    RP.ABFLSHdll("TH",T,H_mass,[1],iFlag).P
H_mass_reverse1=    RP.ABFLSHdll("TH",T,H_mass,[1],iFlag).h
print()
print("ABFLSHdll calculations with mass-based i/o expected:")
print("    Density (reverse calc): {:.3f} {:s}".format(D_mass_reverse1,units_D))
print("    Pressure (reverse calc): {:.3f} {:s}".format(P_mass_reverse1,"kPa"))
print("    Enthalpy (reverse calc): {:,.1f} {:s}".format(H_mass_reverse1,units_H))
print("-->These results do not match original outputs")

print()
print("Reverse calculations using T and H (MOLAR units) as inputs:")
print("    (iFlag still indicates mass-based inputs expected)")
iFlag=              22      #vapor phase (tens digit); mass units except for composition (ones digit)
D_mass_reverse2=    RP.ABFLSHdll("TH",T,H_molar,[1],iFlag).D
P_mass_reverse2=    RP.ABFLSHdll("TH",T,H_molar,[1],iFlag).P
H_mass_reverse2=    RP.ABFLSHdll("TH",T,H_molar,[1],iFlag).h
print()
print("ABFLSHdll calculations with mass-based i/o expected:")
print("    Density (reverse calc): {:.3f} {:s}".format(D_mass_reverse2,units_D))
print("    Pressure (reverse calc): {:.3f} {:s}".format(P_mass_reverse2,"kPa"))
print("    Enthalpy (reverse calc): {:,.1f} {:s}".format(H_mass_reverse2,units_H))
print("-->These correctly match original outputs (despite mismatch between iFlag and input units)")

print()
print("Change iFlag to indicate that molar inputs are expected:")
print("Reverse calculations using T and H (MOLAR units) as inputs:")
iFlag=              20      #vapor phase (tens digit); molar units except for composition (ones digit)
D_molar_reverse=    RP.ABFLSHdll("TH",T,H_molar,[1],iFlag).D
P_molar_reverse=    RP.ABFLSHdll("TH",T,H_molar,[1],iFlag).P
H_molar_reverse=    RP.ABFLSHdll("TH",T,H_molar,[1],iFlag).h
print()
print("ABFLSHdll calculations with molar-based i/o expected:")
print("    Density (reverse calc): {:.3f} {:s}".format(D_molar_reverse,'mol/liter'))
print("    Pressure (reverse calc): {:.3f} {:s}".format(P_molar_reverse,"kPa"))
print("    Enthalpy (reverse calc): {:,.1f} {:s}".format(H_molar_reverse,'J/mol'))
print("-->These correctly match original outputs")

and the output: image

choeges commented 4 years ago

Hello, when using the ABFLSHdll using a mixture as fluid, I receive a similar problem using the Python wrapper. I tried 3 different method-calls your wrapper provides and compared them with values I got from the general REFPROP-GUI: REFPROPdll, TQFLSHdll, ABFLSHdll When using REFPROPdll and TQFLSHdll I get the expected behaviour. When using the ABFLSHdll routine, I receive different results for flag = 0 (molar based) and flag 1 (mass based) - see the code below. As far as I checked, the routine works fine with both flags for pure substances (at least in Version 10.0077 - I had some problems there using 10.0002 as well though).

Versions:

Code used:

import os
from ctREFPROP.ctREFPROP import REFPROPFunctionLibrary
if __name__ == "__main__":
    cwd = os.getcwd()
    root = os.environ["RPPREFIX"]
    RP = REFPROPFunctionLibrary(cwd+r"\\REFPRP64.DLL")
    RP.SETPATHdll(root)
    MOLAR_BASE_SI = RP.GETENUMdll(0, "MOLAR BASE SI").iEnum
    print("Version:", RP.RPVersion())
    mix = 'R407C.mix'
    sm = RP.SETMIXdll(mix, 'HMX.BNC', 'DEF')
    RP.SETREFdll('DEF', 2, sm.z, 0, 0, 0, 0)
    # One state for R407C received by REFPROP GUI
    T = 267.15
    P = 371450
    D = 15.964
    u = 383270
    h = 406530
    s = 1784
    Q = 1.0

    rpdll = RP.REFPROPdll("", "QT", "H;S;M", MOLAR_BASE_SI, 0, 0, Q, T, sm.z)
    Hmolar, Smolar, M = rpdll.Output[0:3]
    print("REFPROPdll", Hmolar / M, Smolar / M)

    tqflsh = RP.TQFLSHdll(T, Q, sm.z, 0)
    print("TQFLSHdll", tqflsh.h / M, tqflsh.s / M)

    abflsh_molar = RP.ABFLSHdll("TQ", T, Q, sm.z, 0)
    print("ABFLSHdll molar", abflsh_molar.h/M, abflsh_molar.s/M)
    abflsh_mass = RP.ABFLSHdll("TQ", T, Q, sm.z, 1)
    print("ABFLSHdll mass", abflsh_mass.h*1000, abflsh_mass.s*1000)

    print("RefProp GUI", h, s)

The results I receive:

Version: 10.0077
REFPROPdll 406534.67572704994 1784.0008962828858
TQFLSHdll 406534.67572704994 1784.0008962828858
ABFLSHdll molar 406534.67572704994 1784.0008962828858
ABFLSHdll mass 460051.10091473174 2018.8476540162453
RefProp GUI 406530 1784

Like I said, ABFLSHdll Routine works fine when using molar base (flag = 0) but I cannot really comprehend how the displayed error occurs.

EricLemmon commented 4 years ago

Can you send me a personal email at eric.lemmon@nist.gov (ignore my out-of-office message)? We are nearing a solution to this, but it has been a huge amount of work and will take quite a bit of testing.

EricLemmon commented 4 years ago

Could you send me your email address again? I can't find a response from you before (far too many emails to sort through).