https://github.com/CoolProp/CoolProp
Tip revision: 84897ce7a110a70ba2229e7a1f5d17ea03efd068 authored by Jorrit Wronski on 16 December 2022, 12:59:36 UTC
use no version number at all to fix the source build
use no version number at all to fix the source build
Tip revision: 84897ce
ideal_curves.py
import numpy as np
import matplotlib.pyplot as plt
import CoolProp, scipy.optimize
class CurveTracer(object):
def __init__(self, backend, fluid, p0, T0):
"""
p0 : Initial pressure [Pa]
T0 : Initial temperatrure [K]
"""
self.P = [p0]
self.T = []
self.AS = CoolProp.AbstractState(backend, fluid)
# Solve for Temperature for first point
T = scipy.optimize.newton(self.objective_T, T0, args=(p0, -1))
self.T.append(T)
def objective_T(self, T, p, rho_guess):
""" Base class function """
if rho_guess < 0:
self.AS.update(CoolProp.PT_INPUTS, p, T)
else:
guesses = CoolProp.CoolProp.PyGuessesStructure()
guesses.rhomolar = rho_guess
self.AS.update_with_guesses(CoolProp.PT_INPUTS, p, T, guesses)
return self.objective()
def TPcoords(self, t, lnT, lnp, rlnT=0.1, rlnp=0.1):
return np.exp(lnT + rlnT * np.cos(t)), np.exp(lnp + rlnp * np.sin(t))
def obj_circle(self, t, lnT, lnp):
T2, P2 = self.TPcoords(t, lnT, lnp)
self.AS.update(CoolProp.PT_INPUTS, P2, T2)
r = self.objective()
return r
def trace(self):
t = self.starting_direction()
for i in range(1000):
try:
lnT = np.log(self.T[-1])
lnp = np.log(self.P[-1])
t = scipy.optimize.brentq(self.obj_circle, t - np.pi / 2, t + np.pi / 2, args=(lnT, lnp))
T2, P2 = self.TPcoords(t, lnT, lnp)
self.T.append(T2)
self.P.append(P2)
if self.T[-1] < self.AS.keyed_output(CoolProp.iT_triple) or self.P[-1] > 1000 * self.AS.keyed_output(CoolProp.iP_critical):
break
except ValueError as VE:
print(VE)
break
return self.T, self.P
class IdealCurveTracer(CurveTracer):
def __init__(self, *args, **kwargs):
CurveTracer.__init__(self, *args, **kwargs)
def objective(self):
""" Z = 1 """
return self.AS.keyed_output(CoolProp.iZ) - 1
def starting_direction(self):
""" Start searching directly up ( or calculate as orthogonal to gradient ) """
return np.pi / 2.0
class BoyleCurveTracer(CurveTracer):
def __init__(self, *args, **kwargs):
CurveTracer.__init__(self, *args, **kwargs)
def objective(self):
""" dZ/dv|T = 0 """
r = (self.AS.p() - self.AS.rhomolar() * self.AS.first_partial_deriv(CoolProp.iP, CoolProp.iDmolar, CoolProp.iT)) / (self.AS.gas_constant() * self.AS.T())
# print self.AS.T(), self.AS.p(), r
return r
def starting_direction(self):
""" Start searching directly up """
return np.pi / 2.0
class JouleInversionCurveTracer(CurveTracer):
def __init__(self, *args, **kwargs):
CurveTracer.__init__(self, *args, **kwargs)
def objective(self):
""" dZ/dT|v = 0 """
r = (self.AS.gas_constant() * self.AS.T() * 1 / self.AS.rhomolar() * self.AS.first_partial_deriv(CoolProp.iP, CoolProp.iT, CoolProp.iDmolar) - self.AS.p() * self.AS.gas_constant() / self.AS.rhomolar()) / (self.AS.gas_constant() * self.AS.T())**2
# print self.AS.T(), self.AS.p(), r
return r
def starting_direction(self):
""" Start searching directly up """
return np.pi / 2.0
class JouleThomsonCurveTracer(CurveTracer):
def __init__(self, *args, **kwargs):
CurveTracer.__init__(self, *args, **kwargs)
def objective(self):
""" dZ/dT|p = 0 """
dvdT__constp = -self.AS.first_partial_deriv(CoolProp.iDmolar, CoolProp.iT, CoolProp.iP) / self.AS.rhomolar()**2
r = self.AS.p() / (self.AS.gas_constant() * self.AS.T()**2) * (self.AS.T() * dvdT__constp - 1 / self.AS.rhomolar())
# print self.AS.T(), self.AS.p(), r
return r
def starting_direction(self):
""" Start searching directly up """
return np.pi / 2.0
backend = 'HEOS'
fluid = 'R125'
kwargs = dict(lw=2)
print('Ideal')
ICT = IdealCurveTracer(backend, fluid, p0=1e5, T0=900)
T, p = ICT.trace()
plt.plot(T, p, '-', label='Ideal Curve', **kwargs)
print('Boyle')
BCT = BoyleCurveTracer(backend, fluid, p0=1e5, T0=800)
T, p = BCT.trace()
plt.plot(T, p, '-', label='Boyle Curve', **kwargs)
print('Joule Inversion')
JIT = JouleInversionCurveTracer(backend, fluid, p0=1e5, T0=1800)
T, p = JIT.trace()
plt.plot(T, p, '-', label='Joule Inversion Curve', **kwargs)
print('Joule-Thomson')
JTCT = JouleThomsonCurveTracer(backend, fluid, p0=1e5, T0=1800)
T, p = JTCT.trace()
plt.plot(T, p, '-', label='Joule-Thomson Curve', **kwargs)
print('Saturation Curve')
Tt = ICT.AS.keyed_output(CoolProp.iT_triple)
Tc = ICT.AS.keyed_output(CoolProp.iT_critical)
Ts = np.linspace(Tt, Tc - 1.e-6)
ps = CoolProp.CoolProp.PropsSI('P', 'T', Ts, 'Q', 0, backend + '::' + fluid)
plt.plot(Ts, ps, '-', label='Saturation Curve', **kwargs)
plt.yscale('log')
plt.xscale('log')
plt.xlabel('T (K)')
plt.ylabel('p (Pa)')
plt.legend(loc='best')
plt.savefig('IdealCurves.png')
plt.show()