search

usnistgov / REFPROP-issues

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

THERM function returns the wrong value for Enthalpy #552

Closed KayvanA closed 1 year ago

KayvanA commented 2 years ago

INTRODUCTION/EXPLANATION:

I am calling THERM function to calculate the enthalpy of saturated liquid for R410A at 20C=293.15K.

I use 'Molar Base SI' for the calculations. I get iUnit=20 by calling: GETENUM(0,'MOLAR BASE SI',iUnit,ierr,here) I have also tried using iUnit=0 as Refprop documentation shows iUnit=0 as the Default for Refprop units and it is supposed to be for 'MOLAR BASE SI'. So I have tried my code both with iUnit=0 and iUnit=20.

I call SETMIXTURE subroutine for R410A, and the composition is returned correct: z = [0.697614699375863 0.302385300624138]

Before I call THERM, I use SATT function to calculate the saturated liquid density. SATT produces the correct saturation pressure (Psat=1447.60377810720 kPa), and also correct saturated liquid density (Dl=14.9212417496128 mol/m3). The pressure and density match with the values received from REFPROP1 function, and also the REFPROP call in Excel. Also, WMOL returns the correct molar density of wmm=72.5854142406598 g/mol.

Later, I use this saturated liquid density (Dl=14.9212417496128 mol/m3) in my call to STN subroutine to calculate the surface tension (sigma=5.920927775333686E-003 N/m). The same saturated liquid density is used in my call to TRNPRP and it produces the correct values for viscosity (eta=125.868581016964 uPa.s) and thermal conductivity (tcx=9.195081001175499E-002 W/(m.K) ). When I compare the values of sigma, eta, or tcx with the values from REFPROP1 call in Fortran, or REFPROP call in Excel, they all match.

I use this saturated liquid density (Dl=14.9212417496128 mol/m3) in my call to THERM. It produces the correct value for specific heat capacity (Molar_Cp=120.337313490060 J/(mol.K) ). I can convert this Cp to Mass_Cp (J/(kg.K)) using the molecular weight (wmm) as: Mass_Cp = Molar_Cp / wmm (1000 g/kg) which gives me: Mass_Cp = 1657.87182933305 J/(kg.K) = 120.337313490060 J/(mol.K) / (72.5854142406598 g/mol) (1000 g/kg) My call to REFPROP1 in Fortran produces the same Molar Cp as THERM . My call to REFPROP in Excel produces the same Mass_Cp as calculated above. So they all match.

PROBLEM/ISSUE:

As mentioned, I use the correct saturated liquid density (Dl=14.9212417496128 mol/m3) calculated from SATT in my call to THERM. The value of saturated liquid enthalpy that I calculate via THERM or REFPROP1 in Fortran is: Molar_h= 16894.6680926860 J/mol

I convert this value to J/kg via Mass_h = Molar_h / wmm (1000 g/kg) Mass_h = 232755.688858798 J/kg = 16894.6680926860 J/mol / (72.5854142406598 g/mol) (1000 g/kg)

However, this value of enthalpy does not match with the Refprop call in Excel with 'MASS BASE SI': Mass_h from Excel= 231541.858128479 J/kg Therefore, there is a discrepancy of 1214 J/kg between that calculated via REFPROP call in Excel and that calculated via REFPROP1 and THERM.

I should also note that if I use 'MOLAR BASE SI' in Excel and call REFPROP for Enthalpy, it produces: Molar_h = 16806.561686308 J/mol This value is consistent with 'MASS BASE SI' value of enthalpy in Excel, and it satisfies the equation: Mass_h = Molar_h / wmm * (1000 g/kg)

QUESTIONS:

For R410A at 20C=293.15K,

1- Could you please let me know what value of saturated liquid enthalpy (x=0) you would get either in Fortran or Python by calling THERM subroutine?

2-Could you please let me know what value of saturated liquid enthalpy (x=0) you would get either in Fortran or Python by calling REFPROP1 subroutine?

3-Could you please let me know what value of saturated liquid enthalpy (x=0) you would get either in Fortran or Python by calling REFPROP subroutine and use 'MOLAR BASE SI' ?

4-Do you get a discrepancy between the values of Questions 1, 2, and 3?

5-Could you please let me know what value of saturated liquid enthalpy (x=0) you would get either in Fortran or Python by calling REFPROP subroutine and use 'MASS BASE SI'?

6-If you convert the enthalpy from REFPROP1 from J/mol to J/kg (from Question 2) and compare it to what you got from REFPROP (from Question 5), do you get a discrepancy? How can you explain the discrepancy?

I have also attached my Fortran file for reference.

Thank you, Kayvan

KayvanA commented 2 years ago 
  subroutine test_call()

  implicit none

  integer, parameter :: ncmax=20,ipropmax=200
  integer, parameter :: n_return_props1=7,n_return_props2=5
  double precision dummy1(ncmax),dummy2(ncmax),dummy3(ncmax)
  double precision x,y,xcalc,z(ncmax),Output(ipropmax) 
  double precision xx(ncmax), yy(ncmax)
  integer iUnit,iMass

  common /Refprop_MC_commons/ iUnit, iMass, z

  double precision P,T,x1,x2,x3,h,hf1,hf2,hg,Tf,hf,error_Tf,Pf
  double precision Dl, Dv
  double precision T_calc_liquid, T_calc_vapor 
  double precision P_calc_liquid, P_calc_vapor
  double precision P_kPa , wmm , WMOL
  double precision rho_sat_l_1, rho_sat_v_1
  double precision rho_sat_l_2, rho_sat_v_2
  double precision error_P, Answer(3)
  double precision sigma
  double precision e_per_mol,h_per_mol,s_per_mol,Cv_per_mol,Cp_per_mol,w_per_mol,hjt_per_mol
  double precision error_h, H_REFPROP1, hf_THERM
  double precision eta_TRNPRP, tcx_TRNPRP, eta_REFPROP1, tcx_REFPROP1, error_eta, error_tcx
  double precision error_sigma, sigma_REFPROP1, sigma_STN
  double precision error_Cp, Cp_REFPROP1

  integer ierr,iFlag,iOut, kph , herr_length

  character*255 hOut,herr,Refrigerant

  herr_length = 255

  Refrigerant='R410A'
  !Refrigerant='R32'

  !iMass=1 
  !call GETENUM(0,'MASS BASE SI',iUnit,ierr,herr)

  iMass = 0 
  call GETENUM(0,'MOLAR BASE SI',iUnit,ierr,herr) ! This returns iUnit = 20

  iUnit = 0 ! Refprop documentation mentions that the default units are 'MOLAR BASE SI' with iUnit=0. Therefore, I am manually setting/overwriting iUnit to 0.

  !Calculate the composition (z)
  call SETMIXTURE (Refrigerant, z, ierr)

  !Calculate the molar weight
  wmm = WMOL(z)   ! g/mol
  !wmm = 72.58_8   ! Molecular weight of R410A

  T = 20.0_8 + 273.15_8  ! Convert from Celsius to Kelvin

  kph = 1
  call SATT (T,z,kph,P_kPa,Dl,Dv,xx,yy,ierr,herr)

  ! Dl and Dv are in mol/L
  rho_sat_l_1 = wmm * Dl     !The result will be in g/L = kg/m3
  rho_sat_v_1 = wmm * Dv     !The result will be in g/L = kg/m3

  !Check SATT result against REFPROP1 result
  x = 0_8
  call REFPROP1('TQ','P',iUnit,iMass,T,x,z,P_calc_liquid,xcalc,ierr,herr)

  error_P = P_calc_liquid - P_kPa    

  !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  ! Viscosity an Thermal conductivity %
  !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  call TRNPRP (T,Dl,z,eta_TRNPRP,tcx_TRNPRP,ierr,herr)   
  call REFPROP1('TQ','VIS',iUnit,iMass,T,x,z,eta_REFPROP1,xcalc,ierr,herr)
  call REFPROP1('TQ','TCX',iUnit,iMass,T,x,z,tcx_REFPROP1,xcalc,ierr,herr)

  error_eta = eta_REFPROP1 - eta_TRNPRP
  error_tcx = tcx_REFPROP1 - tcx_TRNPRP

  !%%%%%%%%%%%%%%%%%%
  ! Surface tension %
  !%%%%%%%%%%%%%%%%%%

  call STN (T,Dl,Dv,xx,yy,sigma_STN,ierr,herr)
  call REFPROP1('TQ','STN',iUnit,iMass,T,x,z,sigma_REFPROP1,xcalc,ierr,herr)

  error_sigma = sigma_REFPROP1 - sigma_STN

  !%%%%%%%%%%%%%%%%%%
  ! Enthalpy and cp %
  !%%%%%%%%%%%%%%%%%%

  call THERM (T,Dl,z,P,e_per_mol,h_per_mol,s_per_mol,Cv_per_mol,Cp_per_mol,w_per_mol,hjt_per_mol)
  hf_THERM = h_per_mol / wmm * 1000_8    ! (J/mol)/(g/mol)*(1000g/kg) = J/kg

  call REFPROP1('TQ','H',iUnit,iMass,T,x,z,H_REFPROP1,xcalc,ierr,herr)
  error_h = H_REFPROP1 - h_per_mol
  !error_h = H_REFPROP1 - hf_THERM

  call REFPROP1('TQ','CP',iUnit,iMass,T,x,z,Cp_REFPROP1,xcalc,ierr,herr)

  error_Cp = Cp_REFPROP1 - Cp_per_mol

end subroutine
ianhbell commented 2 years ago

Some comments: make sure you are always checking ierr after every call. The default units in REFPROP are NOT MOLAR BASE SI, they are DEFAULT, which is a kilo SI variant.

I think there is a bug in the Excel wrapper. I did what you said in Python, and I get the right answer with REFPROPdll.

import os 
root = os.getenv('RPPREFIX')
import ctREFPROP.ctREFPROP as ct
RP = ct.REFPROPFunctionLibrary(root)
RP.SETPATHdll(root)
print(RP.RPVersion())

m = RP.SETMIXTUREdll('R410A')
T = 20.0 + 273.15
s = RP.SATTdll(T, m.z, 1)

print(s.Dl)
print(RP.ENTHALdll(T, s.Dl, m.z))
print(RP.THERMdll(T, s.Dl, m.z).h)
w = RP.WMOLdll(m.z)/1000
print(RP.ENTHALdll(T, s.Dl, m.z)/w)

Dlmass = s.Dl*1000*w # kg/m^3
print(RP.REFPROPdll('','TD&','H',RP.MASS_BASE_SI,0,0,T,Dlmass,m.z).Output[0])

yielding

10.0.1.20
14.921241749612772
16894.652628777152
16894.652628777152
232755.4758144972
232755.47581449724
ianhbell commented 2 years ago

Do you want to try the beta which might fix this bug? If so, please email me at ian.bell@nist.gov

KayvanA commented 2 years ago

Thank you for your reply and for writing the code for me, Ian. Before I ask for the beta version, I would like to check a few things:

Could you please check something for me:

When you call REFPROPdll, can you please use quality instead of density (T=293.15 and Q=0) and use 'MOLAR BASE SI'? When I do that in Excel, I get 16806.561686308 J/mol from REFPROP call in Excel. Could you please let me know what you get in Python when you call REFPROPdll? I ask because I noticed that you used density in your call to REFPROPdll. I wonder if the call via quality instead of density is causing the discrepancy. Also, please note that the value returned from REFPROP in Excel (16806.561686308J/mol) is different from the values returned by THERM, ENTHAL, and REFPROP1 (16894.6680926860 J/mol), by 88.1J/mol, which is significant.

Sidenote 1: MASS BASE SI and MOLAR BASE SI are consistent in Excel

Sidenote 2: REFPROP iUnit

As for iUnit value: GETENUM gives me iUnit = 20. I could not find the value of iUnit=20 in the Refprop documentation pdf file. If I do not change iUnit to 0 and keep iUnit=20, then the unit for viscosity will be wrong, it will be off by 4 decimals: REFPROP1 subroutine would give me eta=1.258685810169644E-004 instead of giving the correct value of eta = 125.868581016964 uPa.s.

Also, Refprop documentation pdf file mentions the below sentence: "In all environments other than FORTRAN, the iUnits variable should be retrieved from the GETENUM function with a call like: GETENUMdll(0,'MOLAR BASE SI',iEnum,ierr,here)" Does that mean that I should not call GETENUM in FORTRAN and manually set iUnit=0 ?

Also, here is a picture of what the Refprop documentation pdf shows for iUnit=0:

image Please note that the document mentions the unit of pressure for Molar SI as MPA. However, if I use iUnit=20 from GETENUM , I still get my result in kPa, which is what iUnit=0 (Default) shows in the above picture.

Sidenote 3: ierr is always 0.

Per your advice, I stepped through the code and checked all ierr after every call, with (iUnit=0) and (iUnit=20). It was always 0.

Thank you, Kayvan

ianhbell commented 2 years ago

Well, that small difference is caused by gas constant, I am 99.999999999% sure. Check your gas constant. I am using the beta which uses the gas constant in the fluid file by default. See https://github.com/usnistgov/REFPROP-issues/issues/239

Never, ever use hard-coded values for unit systems. That is the short answer. The molar base SI are listed here:

image

It is also in the PDF: https://refprop-docs.readthedocs.io/_/downloads/en/latest/pdf/

KayvanA commented 2 years ago

Thank you for the clarification about the iUnits and also thanks for the updated Refprop doc (23 March 2022). My Refprop 10 doc was from 2019. Also thanks for the info about the Universal Gas constant and its modified value.

Could you please call the REFPROPdll in your Python code to calculate the enthalpy with quality(x=0) instead of using density(Dl)? Please use 'MOLAR BASE SI' so we can compare notes.

When I do this call in Excel, I get 16806.561686308 J/mol from REFPROP call. Could you please let me know what you get in Python when you call REFPROPdll? I ask because I noticed that you used density in your call to REFPROPdll. I wonder if the call via quality instead of density is causing the large discrepancy. Also, please note that the value returned from REFPROP in Excel (16806.561686308J/mol) is different from the values returned by THERM, ENTHAL, and REFPROP1 (16894.6680926860 J/mol), by 88.1J/mol, which is significant.

ianhbell commented 2 years ago

The density will be the same, that is not the cause. I used TD& to skip the phase calculation.

Excel and REFPROP use different reference states for reasons I won't get into here. So you need to force one of them to match the other one by a call to SETREF, if you want them to be the same.

KayvanA commented 2 years ago

I was able to call REFPROP in Fortran once with (MOLAR BASE SI) and once with (MASS BASE SI). They are consistent and satisfy the below equation:

h_MASS_BASED= h_MOLAR_BASED / wmm 1000 , unit-wise: J/kg = (J/mol)/(g/mol)(1000g/kg)

The discrepancy between the manually calculated variable (h_MASS_BASED) and the value calculated via REFPROP with MASS BASE SI (in variable H_REFPROP_Mass) was 2.910383045673370D-011 J/kg, which is negligible.

I would agree that the discrepancy is due to the Reference state for enthalpy. To make sure there are no errors, I will calculate the enthalpy of evaporation h_fg = h_g - h_f. In this case, the reference state would have to cancel out and both Excel and Fortran should produce the same value. I will let you know how that goes in my next comment.

Thank you, Kayvan

Sidenote 1: Calling REFPROP in Fortran with (MOLAR BASE SI)

It is worth mentioning that the call to REFPROP in Fortran created the exact same value as THERM, ENTHAL( THERM1), and REFPROP1, which was reassuring: Molar_h = 16894.6680926860 J/mol In my call to REFPROP and REFPROP1 , I used quality(x=0) instead of density.

Side note 2: FORTRAN Code:

For reference, here are my calls to REFPROP in the FORTRAN code:

  ! 1 - Calling REFPROP with MOLAR BASE SI
  iMass = 0 ! The inputs variables are mole fraction based
  call GETENUM(0,'MOLAR BASE SI',iUnit,ierr,herr) ! This returns iUnit = 20
  call REFPROP (Refrigerant,'TQ','H',iUnit,iMass,iFlag,T,x,z_composition,Output,hUnits,iUCode,xx,yy,zz,xcalc,ierr,herr)
  H_REFPROP_Molar = Output(1)

  ! 2 - Calling REFPROP with MASS BASE SI
  iMass = 1 ! The inputs variables are mass fraction based
  call GETENUM(0,'MASS BASE SI',iUnit,ierr,herr) ! This returns iUnit = 21
  call REFPROP (Refrigerant,'TQ','H',iUnit,iMass,iFlag,T,x,z_composition,Output,hUnits,iUCode,xx,yy,zz,xcalc,ierr,herr)
  H_REFPROP_Mass = Output(1)

  ! 3- Error/discrepancy between molar and mass based
  H_MASS_Calc   = H_REFPROP_Molar / wmm * 1000_8    ! J/kg = (J/mol)/(g/mol)*(1000g/kg)
  error_REFPROP = H_REFPROP_Mass - H_MASS_Calc
KayvanA commented 2 years ago

I tried calculating the enthalpy of evaporation (h_fg) with REFPROP subroutine in Fortran; Molar and Mass-based calculations were consistent with the molecular weight. Therefore, I agree that the difference in REFPROP call in Excel and Fortran was in the reference state.

However, I noticed a new issue, which I would be happy to to post in a separate page, if you wish. You may want to consider this bug fix in a future version of REFPROP:

When I use 'MASS BASE SI' and I call REFPROP Subroutine, it changes the composition from [0.697614699375863 0.302385300624138] to [0.5 0.5]. It is probably changing the composition from molar basis to mass basis. This results in the next call to REFPROP subroutine not working: ierr = 813 herr = '[REFPROP error 813] The input composition must match that of the predefined mixture when the mixture name is included in the call to REFPROP. To use a different composition, leave the fluid name blank. '

To fix this, I had to set the composition again before the next REFPROP call, which is obviously time-consuming and redundant: call SETMIXTURE (Refrigerant, z_composition, ierr) ! If I turn this SETMIXTURE Off, then z_composition will be messed up and next REFPROP call won't work.

According to the Refprop documentation, the composition (z) is an INPUT and not an OUTPUT. https://refprop-docs.readthedocs.io/_/downloads/en/latest/pdf/ image

Therefore, one would expect the composition to remain constant after the call to REFPROP. I should mention that I do set iMass=1 for my call to 'MASS BASE SI', in both REFPROP calls. Probably, somewhere in the subroutine, the composition gets altered and it is overwritten on the input value, which could be easily fixed by adding a new variable in the subroutine.

Lastly, this issue happens when I use 'MASS BASE SI'. When I switch back to 'MOLAR BASE SI', the composition remains the same before and after calling REFPROP subroutine.

I will paste the relevant part of my code for reference in the next comment.

ianhbell commented 2 years ago

I think the issue with changing the composition is a bug. Please file a separate issue for that and I will get our team to investigate.

KayvanA commented 2 years ago

will do, thanks.