BioSTEAMDevelopmentGroup / thermosteam

BioSTEAM's Premier Thermodynamic Engine
Other
57 stars 12 forks source link

CoolProp compatibility bug using liquid-vapour-equilibrium (vle) #54

Closed nkawerau closed 2 years ago

nkawerau commented 2 years ago

Describe the bug Running a Thermosteam code in an environment with a clean thermosteam installation (pip install thermosteam) works fine. As soon as I install CoolProp (pip install coolprop) running the same code raises errors.

To Reproduce

import thermosteam as tmo
from thermo import electrochem as ec

def Vl(mol, T, P=None):  # m3
    MWs = mixture.MWs[H2O_index:KOH_index+1]
    wt = MWs * mol[H2O_index:KOH_index+1]
    wt_no_water = list(wt)
    wt_no_water.pop(0)
    rho = ec.Laliberte_density(T, wt_no_water / wt.sum(), CASs_no_water)  # kg / m3
    MW = wt.sum() / mol[H2O_index:KOH_index+1].sum()
    return tmo.functional.rho_to_V(rho, MW)

def mul(mol, T, P=None):  # Pa*s
    MWs = mixture.MWs[H2O_index:KOH_index+1]
    wt = MWs * mol[H2O_index:KOH_index+1]
    wt_no_water = list(wt)
    wt_no_water.pop(0)
    mu = ec.Laliberte_viscosity(T, wt_no_water / wt.sum(), CASs_no_water)
    return mu

def Cnl(mol, T, P=None):  # kJ / kg*K
    MWs = mixture.MWs[H2O_index:KOH_index+1]
    wt = MWs * mol[H2O_index:KOH_index+1]
    wt_no_water = list(wt)
    wt_no_water.pop(0)
    Cp = ec.Laliberte_heat_capacity(
        T, wt_no_water / wt.sum(), CASs_no_water
    )  # J / kg / K
    MW = wt.sum() / mol.sum()
    return (
        Cp * MW / 1000.0 * mol[H2O_index:KOH_index+1].sum()
    )  # Special case, needs to return the total capacity

# according to Gilliam et al., 2006, A review  of specific conductivities of
# potassium hydroxide solutions for various concentrations and temperatures
def electrical_conductivity_KOH(V, mol, T):
    n_KOH = mol[chemicals.index("KOH")]
    molar_concentration = n_KOH / V
    temperature = T

    A = -2.041
    B = -0.0028
    C = 0.005332
    D = 207.2
    E = 0.001043
    F = -0.0000003

    # [Sm/cm]
    kappa = (
        A * molar_concentration
        + B * pow(molar_concentration, 2)
        + C * molar_concentration * temperature
        + D * (molar_concentration / temperature)
        + E * pow(molar_concentration, 3)
        + F * (pow(molar_concentration, 2) * pow(temperature, 2))
    )

    return kappa * 100 # [Sm/m]

"""stream definition"""

chemicals = tmo.Chemicals([])
chemicals.append(tmo.Chemical('H2', phase="g"))
chemicals.append(tmo.Chemical('H2O'))
chemicals.append(tmo.Chemical('KOH', phase="l"))
chemicals.append(tmo.Chemical('O2', phase="g"))

mixture = tmo.Mixture.from_chemicals(chemicals)
thermo = tmo.Thermo(chemicals, mixture, cache=False)

H2O_index = chemicals.index("H2O")
KOH_index = chemicals.index("KOH")
CASs_no_water = list(chemicals.CASs[H2O_index+1:KOH_index+1])

"""implementing new mixture rules and functions"""
mixture.rule = "electrolyte"
tmo.settings.set_thermo(thermo)

replace = object.__setattr__
replace(mixture.V, "l", Vl)
replace(mixture.mu, "l", mul)
replace(mixture.Cn, "l", Cnl)

stream = tmo.MultiStream()

stream.set_flow(1, "kg/hr", ("g", "H2"))
stream.set_flow(0, "kg/hr", ("g", "O2"))
stream.set_flow(3, "kg/hr", ("l", "KOH"))
stream.set_flow(7, "kg/hr", ("l", "H2O"))

stream.set_property("T", 25, "degC")
stream.set_property("P", 1, "bar")

stream.vle(T=stream.T, P=stream.P)

stream.show(flow='kg/hr')

viscosity_fluid = stream.mu #[Pa/s]
specific_enthalpy = (stream.H / stream.F_mass)*1000 #[J/kg]
specific_volume = stream.F_vol / stream.F_mass #[m3/kg]
density = stream.rho #[m3/kg]
specific_entropy = (stream.S / stream.F_mass)*1000 #[J/(K*kg)]

print("thermosteam_testing specific enthalpy:", specific_enthalpy)
print("thermosteam_testing specific volume:", specific_volume)
print("thermosteam_testing density:", density)
print("thermosteam_testing specific entropy:", specific_entropy)
print("thermosteam_testing viscosity:", viscosity_fluid)

Expected behavior no runtime error

Actual behavior MultiStream: s1 phases: ('g', 'l'), T: 298.15 K, P: 100000 Pa flow (kg/hr): (g) H2 1 H2O 0.292 (l) H2O 6.71 KOH 3 Traceback (most recent call last): File "C:\Users\nkawerau\AppData\Roaming\Python\Python38\site-packages\IPython\core\interactiveshell.py", line 2632, in safe_execfile py3compat.execfile( File "C:\Users\nkawerau\AppData\Roaming\Python\Python38\site-packages\IPython\utils\py3compat.py", line 55, in execfile exec(compiler(f.read(), fname, "exec"), glob, loc) File "C:\Users\nkawerau\Documents\python-scripts\ael-model\thermosteam_testing.py", line 102, in viscosity_fluid = stream.mu #[Pa/s] File "C:\Users\nkawerau\Anaconda3\envs\tespy_thermosteam_env\lib\site-packages\thermosteam_multi_stream.py", line 534, in mu self._property_cache['mu'] = mu = self.mixture.xmu(self._imol.iter_composition(), self._thermal_condition) File "C:\Users\nkawerau\Anaconda3\envs\tespy_thermosteam_env\lib\site-packages\thermosteam\mixture\mixture.py", line 368, in xmu return sum([self.mu(phase, mol, T, P) for phase, mol in phase_mol]) File "C:\Users\nkawerau\Anaconda3\envs\tespy_thermosteam_env\lib\site-packages\thermosteam\mixture\mixture.py", line 368, in return sum([self.mu(phase, mol, T, P) for phase, mol in phase_mol]) File "C:\Users\nkawerau\Anaconda3\envs\tespy_thermosteam_env\lib\site-packages\thermosteam\base\phase_handle.py", line 84, in call return getattr(self, phase)(z, T, P) File "C:\Users\nkawerau\Anaconda3\envs\tespy_thermosteam_env\lib\site-packages\thermosteam\mixture\ideal_mixture_model.py", line 56, in call return sum([j i(T, P) for i, j in zip(self.models, mol) if j]) File "C:\Users\nkawerau\Anaconda3\envs\tespy_thermosteam_env\lib\site-packages\thermosteam\mixture\ideal_mixture_model.py", line 56, in return sum([j * i(T, P) for i, j in zip(self.models, mol) if j]) File "C:\Users\nkawerau\Anaconda3\envs\tespy_thermosteam_env\lib\site-packages\thermosteam\thermo\tp_dependent_property.py", line 21, in call return self.TP_dependent_property(T, P) File "C:\Users\nkawerau\AppData\Roaming\Python\Python38\site-packages\thermo\utils\tp_dependent_property.py", line 245, in TP_dependent_property raise RuntimeError("%s method '%s' is not valid at T=%s K and P=%s Pa for component with CASRN '%s'" %(self.name, method_P, T, P, self.CASRN)) RuntimeError: Gas viscosity method 'COOLPROP' is not valid at T=298.15 K and P=99999.99999999999 Pa for component with CASRN '7732-18-5'

Version Python: 3.8.12 Thermosteam: 0.27.19 Thermo: 0.2.13 CoolProp: 6.4.1

yoelcortes commented 2 years ago

Hi! Thanks for posting this issue. I adjusted the ranked methods in thermosteam to use CoolProp only if its the only method available (please update or pull master branch). Note that coolprop uses equations of state, so having water as a gas at STP will lead to errors.

With the new update, you can get your code working by add this line:

chemicals.H2O.mu.g.method_P = None # Before any calculations

Then you should get the following:

MultiStream: s1
 phases: ('g', 'l'), T: 298.15 K, P: 100000 Pa
 flow (kg/hr): (g) H2   1
                   H2O  0.292
               (l) H2O  6.71
                   KOH  3
thermosteam_testing specific enthalpy: 64797.7223270309
thermosteam_testing specific volume: 1.1554499145757928
thermosteam_testing density: 0.8654637361474349
thermosteam_testing specific entropy: 8936.728801346717
thermosteam_testing viscosity: 0.002062542267450246