search

usnistgov / REFPROP-issues

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

Large Error When Calculating Pressure from Density and Temperature for Sat. Liquid Propane - E.W. Lemmon, M.O. McLinden, W. Wagner, J. Chem. Eng. Data 54, 3141 (2009) #568

Closed mattybranham closed 1 year ago

mattybranham commented 1 year ago

This question is not explicitly about REFPROP, but about correctly implementing an equation of state. I am an active REFPROP user, but I am running a model that needs to find fluid properties a few billion times per model run. In order to make that tenable, I am building a home-brewed fluid properties calculator specific to the refrigerant in question - propane.

In E.W. Lemmon, M.O. McLinden, W. Wagner, J. Chem. Eng. Data 54, 3141 (2009) - "Thermodynamic Properties of Propane. III. A Reference Equation of State for Temperatures from the Melting Line to 650 K and Pressures up to 1000 MPa" - the pressure of propane at a given state can be found using Eq. 33 and Eq. 50. After implementing those equations in code, however, I find large discrepancies between my calculated results and the tabular values for propane's pressure in the saturated liquid state. The discrepancies are greatest at low temperatures. The results for saturated vapor seem consistent with the provided data.

I am assuming that the equation of state applies to both liquid and vapor states. If that's true, does the equation of state require a different implementation for liquid vs. vapor states? Or is there another obvious detail that I am missing that would lead to correct results for saturated vapor and incorrect results for the saturated liquid?

In case it would be helpful, I have attached the Python code used to run Eq. 33 and Eq. 50 below. The code also generates a figure showing the magnitude of the error between calculated results and Table A1.

I would appreciate any guidance, I have been banging my head against this for a while.

Python 3.9.7

import numpy as np
import pandas as pd
import os
from matplotlib import pyplot as plt

#Propane test
molarMass=44.09562      #g/mol
RPropane=8.314472       #J/(mol*K)
rhoCritical=5.00        #mol/dm^-3
tempCritical=369.89     #K
pressureCritical=4.2512 #MPa

#Calculate saturated vapor density
satVaporDensityN={'N1':-2.4887,'N2':-5.1069,'N3':-12.174,'N4':-30.495,'N5':-52.192,'N6':-134.89}    #Eq. 7
satLiquidDensityN={'N1':1.82205,'N2':0.65802,'N3':0.21109,'N4':0.083973}                            #Eq. 6

#Propane constants (Table 4)
NkVals=[0.042910051,1.7313671,-2.4516524,0.34157466,-0.46047898,-0.66847295,0.20889705,0.19421381,\
        -0.22917851,-0.60405866,0.066680654,0.017534618,0.33874242,0.22228777,-0.23219062,-0.09220694,\
            -0.47575718,-0.017486824]                       #Table 4
Nk=pd.DataFrame.from_dict(dict(zip(range(1,19),NkVals)),orient='index',columns=['Nk'])
tkVals=[1,0.33,0.8,0.43,0.9,2.46,2.09,0.88,1.09,3.25,4.62,0.76,2.5,2.75,3.05,2.55,8.4,6.75] #Table 4
tk=pd.DataFrame.from_dict(dict(zip(range(1,19),tkVals)),orient='index',columns=['tk'])
dkVals=[4,1,1,2,2,1,3,6,6,2,3,1,1,1,2,2,4,1]                #Table 4
dk=pd.DataFrame.from_dict(dict(zip(range(1,19),dkVals)),orient='index',columns=['dk'])
lkVals=[1,1,1,1,2,2]
lk=pd.DataFrame.from_dict(dict(zip(range(6,12),lkVals)),orient='index',columns=['lk'])
etakVals=[0.963,1.977,1.917,2.307,2.546,3.28,14.6]          #Table 4
etak=pd.DataFrame.from_dict(dict(zip(range(12,19),etakVals)),orient='index',columns=['etak'])
betakVals=[2.33,3.47,3.15,3.19,0.92,18.8,547.8]             #Table 4
betak=pd.DataFrame.from_dict(dict(zip(range(12,19),betakVals)),orient='index',columns=['betak'])
gammakVals=[0.684,0.829,1.419,0.817,1.5,1.426,1.093]        #Table 4
gammak=pd.DataFrame.from_dict(dict(zip(range(12,19),gammakVals)),orient='index',columns=['gammak'])
epsilonkVals=[1.283,0.6936,0.788,0.473,0.8577,0.271,0.948]  #Table 4
epsilonk=pd.DataFrame.from_dict(dict(zip(range(12,19),epsilonkVals)),orient='index',columns=['epsilonk'])
bkVals=[1.062478,3.344237,5.363757,11.762957]               #Equation 22
bk=pd.DataFrame.from_dict(dict(zip(range(3,7),bkVals)),orient='index',columns=['bk'])
vkVals=[3.043,5.874,9.337,7.922]                            #Equation 21
vk=pd.DataFrame.from_dict(dict(zip(range(3,7),vkVals)),orient='index',columns=['vk'])

#Actual saturation temperature, pressure, and density for propane from Table A1
temperatureArray=(np.concatenate([np.array([-187.625,-42.114]),np.arange(-185,100,5)]))+273.15   #K
temperatureArray.sort()
#Pressure in MPa
pressureArray=np.array([1.7203E-10,4.8526E-10,2.9386E-09,1.4596E-08,6.1256E-08,2.2252E-07,7.1368E-07,2.0542E-06,5.3795E-06,1.2965E-05,2.9040E-05,6.0950E-05,1.2073E-04,2.2708E-04,4.0774E-04,7.0215E-04,0.0011644,0.0018661,0.0028994,0.0043795,0.0064475,0.0092716,0.013049,0.018008,0.024404,0.032527,0.042693,0.055249,0.070569,0.089051,0.101325,0.11112,0.13723,0.16783,0.20343,0.24452,0.29162,0.34528,0.40604,0.47446,0.55112,0.6366,0.73151,0.83646,0.95207,1.079,1.2179,1.3694,1.5343,1.7133,1.9072,2.1168,2.343,2.5868,2.8493,3.1319,3.4361,3.7641,4.1195])
#Density in kg/m^-3
densityArray=np.array([733.13,1.07E-08,730.39,2.92E-08,725.2,1.67E-07,720.05,7.89E-07,714.92,3.15E-06,709.82,0.0000109,704.73,0.0000335,699.67,0.0000922,694.61,0.000232,689.56,0.000537,684.51,0.00116,679.46,0.00234,674.4,0.00447,669.33,0.00813,664.26,0.01413,659.16,0.02357,654.05,0.0379,648.91,0.05897,643.74,0.08904,638.55,0.13085,633.32,0.18762,628.06,0.26304,622.76,0.36132,617.41,0.48715,612.02,0.6457,606.57,0.84261,601.08,1.084,595.52,1.3764,589.9,1.727,584.2,2.143,580.88,2.4161,578.43,2.6326,572.58,3.2042,566.64,3.8669,560.6,4.6302,554.45,5.5046,548.19,6.5012,541.8,7.6321,535.27,8.9103,528.59,10.351,521.75,11.969,514.73,13.783,507.5,15.813,500.06,18.082,492.36,20.618,484.39,23.451,476.1,26.618,467.46,30.165,458.4,34.146,448.87,38.63,438.76,43.706,427.97,49.493,416.34,56.152,403.62,63.916,389.47,73.14,373.29,84.406,353.96,98.818,328.83,119,286.51,156.31])
densityArrayMolar=densityArray/molarMass        #mol*dm^-3

pressureCalc=np.zeros(densityArray.size,)
# cv=np.zeros(densityArray.size,)
# cp=np.zeros(densityArray.size,)
storeInd=0
for (density,temp) in zip(densityArrayMolar,np.repeat(temperatureArray,2)):
    Delta=density/rhoCritical
    Tau=tempCritical/temp
    #Calculate pressure
    #DeltadAlphaRdDelta is Equation 50
    DeltadAlphaRdDelta=np.sum(Nk.loc[1:5,'Nk']*Delta**dk.loc[1:5,'dk']*Tau**tk.loc[1:5,'tk']*dk.loc[1:5,'dk'])+\
        np.sum(Nk.loc[6:11,'Nk']*Delta**dk.loc[6:11,'dk']*Tau**tk.loc[6:11,'tk']*np.exp(-Delta**lk.loc[6:11,'lk'])*(dk.loc[6:11,'dk']-lk.loc[6:11,'lk']*Delta**lk.loc[6:11,'lk']))+\
            np.sum(Nk.loc[12:18,'Nk']*Delta**dk.loc[12:18,'dk']*Tau**tk.loc[12:18,'tk']*np.exp(-etak.loc[12:18,'etak']*((Delta-epsilonk.loc[12:18,'epsilonk'])**2)-(betak.loc[12:18,'betak']*(Tau-gammak.loc[12:18,'gammak'])**2))*(dk.loc[12:18,'dk']-2*etak.loc[12:18,'etak']*Delta*(Delta-epsilonk.loc[12:18,'epsilonk'])))   
    pressureCalc[storeInd,]=density*RPropane*temp*(1+DeltadAlphaRdDelta)*1000/1E6   #MPa, Equation 33
    # #Calculate cv
    # cvDifferentialIdeal=-3-Tau**2*np.sum(propaneConstants.loc[3:6,'vk']*propaneConstants.loc[3:6,'bk']**2*(np.exp(Tau*propaneConstants.loc[3:6,'bk'])/(np.exp(Tau*propaneConstants.loc[3:6,'bk'])-1)**2))                                 #Includes tau^2 term
    # cvDifferentialResidual=np.sum(propaneConstants.loc[1:5,'Nk']*Delta**propaneConstants.loc[1:5,'dk']*Tau**propaneConstants.loc[1:5,'tk']*(propaneConstants.loc[1:5,'tk']*(propaneConstants.loc[1:5,'tk']-1)))+\
    #     np.sum(propaneConstants.loc[6:11,'Nk']*Delta**propaneConstants.loc[6:11,'dk']*Tau**propaneConstants.loc[6:11,'tk']*np.exp(-Delta**propaneConstants.loc[6:11,'lk'])*(propaneConstants.loc[6:11,'tk']*(propaneConstants.loc[6:11,'tk']-1)))+\
    #         np.sum(propaneConstants.loc[12:18,'Nk']*Delta**propaneConstants.loc[12:18,'dk']*Tau**propaneConstants.loc[12:18,'tk']*np.exp(-propaneConstants.loc[12:18,'etak']*((Delta-propaneConstants.loc[12:18,'epsilonk'])**2)-(propaneConstants.loc[12:18,'betak']*(Tau-propaneConstants.loc[12:18,'gammak'])**2))*((propaneConstants.loc[12:18,'tk']-2*propaneConstants.loc[12:18,'betak']*Tau*(Tau-propaneConstants.loc[12:18,'gammak']))**2-propaneConstants.loc[12:18,'tk']-2*propaneConstants.loc[12:18,'betak']*Tau**2))  #Includes tau^2 term
    # cv[storeInd,]=-RPropane*(cvDifferentialIdeal+cvDifferentialResidual)     #J/mol*K
    # #Calculate cp
    # DeltaTaud2AlphaRdDeltadTau=np.sum(propaneConstants.loc[1:5,'Nk']*Delta**propaneConstants.loc[1:5,'dk']*Tau**propaneConstants.loc[1:5,'tk']*propaneConstants.loc[1:5,'tk']*propaneConstants.loc[1:5,'dk'])+\
    #     np.sum(propaneConstants.loc[6:11,'Nk']*Delta**propaneConstants.loc[6:11,'dk']*Tau**propaneConstants.loc[6:11,'tk']*np.exp(-Delta**propaneConstants.loc[6:11,'lk'])*(propaneConstants.loc[6:11,'tk']*(propaneConstants.loc[6:11,'dk']-propaneConstants.loc[6:11,'lk']*Delta**propaneConstants.loc[6:11,'lk'])))+\
    #         np.sum(propaneConstants.loc[12:18,'Nk']*Delta**propaneConstants.loc[12:18,'dk']*Tau**propaneConstants.loc[12:18,'tk']*np.exp(-propaneConstants.loc[12:18,'etak']*((Delta-propaneConstants.loc[12:18,'epsilonk'])**2)-(propaneConstants.loc[12:18,'betak']*(Tau-propaneConstants.loc[12:18,'gammak'])**2))*(propaneConstants.loc[12:18,'tk']-2*propaneConstants.loc[12:18,'betak']*Tau*(Tau-propaneConstants.loc[12:18,'gammak']))*(propaneConstants.loc[12:18,'dk']-2*propaneConstants.loc[12:18,'etak']*Delta*(Delta-propaneConstants.loc[12:18,'epsilonk'])))  #Includes tau*delta term
    # Delta2d2AlphaRdDelta2=np.sum(propaneConstants.loc[1:5,'Nk']*Delta**propaneConstants.loc[1:5,'dk']*Tau**propaneConstants.loc[1:5,'tk']*(propaneConstants.loc[1:5,'dk']*(propaneConstants.loc[1:5,'dk']-1)))+\
    #     np.sum(propaneConstants.loc[6:11,'Nk']*Delta**propaneConstants.loc[6:11,'dk']*Tau**propaneConstants.loc[6:11,'tk']*np.exp(-Delta**propaneConstants.loc[6:11,'lk'])*((propaneConstants.loc[6:11,'dk']-propaneConstants.loc[6:11,'lk']*Delta**propaneConstants.loc[6:11,'lk'])*(propaneConstants.loc[6:11,'dk']-1-propaneConstants.loc[6:11,'lk']*Delta**propaneConstants.loc[6:11,'lk'])-propaneConstants.loc[6:11,'lk']**2*Delta**propaneConstants.loc[6:11,'lk']))+\
    #         np.sum(propaneConstants.loc[12:18,'Nk']*Delta**propaneConstants.loc[12:18,'dk']*Tau**propaneConstants.loc[12:18,'tk']*np.exp(-propaneConstants.loc[12:18,'etak']*((Delta-propaneConstants.loc[12:18,'epsilonk'])**2)-(propaneConstants.loc[12:18,'betak']*(Tau-propaneConstants.loc[12:18,'gammak'])**2))*((propaneConstants.loc[12:18,'dk']-2*propaneConstants.loc[12:18,'etak']*Delta*(Delta-propaneConstants.loc[12:18,'epsilonk']))**2-propaneConstants.loc[12:18,'dk']-2*propaneConstants.loc[12:18,'etak']*Delta**2))  #Includes delta^2 term
    # cp[storeInd,]=cv[storeInd,]+RPropane*((1+DeltadAlphaRdDelta-DeltaTaud2AlphaRdDeltadTau)**2/(1+2*DeltadAlphaRdDelta+Delta2d2AlphaRdDelta2))
    storeInd+=1
# cvMass=cv/molarMass
# cpMass=cp/molarMass

pressureErrorSatLiquid=100*np.abs(pressureArray-pressureCalc[np.arange(0,pressureCalc.size,2)])/pressureArray
pressureErrorSatVapor=100*np.abs(pressureArray-pressureCalc[np.arange(1,pressureCalc.size,2)])/pressureArray
###Code below used to visualize accuracy of equations
propCheckFig2,propCheckAx2=plt.subplots()
propCheckAx2.scatter(temperatureArray,pressureErrorSatLiquid,label="Liquid",c="b")
propCheckAx2.scatter(temperatureArray,pressureErrorSatVapor,label="Vapor",c="r")
propCheckAx2.set_ylabel('Error (%)')
propCheckAx2.set_xlabel('Temperature (K)')
propCheckAx2.tick_params(direction='in')
propCheckAx2.set_title('Comparing Calculated Saturation Pressure Against Tabular Data for Propane')
propCheckAx2.legend(fontsize='small')
propCheckAx2.set_yscale('log')
ianhbell commented 1 year ago

I would really, really recommend to not implement the EOS yourself. If you want the highest speed with this EOS, look no further than teqp: https://teqp.readthedocs.io/en/latest/ . It is the fastest implementation available. The implementation in CoolProp is slower, but not too much slower, should be roughly similar to REFPROP. In fact the speed of that EOS was tested in the paper about teqp: https://pubs.acs.org/doi/10.1021/acs.iecr.2c00237

For phase equilibria you can get a massive speed boost with the superancillary EOS I developed for that EOS: https://link.springer.com/article/10.1007/s10765-021-02824-x . Available in C++ so far. Work is underway to develop the same for ALL multiparameter EOS.

nist-aharvey commented 1 year ago

To answer another part of your question, the same EOS, expressed as Eq. (32), with the same parameters, is valid for vapor and liquid.

I have a guess for where your differences are coming from. You say you used Eq. (6) to get saturated liquid densities. However, Eqs. (6) and (7) are ancillary equations, intended to give pretty good density values and to be starting points for iterative solution. The points in Table A1 were generated by rigorous solution of the EOS using the Maxwell criteria. Because it is an approximation to the rigorous solution, a density from Eq. (6) might be off by 0.001% or something from the rigorous value (maybe not enough to show up in the saturated liquid densities in Table A1 at the number of digits printed). For the vapor such an error would make negligible difference in the pressure, but pressure in the liquid is much more sensitive to density. Just making up numbers, let's say that 0.001% in density changes the pressure by 1 Pa. That will be negligible if the saturation pressure is 5 MPa, but will look like a relatively big error at lower T where the saturation pressure might be 100 Pa.

So your p(T,rho) calculations may well be correct -- of course you can check that against REFPROP for some points that are clearly in the liquid phase. It may just be a problem of using a rho that does not quite meet the saturation condition.

As Ian says, there are ways to get these properties pretty fast without going the dangerous route of programming things yourself, especially for saturation calculations if you can use his "superancillary" equations that basically replace Eqs. (6) and (7) but give values essentially matching (in double precision) what the rigorous Maxwell calculation would give.

EricLemmon commented 1 year ago

In case this is helpful, here is the quick and dirty code I wrote with the use of only the proof as my guide in order to test all the equations in the proof before we published the propane document. It is very ineloquent and of little use, but at least we knew there were no mistakes in the document.

DECLARE SUB setup ()
DEFDBL A-Z
COMMON SHARED tc, dc, pc, t0, p0, h00, s00, r
COMMON SHARED n(), t(), d(), l(), eta(), be(), g(), e()
DIM n(20), t(20), d(20), l(20), eta(20), be(20), g(20), e(20)
CALL setup
t = 369.89
d = 5
t = 300
d = 12
tau = tc / t
del = d / dc

a00 = r * t * a0(tau, del, 0)
arr = r * t * ar(tau, del, 0, 0)
a = a00 + arr
a01 = ar(tau, del, 0, 1)
a02 = ar(tau, del, 0, 2)
a03 = ar(tau, del, 0, 3)
a10 = ar(tau, del, 1, 0)
a20 = ar(tau, del, 2, 0)
a11 = ar(tau, del, 1, 1)
ai1 = a0(tau, del, 1)
ai2 = a0(tau, del, 2)
PRINT r * cp0(t), a00, a
p = d * r * t * (1 + a01)
z = p / d / r / t
u = r * t * (ai1 + a10)
h = r * t * ((ai1 + a10) + a01 + 1)
s = r * (ai1 + a10) - a / t
g = r * t * (1 + a01) + a
dpdd = r * t * (1 + 2 * a01 + a02)
dpdd2 = r * t / d * (2 * a01 + 4 * a02 + a03)
dpdt = r * d * (1 + a01 - a11)
cv = -r * (ai2 + a20)
cp = cv + r * (1 + a01 - a11) ^ 2 / (1 + 2 * a01 + a02)
w = 1 + 2 * a01 + a02 - (1 + a01 - a11) ^ 2 / (ai2 + a20)
m = 44.09562
w = SQR(w * 1000 / m * r * t)

FUNCTION a0 (tau, del, itau)
DIM a(6), b(6)
  a(1) = -4.970583
  a(2) = 4.29352
  a(3) = 3.043
  b(3) = 1.062478
  a(4) = 5.874
  b(4) = 3.344237
  a(5) = 9.337
  b(5) = 5.363757
  a(6) = 7.922
  b(6) = 11.762957

  IF itau = 0 THEN
    a = LOG(del) + 3 * LOG(tau) + a(1) + a(2) * tau
    FOR i = 3 TO 6
      a = a + a(i) * LOG(1 - EXP(-b(i) * tau))
    NEXT
  ELSEIF itau = 1 THEN
    a = 3 + a(2) * tau
    FOR i = 3 TO 6
      a = a + tau * a(i) * b(i) / (EXP(b(i) * tau) - 1)
    NEXT
  ELSEIF itau = 2 THEN
    a = -3
    FOR i = 3 TO 6
      a = a - tau ^ 2 * a(i) * b(i) ^ 2 * EXP(b(i) * tau) / (EXP(b(i) * tau) - 1) ^ 2
    NEXT
  END IF
  a0 = a
END FUNCTION

FUNCTION ar (tau, del, itau, idel)
  a = 0
  n18 = 18
  IF idel = 0 AND itau = 0 THEN
    FOR k = 1 TO 5
      a = a + n(k) * del ^ d(k) * tau ^ t(k)
    NEXT
    FOR k = 6 TO 11
      a = a + n(k) * del ^ d(k) * tau ^ t(k) * EXP(-del ^ l(k))
    NEXT
    FOR k = 12 TO n18
      a = a + n(k) * del ^ d(k) * tau ^ t(k) * EXP(-eta(k) * (del - e(k)) ^ 2 - be(k) * (tau - g(k)) ^ 2)
    NEXT
  ELSEIF idel = 1 AND itau = 0 THEN
    FOR k = 1 TO 5
      x = d(k)
      a = a + n(k) * del ^ d(k) * tau ^ t(k) * x
    NEXT
    FOR k = 6 TO 11
      x = d(k) - l(k) * del ^ l(k)
      a = a + n(k) * del ^ d(k) * tau ^ t(k) * EXP(-del ^ l(k)) * x
    NEXT
    FOR k = 12 TO n18
      x = d(k) - 2 * eta(k) * del * (del - e(k))
      a = a + n(k) * del ^ d(k) * tau ^ t(k) * EXP(-eta(k) * (del - e(k)) ^ 2 - be(k) * (tau - g(k)) ^ 2) * x
    NEXT
  ELSEIF idel = 2 AND itau = 0 THEN
    FOR k = 1 TO 5
      x = d(k) * (d(k) - 1)
      a = a + n(k) * del ^ d(k) * tau ^ t(k) * x
    NEXT
    FOR k = 6 TO 11
      x = (d(k) - l(k) * del ^ l(k)) * (d(k) - 1 - l(k) * del ^ l(k)) - l(k) ^ 2 * del ^ l(k)
      a = a + n(k) * del ^ d(k) * tau ^ t(k) * EXP(-del ^ l(k)) * x
    NEXT
    FOR k = 12 TO n18
      'x = -2 * eta(k) * del ^ 2 + 4 * eta(k) ^ 2 * del ^ 2 * (del - e(k)) ^ 2 - 4 * d(k) * eta(k) * del * (del - e(k)) + d(k) * (d(k) - 1)
      x = (d(k) - 2 * eta(k) * del * (del - e(k))) ^ 2 - d(k) - 2 * eta(k) * del ^ 2
      a = a + n(k) * del ^ d(k) * tau ^ t(k) * EXP(-eta(k) * (del - e(k)) ^ 2 - be(k) * (tau - g(k)) ^ 2) * x
    NEXT
  ELSEIF idel = 3 AND itau = 0 THEN
    FOR k = 1 TO 5
      x = d(k) * (d(k) - 1) * (d(k) - 2)
      a = a + n(k) * del ^ d(k) * tau ^ t(k) * x
    NEXT
    FOR k = 6 TO 11
      x = d(k) * (d(k) - 1) * (d(k) - 2) + del ^ l(k) * (-2 * l(k) + 6 * d(k) * l(k) - 3 * d(k) ^ 2 * l(k) - 3 * d(k) * l(k) ^ 2 + 3 * l(k) ^ 2 - l(k) ^ 3) + del ^ (2 * l(k)) * (3 * d(k) * l(k) ^ 2 - 3 * l(k) ^ 2 + 3 * l(k) ^ 3) - l(k) ^ 3 * del ^ ( _
3 * l(k))
      a = a + n(k) * del ^ d(k) * tau ^ t(k) * EXP(-del ^ l(k)) * x
    NEXT
    FOR k = 12 TO n18
      'd = d(k)
      'de = del - e(k)
      'n = eta(k)
      'x = d ^ 3 - 3 * d ^ 2 - 6 * d ^ 2 * n * del * de + 2 * d + 6 * d * n * del * de
      'x = x - 6 * d * n * del ^ 2 + 12 * d * n ^ 2 * del ^ 2 * de ^ 2
      'x = x + 12 * n ^ 2 * del ^ 3 * de - 8 * n ^ 3 * del ^ 3 * de ^ 3
      x = (d(k) - 2 * eta(k) * del * (del - e(k))) ^ 3
      x = x - 3 * d(k) ^ 2 + 2 * d(k) - 6 * d(k) * eta(k) * del ^ 2 + 6 * eta(k) * del * (del - e(k)) * (d(k) + 2 * eta(k) * del ^ 2)
      a = a + n(k) * del ^ d(k) * tau ^ t(k) * EXP(-eta(k) * (del - e(k)) ^ 2 - be(k) * (tau - g(k)) ^ 2) * x
    NEXT
  ELSEIF idel = 0 AND itau = 1 THEN
    FOR k = 1 TO 5
      a = a + t(k) * n(k) * del ^ d(k) * tau ^ t(k)
    NEXT
    FOR k = 6 TO 11
      a = a + t(k) * n(k) * del ^ d(k) * tau ^ t(k) * EXP(-del ^ l(k))
    NEXT
    FOR k = 12 TO n18
      x = t(k) - 2 * tau * be(k) * (tau - g(k))
      a = a + n(k) * del ^ d(k) * tau ^ t(k) * EXP(-eta(k) * (del - e(k)) ^ 2 - be(k) * (tau - g(k)) ^ 2) * x
    NEXT
  ELSEIF idel = 0 AND itau = 2 THEN
    FOR k = 1 TO 5
      a = a + t(k) * (t(k) - 1) * n(k) * del ^ d(k) * tau ^ t(k)
    NEXT
    FOR k = 6 TO 11
      a = a + t(k) * (t(k) - 1) * n(k) * del ^ d(k) * tau ^ t(k) * EXP(-del ^ l(k))
    NEXT
    FOR k = 12 TO n18
      x = -2 * be(k) * tau ^ 2 + 4 * be(k) ^ 2 * tau ^ 2 * (tau - g(k)) ^ 2 - 4 * t(k) * be(k) * tau * (tau - g(k)) + t(k) * (t(k) - 1)
      a = a + n(k) * del ^ d(k) * tau ^ t(k) * EXP(-eta(k) * (del - e(k)) ^ 2 - be(k) * (tau - g(k)) ^ 2) * x
    NEXT
  ELSEIF idel = 1 AND itau = 1 THEN
    FOR k = 1 TO 5
      a = a + d(k) * t(k) * n(k) * del ^ d(k) * tau ^ t(k)
    NEXT
    FOR k = 6 TO 11
      a = a + t(k) * n(k) * del ^ d(k) * tau ^ t(k) * EXP(-del ^ l(k)) * (d(k) - l(k) * del ^ l(k))
    NEXT
    FOR k = 12 TO n18
      x = (d(k) - 2 * del * eta(k) * (del - e(k))) * (t(k) - 2 * tau * be(k) * (tau - g(k)))
      a = a + n(k) * del ^ d(k) * tau ^ t(k) * EXP(-eta(k) * (del - e(k)) ^ 2 - be(k) * (tau - g(k)) ^ 2) * x
    NEXT
  END IF
  ar = a
END FUNCTION

FUNCTION cp0 (t)
DIM v(4), u(4)
v(1) = 3.043
v(2) = 5.874
v(3) = 9.337
v(4) = 7.922
u(1) = 393 / t
u(2) = 1237 / t
u(3) = 1984 / t
u(4) = 4351 / t
c = 4
FOR i = 1 TO 4
  c = c + v(i) * u(i) ^ 2 * EXP(u(i)) / (EXP(u(i)) - 1) ^ 2
NEXT
cp0 = c

END FUNCTION

FUNCTION dl (t)
  n1 = 1.82205
  n2 = .65802
  n3 = .21109
  n4 = .083973
  th = 1 - t / tc
  d = 1 + n1 * th ^ .345 + n2 * th ^ .74 + n3 * th ^ 2.6 + n4 * th ^ 7.2
  dl = dc * d
END FUNCTION

FUNCTION dv (t)
n1 = -2.4887
n2 = -5.1069
n3 = -12.174
n4 = -30.495
n5 = -52.192
n6 = -134.89
th = 1 - t / tc
d = n1 * th ^ .3785 + n2 * th ^ 1.07 + n3 * th ^ 2.7 + n4 * th ^ 5.5 + n5 * th ^ 10 + n6 * th ^ 20
dv = dc * EXP(d)
END FUNCTION

FUNCTION pv (t)
n1 = -6.7722
n2 = 1.6938
n3 = -1.3341
n4 = -3.1876
n5 = .94937
th = 1 - t / tc
p = n1 * th + n2 * th ^ 1.5 + n3 * th ^ 2.2 + n4 * th ^ 4.8 + n5 * th ^ 6.2
p = tc / t * p
pv = pc * EXP(p)
END FUNCTION

SUB setup
  tc = 369.89
  pc = 4.2512
  dc = 5
  t0 = 273.15
  p0 = .001
  h00 = 26148.48
  s00 = 157.9105
  r = 8.314472

  n(1) = .042910051
  n(2) = 1.7313671
  n(3) = -2.4516524
  n(4) = .34157466
  n(5) = -.46047898
  n(6) = -.66847295
  n(7) = .20889705
  n(8) = .19421381
  n(9) = -.22917851
  n(10) = -.60405866
  n(11) = .066680654
  n(12) = .017534618
  n(13) = .33874242
  n(14) = .22228777
  n(15) = -.23219062
  n(16) = -.09220694
  n(17) = -.47575718
  n(18) = -.017486824

  t(1) = 1
  t(2) = .33
  t(3) = .8
  t(4) = .43
  t(5) = .9
  t(6) = 2.46
  t(7) = 2.09
  t(8) = .88
  t(9) = 1.09
  t(10) = 3.25
  t(11) = 4.62
  t(12) = .76
  t(13) = 2.5
  t(14) = 2.75
  t(15) = 3.05
  t(16) = 2.55
  t(17) = 8.4
  t(18) = 6.75

  d(1) = 4
  d(2) = 1
  d(3) = 1
  d(4) = 2
  d(5) = 2
  d(6) = 1
  d(7) = 3
  d(8) = 6
  d(9) = 6
  d(10) = 2
  d(11) = 3
  d(12) = 1
  d(13) = 1
  d(14) = 1
  d(15) = 2
  d(16) = 2
  d(17) = 4
  d(18) = 1

  l(6) = 1
  l(7) = 1
  l(8) = 1
  l(9) = 1
  l(10) = 2
  l(11) = 2

  eta(12) = .963
  eta(13) = 1.977
  eta(14) = 1.917
  eta(15) = 2.307
  eta(16) = 2.546
  eta(17) = 3.28
  eta(18) = 14.6

  be(12) = 2.33
  be(13) = 3.47
  be(14) = 3.15
  be(15) = 3.19
  be(16) = .92
  be(17) = 18.8
  be(18) = 547.8

  g(12) = .684
  g(13) = .829
  g(14) = 1.419
  g(15) = .817
  g(16) = 1.5
  g(17) = 1.426
  g(18) = 1.093

  e(12) = 1.283
  e(13) = .6936
  e(14) = .788
  e(15) = .473
  e(16) = .8577
  e(17) = .271
  e(18) = .948
END SUB
mattybranham commented 1 year ago

Thank you @ianhbell, @nist-aharvey , and @EricLemmon for taking the time to give thoughtful feedback. As I have shared with @EricLemmon, I am very thankful for the work that you do in developing and distributing these fundamental data and resources. Were it not for your work, my job would be orders of magnitude more difficult. @EricLemmon, thank you for sharing your code. I ported it over to Python and ran it against my code, confirming that the outputs are consistent. @nist-aharvey, thank you for your insight! You are exactly right; the crux of my problem is that high precision is needed for liquid density to return an accurate pressure, which explains why I was having issues returning accurate liquid saturation pressures using liquid densities from Table A1 in the propane paper. @ianhbell, thank you for the suggestion to use teqp, a resource with which I was not familiar. After reading through the documentation, it is not immediately clear to me how I would use it for the problem that I am hoping to solve (i.e., running REFPROP as an array calculation), likely because I do not have a strong technical programming background. I have a request out for the two papers you cited though and will hopefully be able to move forward after reviewing them.

ianhbell commented 1 year ago

@mattybranham I can send you copies of the papers, shoot me an email at ian.bell@nist.gov.

ianhbell commented 1 year ago

And I re-iterate that I strongly advise against use of your home-brewed code. You are unlikely to achieve the same speed as REFPROP or teqp, and certainly not in Python without Cython or pybind11. If you need speed, parallelism is your friend.