usnistgov / REFPROP-wrappers

Wrappers around NIST REFPROP for languages such as Python, MATLAB, etc.
195 stars 127 forks source link

Error ABFL1dll #554

Closed ArezooHa closed 6 months ago

ArezooHa commented 1 year ago

If you have an issue to report, please continue and fill out the applicable sections below. The details provided will help to resolve this issue as quickly as possible.

Prerequisites

NOTE: Be aware that all issues are publicly accessible and viewable. Thus please do not post any code or other content that is protected intellectual property or under copyright.

Before posting an issue, please:

You can also use markdown to format your issue: GitHub guide on Markdown. If you have code snippets, please use a code block to ensure that the formatting remains legible in the web interface. For instance, surround your code in triple backticks:

print('Hello world!')

Description

I am plotting the Joule-Thomson cooling curve for natural gas samples with different compositions using Python wrapper of REFPROP (version 10). I intend to start at maximum cylinder pressure and a known temperature, and calculate temperature points as pressure is reduced. In order to use JTCRVdll subroutine, I would need to calculate the gas composition density at each pressure point, for which I am trying to use ABFL1dll subroutine. However, an error is returned:

RP.ABFL1dll(T_initial, P_initial, z, kph, 'TP', Dmin, Dmax, tt, pp, D, ierr, herr) TypeError: ABFL1dll() takes 8 positional arguments but 13 were given

The syntax I am using is based on REFPROP Documentation Release 10.0, with 12 positional arguments. I was wondering if there's any way I can resolve this error or if another approach can be used for plotting joule-Thomson cooling curve?

Steps to Reproduce

Here's the line for density calculations. Z is a mixture of mainly methane and some larger hydrocarbons.

RP.ABFL1dll(T_initial, P_initial, z, kph, 'TP', Dmin, Dmax, tt, pp, D, ierr, herr)

Expected behavior: I expect to be able to calculate density using ABFL1dll, to be fed into JTCRVdll subroutine.

Actual behavior: An error is returned as: TypeError: ABFL1dll() takes 8 positional arguments but 13 were given

Versions

REFPROP Version: 10 Operating System and Version:Windows 10
Access Method: Python

Additional Information

If possible, please post examples and/or screenshots of the issue.

ianhbell commented 1 year ago

You can get all of these things from the REFPROPdll function. Can you please provide a working example?

ArezooHa commented 1 year ago

Unfortunately, I still don't have a working example for JT cooling curve. Being able to use REFPROPdll would be great, as I would like to plot both phase boundary and cooling curve on the same plot. However, I am not getting expected results from an isenthaplic caclulation of JT temperatures. I have attached a snippet of the code as well as the resulted plots (in both logarithmic and non-logarithmic P scale).

Dew-point T at max P calculation

tdewr = RP.REFPROPdll(" ", "PQ", "T", MOLAR_BASE_SI, 0, 0, p_Pa, 1, [MoleF1,MoleF2,MoleF3,MoleF4,MoleF5,MoleF6,MoleF7,MoleF8,MoleF9,MoleF10,MoleF11,MoleF12,MoleF13,MoleF14,MoleF15,MoleF16,MoleF17,MoleF18,MoleF19,MoleF20]).Output[0]

Define the initial temperature (in K) for Joule-Thomson cooling curve

min_temperature = tdewr + 11.1

Calculate the enthalpy at the initial conditions

enthalpy = RP.REFPROPdll("", "PT", "H", MOLAR_BASE_SI, 0, 0, initial_pressure, initial_temperature, [MoleF1,MoleF2,MoleF3,MoleF4,MoleF5,MoleF6,MoleF7,MoleF8,MoleF9,MoleF10,MoleF11,MoleF12,MoleF13,MoleF14,MoleF15,MoleF16,MoleF17,MoleF18,MoleF19,MoleF20]).Output[0]

Programatically extract the nodes of the isopleth

Nnodes = RP.REFPROPdll("","SPLINENODES","",0,0,0,0,0,[MoleF1,MoleF2,MoleF3,MoleF4,MoleF5,MoleF6,MoleF7, MoleF8,MoleF9,MoleF10,MoleF11,MoleF12,MoleF13,MoleF14,MoleF15,MoleF16,MoleF17,MoleF18,MoleF19,MoleF20]).iUCode data = [] for i in range(1,Nnodes+1): rho_vap,T,p,rho_eq,z0,z1 = RP.REFPROPdll("","SPLINENODES","",0,0,i,0,0,[MoleF1,MoleF2,MoleF3,MoleF4,MoleF5,MoleF6,MoleF7,MoleF8,MoleF9,MoleF10,MoleF11,MoleF12,MoleF13,MoleF14,MoleF15,MoleF16,MoleF17,MoleF18,MoleF19,MoleF20]).Output[0:6] data.append(dict(rho_vap = rho_vap,T=T,p=p,z0=z0,z1=z1))

Convert the data points into a pandas DataFrame

df = pandas.DataFrame(data)

Calculate JT temperatures in an isenthaplic process

def calculate_JTtemperature(p): JTtemperature = RP.REFPROPdll("", "PH", "T", MOLAR_BASE_SI, 0, 0, p, enthalpy, [MoleF1,MoleF2,MoleF3,MoleF4,MoleF5,MoleF6,MoleF7,MoleF8,MoleF9,MoleF10,MoleF11,MoleF12,MoleF13,MoleF14,MoleF15,MoleF16,MoleF17,MoleF18,MoleF19,MoleF20]).Output[0] return JTtemperature

Apply the calculate_JTtemperature function to the 'p' column and store results in a new column 'calculated_JTtemperature'

df['calculated_JTtemperature'] = df['p'].apply(calculate_JTtemperature)

Plot the data

Create a plot with 'T' on the x-axis

ax1 = df.plot(x='T', y='p', label='Phase Boundary Data', color='black', logy=True) #Set logy=True for logarithmic y-axis

Create a second plot with 'calculated_JTtemperature' on the x-axis

ax2 = df.plot(x='calculated_JTtemperature', y='p', ax=ax1, label='Joule-Thomson Cooling Curve', color='blue', logy=True)

Set the x-axis range to span the full range of T

ax1.set_xlim(df['T'].min(), df['T'].max()+10)

image image

ianhbell commented 1 year ago

Please provide a runnable example

ianhbell commented 1 year ago

I cleaned up your example to follow python guidelines a bit better:

z = [MoleF1,MoleF2,MoleF3,MoleF4,MoleF5,MoleF6,MoleF7,MoleF8,MoleF9,MoleF10,MoleF11,MoleF12,MoleF13,MoleF14,MoleF15,MoleF16,MoleF17,MoleF18,MoleF19,MoleF20]

# Dew-point T at max P calculation
tdewr = RP.REFPROPdll(" ", "PQ", "T", MOLAR_BASE_SI, 0, 0, p_Pa, 1, z).Output[0]
# Define the initial temperature (in K) for Joule-Thomson cooling curve
min_temperature = tdewr + 11.1
# Calculate the enthalpy at the initial conditions
enthalpy = RP.REFPROPdll("", "PT", "H", MOLAR_BASE_SI, 0, 0, initial_pressure, initial_temperature, z).Output[0]
# Programatically extract the nodes of the isopleth
Nnodes = RP.REFPROPdll("","SPLINENODES","",0,0,0,0,0,z).iUCode
data = []
for i in range(1,Nnodes+1):
    rho_vap,T,p,rho_eq,z0,z1 = RP.REFPROPdll("","SPLINENODES","",0,0,i,0,0,z).Output[0:6]
    data.append(dict(rho_vap = rho_vap,T=T,p=p,z0=z0,z1=z1))
# Convert the data points into a pandas DataFrame
df = pandas.DataFrame(data)
# Calculate JT temperatures in an isenthaplic process
def calculate_JTtemperature(p):
    JTtemperature = RP.REFPROPdll("", "PH", "T", MOLAR_BASE_SI, 0, 0, p, enthalpy, z).Output[0]
    return JTtemperature
# Apply the calculate_JTtemperature function to the 'p' column and store results in a new column 'calculated_JTtemperature'
df['calculated_JTtemperature'] = df['p'].apply(calculate_JTtemperature)
# Plot the data
# Create a plot with 'T' on the x-axis
ax1 = df.plot(x='T', y='p', label='Phase Boundary Data', color='black', logy=True) #Set logy=True for logarithmic y-axis
# Create a second plot with 'calculated_JTtemperature' on the x-axis
ax2 = df.plot(x='calculated_JTtemperature', y='p', ax=ax1, label='Joule-Thomson Cooling Curve', color='blue', logy=True)
# Set the x-axis range to span the full range of T
ax1.set_xlim(df['T'].min(), df['T'].max()+10)
ArezooHa commented 1 year ago

Thank you for your feedback. Below is a runnable code:

# Import the main class from the Python library
from ctREFPROP.ctREFPROP import REFPROPFunctionLibrary
# Imports from the standard library
import glob
# Imports from conda-installable packages
import pandas
# Import numpy
import numpy as np
# Import matplotlib for plotting
import matplotlib.pyplot as plt
# explicitly state which path we want to use.
import os
RP = REFPROPFunctionLibrary(os.environ['RPPREFIX'])
RP.SETPATHdll(os.environ['RPPREFIX'])
#Get the unit system
MOLAR_BASE_SI = RP.GETENUMdll(0, "MOLAR BASE SI").iEnum
#Specify calculation model
RP.FLAGSdll("GERG",1)
#Obtain the isopleth data
#Set up the fluids
RP.SETFLUIDSdll('Hydrogen * Helium * Water * Nitrogen * Oxygen * Hydrogen sulfide * Carbon dioxide * Methane * Ethane * Propane * Iso-butane * N-butane * 2,3-Dimethyl butane * Iso-pentane * N-pentane * N-hexane * N-heptane * N-octane * N-nonane * N-decane')
#Import win32com
import win32com.client
#Get maximum cylinder pressure in Pa
p_Pa = 750 *0.0068947573
#Get mole fraction values
MoleF1 = 0.000000 
MoleF2 = 0.000000 
MoleF3 = 0.000000 
MoleF4 = 0.027680 
MoleF5 = 0.000000 
MoleF6 = 0.000000 
MoleF7 = 0.004917 
MoleF8 = 0.899237 
MoleF9 = 0.049770 
MoleF10 = 0.010040 
MoleF11 = 0.003005 
MoleF12 = 0.003017 
MoleF13 = 0.000000 
MoleF14 = 0.001013 
MoleF15 = 0.001009 
MoleF16 = 0.000312 
MoleF17 = 0.000000 
MoleF18 = 0.000000 
MoleF19 = 0.000000 
MoleF20 = 0.000000 
#Construct the saturation spline
z = [MoleF1,MoleF2,MoleF3,MoleF4,MoleF5,MoleF6,MoleF7,MoleF8,MoleF9,MoleF10,MoleF11,MoleF12,MoleF13,MoleF14,MoleF15,MoleF16,MoleF17,MoleF18,MoleF19,MoleF20]
RP.SATSPLNdll(z)
#Dew-point T at max P calculation
tdewr = RP.REFPROPdll(" ", "PQ", "T", MOLAR_BASE_SI, 0, 0, p_Pa, 1, z).Output[0]
#Define the initial temperature (in K) for Joule-Thomson cooling curve
min_temperature = tdewr + 11.1
#Set the initial for Joule-Thomson cooling curve using maximum cylinder pressure and minimum use temperature
initial_temperature = min_temperature
initial_pressure = p_Pa
#Calculate the enthalpy at the initial conditions
enthalpy = RP.REFPROPdll("", "PT", "H", MOLAR_BASE_SI, 0, 0, initial_pressure, initial_temperature, z).Output[0]
#Programatically extract the nodes of the isopleth
Nnodes = RP.REFPROPdll("","SPLINENODES","",0,0,0,0,0,z).iUCode
data = []
for i in range(1,Nnodes+1):
    rho_vap,T,p,rho_eq,z0,z1 = RP.REFPROPdll("","SPLINENODES","",0,0,i,0,0,z).Output[0:6]
    data.append(dict(rho_vap = rho_vap,T=T,p=p,z0=z0,z1=z1))
#Convert the data points into a pandas DataFrame
df = pandas.DataFrame(data)
# Calculate JT temperatures in an isenthaplic process
def calculate_JTtemperature(p):
    JTtemperature = RP.REFPROPdll("", "PH", "T", MOLAR_BASE_SI, 0, 0, p, enthalpy, z).Output[0]
    return JTtemperature
#Apply the calculate_JTtemperature function to the 'p' column and store results in a new column 'calculated_JTtemperature'
df['calculated_JTtemperature'] = df['p'].apply(calculate_JTtemperature)
df['p'] = df['p'] / 1000
#Plot the data
#Create a plot with 'T' on the x-axis
ax1 = df.plot(x='T', y='p', label='Phase Boundary Data', color='black', logy=True) #Set logy=True for logarithmic y-axis
#Create a second plot with 'calculated_JTtemperature' on the x-axis
ax2 = df.plot(x='calculated_JTtemperature', y='p', ax=ax1, label='Joule-Thomson Cooling Curve', color='blue', logy=True)
ax1.set_xlabel('Temperature (K)')
ax1.set_ylabel('Pressure (MPa)')
#Show legends
ax1.legend(loc='lower left')
ax2.legend(loc='lower left')
#Set the x-axis range to span the full range of T
ax1.set_xlim(df['T'].min(), df['T'].max()+10)
plt.show()
ianhbell commented 1 year ago

Have you tried to enable the saturation splines in your call to REFPROPdll with the J-T calculations? Otherwise your code looks fine, and it could just be a flash calculation failure with REFPROP. You might have to do the J-T calculations yourself. It is not very convenient, but might be your only option in this case.

ianhbell commented 1 year ago

I changed this line:

JTtemperature = RP.REFPROPdll("", "PH", "T", MOLAR_BASE_SI, 1, 0, p, enthalpy, z).Output[0]

and success:

image
ArezooHa commented 1 year ago

Thank you for the new edit. I only realized that the curve seems to deviate from what is anticipated. Checking on enthalpy calculations, the line below (enthalpy calculation at minimum use temperature (derived from dew point calculations) and the gas cylinder's maximum pressure):

enthalpy = RP.REFPROPdll("", "PT", "H", MOLAR_BASE_SI, 0, 0, initial_pressure, initial_temperature, z).Output[0]

returns a value of -2701. Enabling the saturation splines in this call to REFPROP changes the returned enthalpy value to -2606. I am not able to explain the negative value. And my apologies, this line needs to be changed too:

p_Pa = 750 *6894.7573

ianhbell commented 1 year ago

There is in principle nothing wrong with negative enthalpy. It all comes down to the selected reference state.

ArezooHa commented 1 year ago

You are right. I'm sorry I was looking at the wrong Y-axis scaling when checking the cooling curve.