Using REFPROP for a binary mixture VLE (PH-flash) for species that aren't in FLD files?

[wandadars]

Can REFPROP python interface be used to perform a PH-flash calculation for a mixture of two species if those species don't have available FLD files in the refprop database? I saw that the Peng-Robinson equation of state can be used, and that only requires critical properties and acentric factors, but I'm getting lost in the documentation. Can something like this be done?

[ianhbell]

Yes, you need to represent each species with its own PR model. Here is a template for a PR model:

xxxxxxxxxxxxxxx    !short name
xxx-xx-x           !CAS number
x                  !full name
x                  !chemical formula
x                  !synonym
100.               !molecular weight [g/mol]
200.               !triple point temperature [K]
300.               !normal boiling point [K]
400.               !critical temperature [K]
5000.              !critical pressure [kPa]
6.                 !critical density [mol/L]
0.2                !acentric factor
-1.                !dipole moment [Debye] (use -1 for unknown)
NBP                !default reference state
8.0                !version number


! Add the name of the fluid you are using on the first line, and update
! the other lines with the properties of the fluid.  Below add the molar mass,
! critical point, and acentric factor. A Cp0 equation is required if you are
! using caloric properties (H, Cp, Snd, etc.)


#EOS
PRT  translated Peng-Robinson equation
?LITERATURE REFERENCES \
!end of info section
200.               !lower temperature limit [K]
700.               !upper temperature limit [K]
100000.            !upper pressure limit [kPa]
50.                !maximum density [mol/L]
CPP                !pointer to Cp0 model
100.               !molecular weight [g/mol]       |update these
0.2                !acentric factor                |
400.               !critical temperature [K]       |
5000.              !critical pressure [kPa]        |
6.                 !critical density [mol/L]       |
8.3144598          !gas constant [J/mol-K]
1                  !Number of parameters for translation
0.0                !value of translation parameter


#AUX               !auxiliary model specification
CPP  ideal gas heat capacity function
?LITERATURE REFERENCE \
?\
!end of info section
50.0               !lower temperature limit [K]
5000.0             !upper temperature limit [K]
0.0                !upper pressure limit [kPa]
0.0                !maximum density [mol/L]
1.0          8.3144598                 !reducing parameters for T, Cp0
  1  0    0  0    0  0  0              !Nterms:  polynomial, exponential, cosh, sinh
 3.5             0.0


@END
c        1         2         3         4         5         6         7         8
c2345678901234567890123456789012345678901234567890123456789012345678901234567890

[wandadars]

I'm trying to write (what I thought) was a straight-forward Python PH-query of a mixture of two components that are governed by the Peng-Robinson equation of state. Below is my script that tries to compute the pure-species enthalpy for the MMH and NTO species at different temperatures, but the same pressure. I thought that I could prescribe a mixture enthalpy that is the mole-weighted sum of the pure species enthalpies of the pure species, and then perform a VLE that would give the composition of the binary mixture in an equilibrium state with the P and enthalpy fixed.

One thing that is perplexing me is that a call to try and get the saturation temperature at P=1atm, the call seems to fail.

t_sat = rp.REFPROPdll('MMH', 'PQ', 'T', units, 0, 0, P, 0.0, [1.0]).Output[0]

The output I see is:
MMH Saturation temperature at P= 101.325 kPa: -9999990.0 Kelvin

Full script and fluid files provided below for reproducing.

import numpy as np
import os
import pandas
import matplotlib.pyplot as plt

from ctREFPROP.ctREFPROP import REFPROPFunctionLibrary

def plot_mixture_phase_boundaries(mixture_name, mixture_composition, P_in, T_in, rp):
    '''
    This method only works for mixtures with multiple components and can't be used for 
    pure species
    '''
    rp.SETFLUIDSdll(mixture_name)
    # Construct the saturation spline
    rp.SATSPLNdll(mixture_composition)
    # Now we programatically extract the nodes of the isopleth...
    Nnodes = rp.REFPROPdll("","SPLINENODES","",0,0,0,0,0,mixture_composition).iUCode
    print('Nodes: ' + str(Nnodes))
    data = []
    for i in range(1,Nnodes+1):
        rho_vap,T,p,rho_eq,z0,z1 = rp.REFPROPdll("","SPLINENODES","",0,0,i,0,0,mixture_composition).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)
    # And plot them...
    plt.plot(df['T'],df.p,'o-')

    # Plot the specified state point at STP
    plt.plot(T_in, P_in,'o-',ms=8)

    plt.yscale('log')
    plt.xlabel('T [K]')
    plt.ylabel('p [kPa]')
    plt.show()


def print_fluid_info(fluid_name, rp):
    # Print out a range of properties for a fluid 
    z = [1.0] + [0.0]*19
    MOLAR_SI = rp.GETENUMdll(0,"MOLAR SI").iEnum
    print('Numerical outputs, stored in Output')
    print('-----------------------------------')
    for k in ["TC","PC","DC","TTRP","PTRP","MM","ACF",
              "DIPOLE","TNBP","REOS","ODP","GWP","TMIN",
              "TMAX","PMAX","HEATCOMB","HFRM"]:
        r = rp.REFPROPdll(fluid_name,"",k,MOLAR_SI,0,0,0,0,z)
        print(k + ' [' + r.hUnits+']:', r.Output[0], r.herr)

    print('\nString Outputs, stored in hUnits')
    print('--------------------------------')
    for k in ["ALTID","CAS#","CHEMFORM","FULLCHEMFORM","SYNONYM",
              "FAMILY","HASH","INCHI","INCHIKEY","LONGNAME",
              "SAFETY","NAME","UNNUMBER","DOI_EOS(1)" ,
              "WEB_EOS(1)","REFSTATE","FLDNAME","FDIR(1)"]:
        r = rp.REFPROPdll(fluid_name,"",k,MOLAR_SI,0,0,0,0,z)
        print(k + ':', r.hUnits)



if __name__ == "__main__":
    #Refprop setup
    rp = REFPROPFunctionLibrary(os.environ['RPPREFIX'])
    rp.SETPATHdll(os.environ['RPPREFIX'])
    #units = rp.GETENUMdll(0,'MASS BASE SI').iEnum
    units = rp.GETENUMdll(0,'MOLAR BASE SI').iEnum
    print('REFPROP Version: ' + str(rp.RPVersion()) )

    #Setting Peng-Robinson EOS to be used
    flags=rp.FLAGSdll('PR',3)
    print(flags)


    #Here we mix mmh and nto as liquids together. Starting first by determining the enthalpy of the liquids at their respective states
    P = 101.325 # kPa, pressure of the constant-pressure system
    Z_mix = 0.3  #Percentage of mmh and nto by the rule:  h_mix = Z_mix*h_mmh + (1-Z_mix)*h_nto
    mix = {'MMH': Z_mix, 'NTO': 1-Z_mix}

    #Determine input string based on Z_mix
    mixture_fluids = ' * '.join([str(key) for key in mix])
    mixture_composition = [mix[key] for key in mix]
    print('Mixture REFPROP name: ' + mixture_fluids)
    print('Mixture Composition: ' + str(mixture_composition))


    #-----Component definitions-----

    #Get pure mmh enthalpy 
    # @1atm, Tboil = 360.65 K
    T_mmh=290 #Kelvin
    print_fluid_info('MMH', rp)
    mmh_enthalpy = rp.REFPROPdll('MMH', 'PT', 'H', units, 0, 0, P, T_mmh, [1.0])
    print('Enthalpy of MMH at T=' + str(T_mmh) + ', and P=' + str(P) + ': ' + str(mmh_enthalpy.Output[0]) + ' J/mol')
    print(mmh_enthalpy)
    t_sat = rp.REFPROPdll('MMH', 'PQ', 'T', units, 0, 0, P, 0.0, [1.0]).Output[0]
    print('MMH Saturation temperature at P= ' + str(P) + ' kPa: ' + str(t_sat) + ' Kelvin')


    #Get pure nto enthalpy 
    # @1atm, Tboil = 294.25 K
    T_nto=250
    print_fluid_info('NTO', rp)
    nto_enthalpy = rp.REFPROPdll('NTO', 'PT', 'H', units, 0, 0, P, T_nto, [1.0])
    print('Enthalpy of NTO at T=' + str(T_nto) + ', and P=' + str(P) + ': ' + str(nto_enthalpy.Output[0]) + ' J/mol')
    t_sat = rp.REFPROPdll('NTO', 'PQ', 'T', units, 0, 0, P, 0.0, [1.0]).Output[0]
    print(nto_enthalpy)
    print('NTO Saturation temperature at P= ' + str(P) + ' kPa: ' + str(t_sat) + ' Kelvin')

    mixture_enthalpy = Z_mix*mmh_enthalpy.Output[0] + (1-Z_mix)*nto_enthalpy.Output[0]  #J/mol
    plot_mixture_phase_boundaries(mixture_fluids, mixture_composition, P, 300, rp)

    #Make a PH-flash calculation
    o = rp.REFPROPdll(mixture_fluids, 'PH', "H;T;Qmole", units, 0, 0, P, mixture_enthalpy, mixture_composition)
    print(o.ierr)
    H, T, Qmole = o.Output[0:3]
    print(H, T, Qmole)
    print('Mixture Enthalpy: ' + str(mixture_enthalpy) + ' J/mol')

MMH.FLD

MMH                !short name
60-34-4           !CAS number
Monomethylhydrazine             !full name
CH3NHNH2                  !chemical formula
x                  !synonym
100.               !molecular weight [g/mol]
200.               !triple point temperature [K]
300.               !normal boiling point [K]
400.               !critical temperature [K]
5000.              !critical pressure [kPa]
6.                 !critical density [mol/L]
0.2                !acentric factor
-1.                !dipole moment [Debye] (use -1 for unknown)
NBP                !default reference state
8.0                !version number


! Add the name of the fluid you are using on the first line, and update
! the other lines with the properties of the fluid.  Below add the molar mass,
! critical point, and acentric factor. A Cp0 equation is required if you are
! using caloric properties (H, Cp, Snd, etc.)


#EOS
PRT  translated Peng-Robinson equation
?LITERATURE REFERENCES \
!end of info section
200.               !lower temperature limit [K]
700.               !upper temperature limit [K]
100000.            !upper pressure limit [kPa]
50.                !maximum density [mol/L]
CPP                !pointer to Cp0 model
46.07              !molecular weight [g/mol]       |update these
0.425              !acentric factor                |
567                !critical temperature [K]       |
8349.18            !critical pressure [kPa]        |
3.690              !critical density [mol/L]       |
8.3144598          !gas constant [J/mol-K]
1                  !Number of parameters for translation
0.0                !value of translation parameter


#AUX               !auxiliary model specification
CPP  ideal gas heat capacity function
?LITERATURE REFERENCE \
?\
!end of info section
50.0               !lower temperature limit [K]
5000.0             !upper temperature limit [K]
0.0                !upper pressure limit [kPa]
0.0                !maximum density [mol/L]
1.0          8.3144598                 !reducing parameters for T, Cp0
  1  0    0  0    0  0  0              !Nterms:  polynomial, exponential, cosh, sinh
 3.5             0.0


@END
c        1         2         3         4         5         6         7         8
c2345678901234567890123456789012345678901234567890123456789012345678901234567890

NTO.FLD

NTO         !short name
10544-72-6           !CAS number
Nitrogen Tetroxide               !full name
N2O4                 !chemical formula
x                  !synonym
100.               !molecular weight [g/mol]
200.               !triple point temperature [K]
300.               !normal boiling point [K]
400.               !critical temperature [K]
5000.              !critical pressure [kPa]
6.                 !critical density [mol/L]
0.2                !acentric factor
-1.                !dipole moment [Debye] (use -1 for unknown)
NBP                !default reference state
8.0                !version number


! Add the name of the fluid you are using on the first line, and update
! the other lines with the properties of the fluid.  Below add the molar mass,
! critical point, and acentric factor. A Cp0 equation is required if you are
! using caloric properties (H, Cp, Snd, etc.)


#EOS
PRT  translated Peng-Robinson equation
?LITERATURE REFERENCES \
!end of info section
200.               !lower temperature limit [K]
700.               !upper temperature limit [K]
100000.            !upper pressure limit [kPa]
50.                !maximum density [mol/L]
CPP                !pointer to Cp0 model
92.011             !molecular weight [g/mol]       |update these
0.834              !acentric factor                |
431.4              !critical temperature [K]       |
10233.8            !critical pressure [kPa]        |
5.9558             !critical density [mol/L]       |
8.3144598          !gas constant [J/mol-K]
1                  !Number of parameters for translation
0.0                !value of translation parameter


#AUX               !auxiliary model specification
CPP  ideal gas heat capacity function
?LITERATURE REFERENCE \
?\
!end of info section
50.0               !lower temperature limit [K]
5000.0             !upper temperature limit [K]
0.0                !upper pressure limit [kPa]
0.0                !maximum density [mol/L]
1.0          8.3144598                 !reducing parameters for T, Cp0
  1  0    0  0    0  0  0              !Nterms:  polynomial, exponential, cosh, sinh
 3.5             0.0


@END
c        1         2         3         4         5         6         7         8
c2345678901234567890123456789012345678901234567890123456789012345678901234567890

[wandadars]

With thanks to our AI overlords, Chat GPT-4 found the error, I was using kPa for pressure inputs and not Pascals.

[ianhbell]

And does that work after the fix?

[wandadars]

No, but this is an issue with my understanding of how to go about:

1.) Getting REFPROPdll to perform a flash calculation.
2.) How to get from something like Qmole to molar fractions of each component in the liquid and vapor phases.

Liquid phase = {'MMH': 0.2, 'NTO': 0.8}
Vapor phase = {'MMH': 0.3, 'NTO': 0.7}
(I'm not entirely sure that the above situation is even something that is computable)

[wandadars]

I think I was able to get something to work. Seems like the best approach is to use the ABFLSHdll to perform a PH-flash, and depending on what the value of q is from that, I store the liquid and vapor component mole fractions. I also included a call to the REFPROPdll just as a sanity check, and it seems to give similar results. It does give a different mixture temperature and quality, but the value from the two are similar.

The only thing that worries me is that I saw a post on one of the REFPROP-issues issues page that talked about the difference between pure fluid enthalpies being computed with one reference state, and mixture enthalpies using a different reference state. I still haven't figured out how to handle that as I get my desired enthalpy for the PH-flash by simply weighting the pure-fluid enthalpies by a mole-fraction factor between zero and one, meant to represent an initial combination of the two components.

import numpy as np
import os
import pandas
import matplotlib.pyplot as plt

import ctREFPROP.ctREFPROP as ct
import refprop_vle_plots as refplots

def print_fluid_info(fluid_name, rp):
    # Print out a range of properties for a fluid 
    z = [1.0] + [0.0]*19
    MOLAR_SI = rp.GETENUMdll(0,"MOLAR SI").iEnum
    print('Numerical outputs, stored in Output')
    print('-----------------------------------')
    for k in ["TC","PC","DC","TTRP","PTRP","MM","ACF",
              "DIPOLE","TNBP","REOS","ODP","GWP","TMIN",
              "TMAX","PMAX","HEATCOMB","HFRM"]:
        r = rp.REFPROPdll(fluid_name,"",k,MOLAR_SI,0,0,0,0,z)
        print(k + ' [' + r.hUnits+']:', r.Output[0], r.herr)

    print('\nString Outputs, stored in hUnits')
    print('--------------------------------')
    for k in ["ALTID","CAS#","CHEMFORM","FULLCHEMFORM","SYNONYM",
              "FAMILY","HASH","INCHI","INCHIKEY","LONGNAME",
              "SAFETY","NAME","UNNUMBER","DOI_EOS(1)" ,
              "WEB_EOS(1)","REFSTATE","FLDNAME","FDIR(1)"]:
        r = rp.REFPROPdll(fluid_name,"",k,MOLAR_SI,0,0,0,0,z)
        print(k + ':', r.hUnits)


def setup_refprop():
    #Refprop setup
    rp = ct.REFPROPFunctionLibrary(os.environ['RPPREFIX'])  
    rp.SETPATHdll(os.environ['RPPREFIX'])
    units = rp.GETENUMdll(0,'MOLAR BASE SI').iEnum
    print('REFPROP Version: ' + str(rp.RPVersion()) )

    #Setting Peng-Robinson EOS for all fluids
    flags=rp.FLAGSdll('PR',3)
    print('Flags: ', flags.herr)
    return rp, units

def get_saturation_temperature(fluid_name, units,  P, rp):
    '''
    Returns the saturation temperature of a fluid at a given pressure in Kelvin 
    '''
    t_sat = rp.REFPROPdll(fluid_name, 'PQ', 'T', units, 0, 0, P, 0.0, [1.0]).Output[0]
    return t_sat

def get_enthalpy(fluid_name, units, T, P, rp):
    #Get pure fluid enthalpy 
    enthalpy = rp.REFPROPdll(fluid_name, 'PT', 'H', units, 0, 0, P, T,[1.0]+[0.0]*19) 
    return enthalpy.Output[0]

def get_mixture_saturation_temperature(mixture_fluids, mixture_composition, P, units, rp):
    '''
    Returns the saturation temperature of a mixture at a given pressure in Kelvin 
    '''
    t_sat = rp.REFPROPdll(mixture_fluids, 'PQ', 'T', units, 0, 0, P, 0.0, mixture_composition).Output[0]
    return t_sat

def get_mixture_saturation_pressure(mixture_fluids, mixture_composition, T, units, rp):
    '''
    Returns the saturation pressure of a mixture at a given temperature in Pa
    '''
    p_sat = rp.REFPROPdll(mixture_fluids, 'TQ', 'P', units, 0, 0, T, 0.0, mixture_composition).Output[0]
    return p_sat

def get_mixture_critical_properties(mixture_fluids, mixture_composition, units, rp):
    # Set up the mixture
    #ncc, hFiles, z, ierr, herr = rp.SETMIXdll(mixture_fluids, "HMX.BNC", "DEF")
    #if ierr != 0:
     #   print(f"Error in setting up the mixture: {herr}")

    # Set up splines for the mixture
    ierr, herr = rp.SATSPLNdll(mixture_composition)
    if ierr != 0:
        print(f"Error in setting up splines for the mixture: {herr}")

    # Get the critical properties of the mixture
    z = mixture_composition + [0.0] * (20 - len(mixture_composition))
    #output, _, _, ierr, herr = rp.ALLPROPSdll('TC;PC;DC', units, 0, 0, 1, 1, z) #Doesn't work
    output, _, _, ierr, herr = rp.ALLPROPS20dll('TC;PC;DC', units, 0, 0, 1, 1, z)

    if ierr != 0:
        print(f"Error in getting critical properties: {herr}")

    T_critical, P_critical, rho_critical = output[:3]

    return T_critical, P_critical

def get_vapor_liquid_composition(mixture_fluids, mixture_composition, P, h, rp):
    # Set the fluid mixture in REFPROP
    rp.SETMIXdll(mixture_fluids, "HMX.BNC", "DEF")

    # Calculate the properties at given temperature and pressure
    P_kPa = P/1000
    T,P,D,Dl,Dv,x,y,q,e,h,s,Cv,Cp,w,ierr,herr = rp.ABFLSHdll('PH',P_kPa,h,mixture_composition,0)
    print("T: ", T)
    print("D: ", D)
    print("Dl: ", Dl)
    print("Dv: ", Dv)
    print("x: ", x)
    print("y: ", y)
    print("q: ", q)
    print("e: ", e)
    print("s: ", s)
    print("Cv: ", Cv)
    print("Cp: ", Cp)

    # x and y are the liquid and vapor phase mole fractions, respectively
    print('Liquid mole fractions: ', x)
    print('Vapor mole fractions: ', y)
    print('quality: ', q)
    return x, y

def compute_equilibrium_composition(Z_mix, T_mmh, T_nto, P):
    '''
        Main function that calculates the enthalpy of a mixture of MMH and NTO at a given temperature and pressure
        Arguments are: 
        Z_mix: Mole fraction of MMH in the mixture
        T_mmh: Temperature of MMH in K
        T_nto: Temperature of NTO in K
        P: Pressure of the system in Pa
    '''
    rp, units = setup_refprop()

    #Here we mix mmh and nto as liquids together. Starting first by determining the enthalpy of the liquids at their respective states
    mix = {'MMH': Z_mix, 'NTO': 1-Z_mix} 

    #Determine REFPROP input string format based on Z_mix
    mixture_fluids = ' * '.join([str(key) for key in mix])
    mixture_composition = [mix[key] for key in mix]
    print('Mixture REFPROP name: ' + mixture_fluids)
    print('Mixture Composition: ' + str(mixture_composition))

    #-----Individual Component definitions-----
 
    #Get enthalpy of MMH at T and P
    # Empirical data to test with: @1atm, Tboil = 360.65 K
    #print_fluid_info('MMH', rp)
    #refplots.plot_pure_fluid_saturation_curve('MMH', rp)
    #refplots.plot_pure_fluid_enthalpy_saturation_curve('MMH', rp)
    mmh_enthalpy = get_enthalpy('MMH', units, T_mmh, P, rp) 
    print('Enthalpy of MMH at T=' + str(T_mmh) + ', and P=' + str(P) + ': ' + str(mmh_enthalpy) + ' J/mol')

    t_sat = get_saturation_temperature('MMH', units, P, rp)
    print('MMH Saturation temperature at P= ' + str(P) + ' kPa: ' + str(t_sat) + ' Kelvin')


    #Get enthalpy of NTO at T and P
    # Empirical data to test with: @1atm, Tboil = 294.25 K
    #print_fluid_info('NTO', rp)
    #refplots.plot_pure_fluid_saturation_curve('NTO', rp)
    #refplots.plot_pure_fluid_enthalpy_saturation_curve('NTO', rp)
    nto_enthalpy = get_enthalpy('NTO', units, T_nto, P, rp)
    print('Enthalpy of NTO at T=' + str(T_nto) + ', and P=' + str(P) + ': ' + str(nto_enthalpy) + ' J/mol')
    
    t_sat = get_saturation_temperature('NTO', units, P, rp)
    print('NTO Saturation temperature at P= ' + str(P) + ' kPa: ' + str(t_sat) + ' Kelvin')

    mixture_enthalpy = Z_mix*mmh_enthalpy + (1-Z_mix)*nto_enthalpy  #J/mol
    print('Mixture Enthalpy: ' + str(mixture_enthalpy) + ' J/mol')

    #refplots.plot_mixture_phase_boundaries(mixture_fluids, mixture_composition, P, 300, rp)

    #Make a PH-flash calculation
    P_kPa = P/1000
    T_mixture,P_mixture,D,Dl,Dv,x,y,q,e,h,s,Cv,Cp,w,ierr,herr = rp.ABFLSHdll('PH',P_kPa,mixture_enthalpy,mixture_composition,0)
    print('Results of P-H flash calculation(ABFLSHdll):')
    print('Mixture Enthalpy: ' + str(mixture_enthalpy) + ' J/mol')
    print('Mixture Temperature: ' + str(T_mixture) + ' K')
    print('Mixture Quality: ' + str(q))


    rp.REFPROPdll(mixture_fluids, 'SETREF', "H;T;Qmole", units, 0, 0, 2, 0, mixture_composition) #Set reference state for mixtures to be same as pure fluids
    o = rp.REFPROPdll(mixture_fluids, 'PH', "H;T;Qmole", units, 0, 1, P, mixture_enthalpy, mixture_composition)
    if o.ierr != 0:
        print(f"Error {o.ierr}: {o.herr}")
    H, T_mixture_refprop, Qmole = o.Output[0:3]
    print('Results of P-H flash calculation(REFPROPdll):')
    print('Mixture Enthalpy: ' + str(mixture_enthalpy) + ' J/mol')
    print('Mixture Temperature: ' + str(T_mixture_refprop) + ' K')
    print('Mixture Quality: ' + str(Qmole))

    # ABFLSHdll returns the following:
    #q > 0 and q < 1: indicates a 2-phase state
    #q < 0:	Subcooled (compressed) liquid
    #q = 0:	Saturated liquid
    #q = 1:	Saturated vapor
    #q > 1:	Superheated vapor
    #q = -998:	Subcooled liquid, but quality not defined (usually P > Pc)
    #q = 998:	Superheated vapor, but quality not defined (usually T > Tc)
    #q = 999:	Supercritical state (T>Tc and P>Pc)
    if 0 < q < 1:  # Two-phase region
        print("Mixture is two-phase")
        liquid_composition, vapor_composition = get_vapor_liquid_composition(mixture_fluids, mixture_composition, P, mixture_enthalpy, rp)
        
        #create a dictionary with the keys being the values in mixture_fluids
        liquid_composition = dict(zip(mix.keys(), liquid_composition))
        vapor_composition = dict(zip(mix.keys(), vapor_composition))
    elif q <= 0 or q == -998:
        if q< 0:
            print("Mixture is subcooled liquid")
        elif q == 0:
            print("Mixture is saturated liquid")
        elif q == -998:
            print("Mixture is subcooled liquid, but quality not defined")

        liquid_composition = dict(zip(mix.keys(), mixture_composition))
        vapor_composition = dict(zip(mix.keys(), [0]*len(mixture_composition)))
    elif q >=1 or q == 998:
        if q > 1:
            print("Mixture is superheated vapor")
        elif q == 1:
            print("Mixture is saturated vapor")
        elif q == 998:
            print("Mixture is superheated vapor, but quality not defined")

        liquid_composition = dict(zip(mix.keys(), [0]*len(mixture_composition)))
        vapor_composition = dict(zip(mix.keys(), mixture_composition))

    print(f"Liquid phase composition: {liquid_composition}")
    print(f"Vapor phase composition: {vapor_composition}")

    return liquid_composition, vapor_composition

if __name__ == "__main__":
# Loop and call main() and pass it a series of Z_mix values from 0 to 1.0
    T_mmh = 350 # Kelvin
    T_nto = 300 # Kelvin
    P = 101325 # Pa, pressure of the constant-pressure system
    num_z_mix = 30 # Number of Z_mix values to loop over
    results = []
    for Z_mix in np.linspace(0, 1.0, num_z_mix):
        print('\n\nSolving for Z_mix = ' + str(Z_mix))
        liquid_composition, vapor_composition = compute_equilibrium_composition(Z_mix, T_mmh, T_nto, P)
        print('Liquid Composition: ' + str(liquid_composition))
        print('Vapor Composition: ' + str(vapor_composition))
        
        # Store the results dictionary and keep them for all Z_mix values.
        results.append({
            'Z_mix': Z_mix,
            'T_mmh': T_mmh,
            'T_nto': T_nto,
            'P': P,
            'liquid_composition': liquid_composition,
            'vapor_composition': vapor_composition
        })

# Print out a table summarizing the results, and make it look nice.
header = "{:<8}{:<8}{:<8}{:<12}{:<20}{:<20}".format("Z_mix", "T_mmh", "T_nto", "P", "Liquid Composition", "Vapor Composition")
print(header)
print("-" * len(header))

for result in results:
    row = "{:<8.2f}{:<8}{:<8}{:<12}{:<20}{:<20}".format(
        result['Z_mix'],
        result['T_mmh'],
        result['T_nto'],
        result['P'],
        str(result['liquid_composition']),
        str(result['vapor_composition'])
    )
    print(row)

[ianhbell]

If the inputs are homogeneous, then the bulk composition is the only one that is meaningful. If the inputs are not homogeneous, you can query the output compositions for each phase from the call to ABFLSHdll.

[wandadars]

Can you help me understand what you mean by homogeneous inputs and bulk composition?

[ianhbell]

At a given statepoint, there are two "normal" things that can happen when you specify a set of two state variables (more things that can happen when multiple phases can be in equilibria). But let's assume we have either a homogenous phase in which you cannot distinguish liquid and vapor phases, or a phase split, in which the liquid and vapor phases have different compositions and densities (but the same pressure and temperature). In every calculation in REFPROP (with few exceptions), what you specify is the composition in the bulk -- what you would get if you massed out a mixture of given composition into a cylinder.