usnistgov / REFPROP-wrappers

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

gamma confusion called from Mathematica #416

Closed OxyCombustion closed 2 years ago

OxyCombustion commented 2 years ago

Description

I am having problems estimating gamma by using a few different approaches in Mathematica. I seem to get different answers in different ways. I am attaching my example so you can have a look at what I am doing. I am sure I am missing something but I can't figure out what it is. I tried to attach it as a notebook (.nb file) but GitHub does not accept that file type. So, I am dropping it here as a PDF. Sorry for the inconvenience of the file format

Steps to Reproduce

Example 01.pdf

Versions

REFPROP Version: [10.0 with the Mathematica wrapper Operating System and Version: Win 10 Access Method: Mathematica wrapper

ianhbell commented 2 years ago

@henningjp Can you please look into this?

henningjp commented 2 years ago

A couple of issues that I note immediately.

  1. When you set satcomp the returned x is the liquid mole fraction and y is the vapor mole fraction. Notice that your water fraction is almost 1.0 in the x array.
  2. In the two phase region, H and E are the mole weighted values between the sat liquid and sat vapor values. Cp and Cv are LARGE because they are undefined in a two phase mixture.
  3. With a mixture of more than a few components with widely varying comps and pure fluid phases at this state point, you probably want to set the option to call SATSPLN. (See help on RefProp function for setting this option as well as the mixtures tutorial.) The REFPROP GUI does this by default, and it can give more stable behavior along the saturation curves (but not always).

General MM notes:

See what results you get, particularly with calling SATSPLN.

henningjp commented 2 years ago

Note also that “CP/CV” (gamma) is a returnable variable from the RefProp function. (See REFPROP 10 DLL documentation - table in ALLPROPS).

henningjp commented 2 years ago

@OxyCombustion - See the attached (save the .txt file and rename the extension to .nb to open it in Wolfram). What is the exact DLL version number you have (10.0.X.XX) ?

@ianhbell - I've also attached the PDF. Two-phase is not my strength. Can you explain why we're not getting exactly the CP/CV values returned from the GUI? SATSPLN didn't seem to make a difference. Are we making the saturated vapor liquid calls incorrectly, or with the wrong compositions? The numbers are really close, but not exactly what the GUI returns.

GammaCalc.txt

GammaCalc.pdf

ianhbell commented 2 years ago

If the inputs are two-phase, stop, you cannot calculate cp/cv because cp is not defined. So your first check should always be q. If it is 0 <= q <= 1, there are two phases in equilibrium.

ianhbell commented 2 years ago

Couple of notes:

import os
from ctREFPROP.ctREFPROP import REFPROPFunctionLibrary
RPPREFIX = os.getenv("RPPREFIX")
RP = REFPROPFunctionLibrary(RPPREFIX)
RP.SETPATHdll(RPPREFIX)
print('REFPROP: ', RP.RPVersion())

names = "CO2;H2O;Argon;Nitrogen;Oxygen;SO2"
zbulk = [0.5, 0.4, 0.04, 0.01, 0.005, 0.045]
p001 = 101325;
r = RP.REFPROPdll(names, "TP","D;QMOLE;CP;CV;H;E;CP/CV", RP.MOLAR_BASE_SI, 0, 0, 273.15+30, p001, zbulk)
print(r)

yielding

REFPROP:  10.0.0.02
REFPROPdlloutput(
z=array('d', [0.5, 0.4, 0.04, 0.01, 0.005, 0.045, 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', [64.54535799860605, 0.6256724410483732, -9999990.0, -9999990.0, 14853.447481378897, 13283.621189559191, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0, -9999990.0]), 
hUnits='mol/m^3', 
iUCode=100, 
x=array('d', [0.0005203367609868795, 0.9983167195759102, 1.834947183754402e-06, 1.879458926346061e-07, 1.8431967541143734e-07, 0.0011607364503512234, 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]), 
y=array('d', [0.7988288932351826, 0.042038833061723874, 0.06393011822876085, 0.01598269156608044, 0.007991291730362948, 0.07122817217788935, 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]), 
x3=array('d', [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.0]), q=0.6256724410483732, 
ierr=-810, 
herr='[REFPROP warning -810] One or more of the input properties is not thermophysically defined for two-phase states and cannot be calculated.'
)
henningjp commented 2 years ago

Thanks Ian.

Yes, he really wants the saturated liquid and vapor states. These are calculated further down in the file, but I suspect with inaccurate mole fractions than what the GUI is using.

  1. Given. H/E shouldn’t be used.
  2. I checked the Temp multiplier in Mathematica first. While not usually a good idea, his multiplier does return a Temp of 30°C. I agree, though, this is bad practice as most codes will multiply the Kelvin value by 30 and be wrong.
  3. The VAP and LIQ suffixes should work as the Mathematica wrapper just passes these through to REFPROPdll(). This convention, with the bulk composition, should give values that agree more exactly with the GUI. I’ll check that out later.
OxyCombustion commented 2 years ago

As you asked, It is version 10.0. No additional digits.

Sent from my iPhone

On Nov 1, 2021, at 6:09 AM, Jeff Henning @.***> wrote:



Thanks Ian.

Yes, he really wants the saturated liquid and vapor states. These are calculated further down in the file, but I suspect with inaccurate mole fractions than what the GUI is using.

  1. Given. H/E shouldn’t be used.
  2. I checked the Temp multiplier in Mathematica first. While not usually a good idea, his multiplier does return a Temp of 30°C. I agree, though, this is bad practice as most codes will multiply the Kelvin value by 30 and be wrong.
  3. The VAP and LIQ suffixes should work as the Mathematical wrapper just passes these through to REFPROPdll(). This convention, with the bulk composition, should give values that agree more exactly with the GUI. I’ll check that out later.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://na01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fusnistgov%2FREFPROP-wrappers%2Fissues%2F416%23issuecomment-956219572&data=04%7C01%7C%7C60dc8aa4551540b80f0808d99d38e095%7C84df9e7fe9f640afb435aaaaaaaaaaaa%7C1%7C0%7C637713689876225523%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=L%2FZChdE%2FcR4IXEjFeR5RDPU%2BC%2FVT8eLT3ryRHup7NM4%3D&reserved=0, or unsubscribehttps://na01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FAWDHHX5PJIIGT7FHTRMYBTTUJ2GRTANCNFSM5HC2KN4Q&data=04%7C01%7C%7C60dc8aa4551540b80f0808d99d38e095%7C84df9e7fe9f640afb435aaaaaaaaaaaa%7C1%7C0%7C637713689876235479%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=TJmlNxH%2FecIjGdyC%2BjlCC%2BMbTlWPT55nDJVOEBZK7SQ%3D&reserved=0. Triage notifications on the go with GitHub Mobile for iOShttps://na01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fapps.apple.com%2Fapp%2Fapple-store%2Fid1477376905%3Fct%3Dnotification-email%26mt%3D8%26pt%3D524675&data=04%7C01%7C%7C60dc8aa4551540b80f0808d99d38e095%7C84df9e7fe9f640afb435aaaaaaaaaaaa%7C1%7C0%7C637713689876235479%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=I%2FIZ53f6wzrkdU%2FC26xIbLHBp1okXvxU9U56S65Zl9A%3D&reserved=0 or Androidhttps://na01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fplay.google.com%2Fstore%2Fapps%2Fdetails%3Fid%3Dcom.github.android%26referrer%3Dutm_campaign%253Dnotification-email%2526utm_medium%253Demail%2526utm_source%253Dgithub&data=04%7C01%7C%7C60dc8aa4551540b80f0808d99d38e095%7C84df9e7fe9f640afb435aaaaaaaaaaaa%7C1%7C0%7C637713689876245429%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=ffWWgVgZ5TFHrbr8ia2fdMpQ2HuFOpSilxWjBiLuZlw%3D&reserved=0.

ianhbell commented 2 years ago

Thanks for looking into this @henningjp

@OxyCombustion Can you try to get the value as described with the LIQ and VAP suffix?

henningjp commented 2 years ago

@OxyCombustion - Mystery Solved!

Using the GUI, I get exactly the values in the notebook I sent you (and in the PDF). In the REFPROP GUI, you have to specify if you are entering mass or mole fractions for the composition. If I enter the bulk composition numbers as mole fractions, I get the saturated vapor/liquid gamma values calculated through the Mathematica wrapper; gamma = 1.01646 (liquid) and 1.30463 (vapor). If I enter your composition numbers as mass fractions, I get the saturated vapor/liquid gamma values that you reported getting from the GUI; 1.01634 (liquid) and 1.30737 (vapor).

The RefProp function in RefpropLink, by default, takes the bulk composition as mole fractions. You can add the option iMass -> 1 to specify mass-based composition in z-bulk (default is iMass -> 0 for mole fractions).

Additionally, with TP inputs, as @ianhbell suggested, you can request output properties of CPLIQ, CPVAP, CVLIQ, CVVAP, CP/CVLIQ, and CP/CVVAP to get the values along the saturation curves for the saturated liquid and saturated vapor phases; or all of them at once.

Please close this issue if this is satisfactory.

OxyCombustion commented 2 years ago

You guys are incredible. I will go over all the notes and close it out later today.

Thank you for ferreting this out.

Tom

OxyCombustion commented 2 years ago

Great explanation. Outstanding effort. Thank you both.

ianhbell commented 2 years ago

Glad to help. Thermodynamics can be confusing, so don't hesitate to ask.