Skip to content
Snippets Groups Projects
Commit 29b34603 authored by Ian Bell's avatar Ian Bell
Browse files

Updated isochoric to to make it more like a reference

parent f96ffb72
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 110 additions and 307 deletions
notebooks/isochoric.py
+ 110
− 307
View file @ 29b34603
# Standard libraries
import timeit, sys
# Scipy stack packages
import numpy as np
import pandas
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
from matplotlib.backends.backend_pdf import PdfPages
import VLEIsoTracer as vle
# print('VLEIsoTracer is located at:', vle.__file__)
IMPOSED_T = vle.VLEIsolineTracer.imposed_variable.IMPOSED_T
# Some integration tools from PDSim
from PDSIMintegrators import AbstractRK45ODEIntegrator, AbstractSimpleEulerODEIntegrator
sys.path.append('../bld/Release')
# Special packages
import teqp
import CoolProp.CoolProp as CP
def get_drhovecdp_sat(model, T, rhovecL, rhovecV):
tic = timeit.default_timer()
Hliq = teqp.build_Psi_Hessian_autodiff(model, T, rhovecL)
Hvap = teqp.build_Psi_Hessian_autodiff(model, T, rhovecV)
Hvap[~np.isfinite(Hvap)] = 1e20
Hliq[~np.isfinite(Hliq)] = 1e20
N = len(rhovecL)
A = np.zeros((N, N))
b = np.ones((N, 1))
assert(len(rhovecL)==len(rhovecV))
if np.all(rhovecL != 0) and np.all(rhovecV != 0):
A[0,0] = np.dot(Hliq[0,:], rhovecV)
A[0,1] = np.dot(Hliq[1,:], rhovecV)
A[1,0] = np.dot(Hliq[0,:], rhovecL)
A[1,1] = np.dot(Hliq[1,:], rhovecL)
drhodp_liq = np.linalg.solve(A, b).squeeze()
drhodp_vap = np.linalg.solve(Hvap, np.dot(Hliq,drhodp_liq))
return drhodp_liq, drhodp_vap
else:
# Special treatment for infinite dilution
murL = teqp.build_Psir_gradient_autodiff( model, T, rhovecL)
murV = teqp.build_Psir_gradient_autodiff( model, T, rhovecV)
RL = model.get_R(rhovecL/rhovecL.sum())
RV = model.get_R(rhovecV/rhovecV.sum())
# First, for the liquid part
for i in range(N):
for j in range(N):
if rhovecL[j] == 0:
# Analysis is special if j is the index that is a zero concentration. If you are multiplying by the vector
# of liquid concentrations, a different treatment than the case where you multiply by the vector
# of vapor concentrations is required
# ...
# Initial values
Aij = Hliq[j,:]*(rhovecV if i == 0 else rhovecL) # coefficient-wise product
# A throwaway boolean for clarity
is_liq = (i==1)
# Apply correction to the j term (RT if liquid, RT*phi for vapor)
Aij[j] = RL*T if is_liq else RL*T*np.exp(-(murV[j] - murL[j])/(RL*T))
# Fill in entry
A[i, j] = Aij.sum()
else:
# Normal
A[i, j] = np.dot(Hliq[j,:], rhovecV if (i == 0) else rhovecL)
drhodp_liq = np.linalg.solve(A, b)
# Then, for the vapor part, also requiring special treatment
# Left-multiplication of both sides of equation by diagonal matrix with liquid concentrations along diagonal, all others zero
diagrhovecL = np.diagflat(rhovecL)
PSIVstar = diagrhovecL@Hvap
PSILstar = diagrhovecL@Hliq
for j in range(N):
if rhovecL[j] == 0:
PSILstar[j, j] = RL*T
PSIVstar[j, j] = RV*T/np.exp(-(murV[j] - murL[j])/(RV*T))
drhodp_vap = np.linalg.solve(PSIVstar, PSILstar@drhodp_liq)
toc = timeit.default_timer()
# print(toc-tic)
return drhodp_liq.squeeze(), drhodp_vap.squeeze()
def get_dxdp(*, rhovec, drhovecdp):
rhomolar = np.sum(rhovec)
molefrac = rhovec/rhomolar
......@@ -90,145 +28,139 @@ def get_drhovecdx(*, i, model, T, rhovec):
dpdxL = 1.0/get_dxdp(rhovec=rhovecL, drhovecdp=drhovecdpL)
return np.array((drhovecdpL*dpdxL[i]).tolist()+(drhovecdpV*dpdxL[i]).tolist())
def traceT(model, T, rhovec0, *, kwargs={}):
def traceT_x(model, T, rhovec0, *, kwargs={}):
R = model.get_R(np.array([0.0, 1.0]))
def rhovecprime(x, rhovec):
print(x)
return get_drhovecdx(i=0, model=model, T=T, rhovec=rhovec)
tic = timeit.default_timer()
sol = solve_ivp(rhovecprime, [0.0, 0.5], y0=rhovec0, method='RK45',dense_output=True, t_eval = np.linspace(0,0.49,1000), rtol=1e-8)
sol = solve_ivp(rhovecprime, [1.0, 0.2], y0=rhovec0, method='RK45',dense_output=True, t_eval = np.linspace(1,0.21,1000), rtol=1e-8)
p = [teqp.get_pr(model, T, sol.y[0:2,j]) + R*T*sol.y[0:2,j].sum() for j in range(len(sol.t))]
toc = timeit.default_timer()
print(toc-tic)
x = sol.y[0,:]/(sol.y[0,:] + sol.y[1,:])
y = sol.y[2,:]/(sol.y[2,:] + sol.y[3,:])
plt.plot(sol.t, p, 'r',**kwargs)
plt.plot(x, p, 'r',**kwargs)
plt.plot(y, p, 'g',**kwargs)
def traceT_arclength(model, T, rhovec0, *, kwargs={}):
R = model.get_R(np.array([0.0, 1.0]))
print(T)
c = 1.0
def norm(x): return (x*x).sum()**0.5
def rhovecprime(t, rhovec):
tic = timeit.default_timer()
drhovecdpL, drhovecdpV = teqp.get_drhovecdp_Tsat(model, T, rhovec[0:2], rhovec[2::])
# drhovecdpL, drhovecdpV = tracer.get_drhovecdp_sat(rhovec[0:2], rhovec[2::], T)[0:2]
# print(, drhovecdpL, drhovecdpV)
toc = timeit.default_timer()
# print(toc-tic)
def norm(x):
return (x*x).sum()**0.5
dpdt = (norm(drhovecdpL) + norm(drhovecdpV))**-0.5
der = np.zeros_like(rhovec)
der[0:2] = (drhovecdpL*dpdt).squeeze()
der[2:4] = (drhovecdpV*dpdt).squeeze()
# print(rhovec, rhovec[0]+rhovec[1], rhovec[2]+rhovec[3])
# print(t, rhovec, rhovec[0]+rhovec[1], rhovec[2]+rhovec[3], rhovec[0]/(rhovec[0]+rhovec[1]))
# if any(~np.isfinite(rhovec)):
# raise ValueError()
return c*der
class TestIntegrator(object):
"""
Implements the functions needed to satisfy the ABC requirements
"""
def __init__(self, Nmax):
self.x, self.y = [], []
self.Itheta = 0
self.Nmax = Nmax
self.termination_reason = None
def post_deriv_callback(self): pass
def premature_termination(self):
if self.Itheta > self.Nmax:
self.termination_reason = 'too many steps'
return True
if np.any(self.xold<0):
self.termination_reason = 'negative x'
return True
if np.any(~np.isfinite(self.xold)):
self.termination_reason = 'nan value of x'
return True
return False
def get_initial_array(self):
return rhovec0.copy()
def pre_step_callback(self):
if self.Itheta == 0:
self.x.append(0)
self.y.append(self.get_initial_array().tolist())
def post_step_callback(self):
self.x.append(self.t0_cache)
self.y.append(self.xold_cache.tolist())
# Polish the solution...
def derivs(self, t0, xold):
self.t0_cache = t0
self.xold_cache = xold.copy()
return rhovecprime(t0, xold)
class EulerIntegrator(TestIntegrator, AbstractSimpleEulerODEIntegrator):
""" Mixin class using the functions defined in TestIntegrator """
pass
class RK45Integrator(TestIntegrator, AbstractRK45ODEIntegrator):
""" Mixin class using the functions defined in TestIntegrator """
pass
der_init = rhovecprime(0, rhovec0)
if np.any(rhovec0 + der_init*1e-6 < 0):
c *= -1
# print('init', der_init, T)
print('flip c', der_init, rhovec0)
events = [lambda t,x: x[0], lambda t,x: x[1],
lambda t,x: x[2], lambda t,x: x[3]]
events.append(lambda t, z: ((z[0]+z[1])/(z[2]+z[3])-1)-0.2)
# class StepCounter():
# def __init__(self, Nmax):
# self.counter = 0
# self.Nmax = Nmax
# def __call__(self, t, z):
# self.counter += 1
# v = self.Nmax-self.counter+0.5
# return v
# events.append(StepCounter(200))
for e in events:
e.direction = -1
e.terminal = True
# tic = timeit.default_timer()
# sol = solve_ivp(rhovecprime, [0.0, 200000], y0=rhovec0.copy(), method='RK45',dense_output=True, events=events, rtol=1e-8, atol=1e-10)
# # print(sol.t_events, sol.y_events)
# # print(sol.message)
# print(sol.t, sol.nfev)
# toc = timeit.default_timer()
# print(toc-tic,'s for parametric w/ solve_ivp')
# t = np.linspace(np.min(sol.t), np.max(sol.t), 10000)
# soly = sol.sol(t)
tic = timeit.default_timer()
sol = solve_ivp(rhovecprime, [0.0, 3500000], y0=rhovec0.copy(), method='RK45',dense_output=True, events=events)#, rtol=1e-8)
# print(sol.t_events, sol.y_events)
# print(sol.message)
i = RK45Integrator(Nmax=1000)
i.do_integration(0, 5000000, atol=1e-6, rtol=1e-6, hmin=4.46772051e-03)
t = i.x
soly = np.array(i.y).T
toc = timeit.default_timer()
print(toc-tic,'s for parametric')
t = np.linspace(np.min(sol.t), np.max(sol.t), 1000)
soly = sol.sol(t)
print(toc-tic,'s for parametric w/ my RK45')
reason = i.termination_reason
if reason:
print('reason:', reason)
print(i.Itheta)
x = soly[0,:]/(soly[0,:] + soly[1,:])
y = soly[2,:]/(soly[2,:] + soly[3,:])
p = [teqp.get_pr(model, T, soly[0:2,j]) + R*T*soly[0:2,j].sum() for j in range(len(t))]
plt.plot(x, p, 'r', **kwargs)
plt.plot(y, p, 'g', **kwargs)
# print(x,p)
def drhovecdT_isopleth(model, T, rhovec):
assert(len(rhovec)==4)
rhovecL = rhovec[0:2]
rhovecV = rhovec[2::]
Hliq = teqp.build_Psi_Hessian_autodiff(model, T, rhovecL)
Hvap = teqp.build_Psi_Hessian_autodiff(model, T, rhovecV)
rhoL = sum(rhovecL)
rhoV = sum(rhovecV)
xL = rhovecL/rhoL
xV = rhovecV/rhoV
R = model.get_R(xL)
deltas = teqp.get_dchempotdT_autodiff(model, T, rhovecV) - teqp.get_dchempotdT_autodiff(model, T, rhovecL)
deltarho = rhovecV-rhovecL
dpdTV = R*rhoV*(1+model.get_Ar01(T, rhoV, xV)-model.get_Ar11(T, rhoV, xV))
dpdTL = R*rhoL*(1+model.get_Ar01(T, rhoL, xL)-model.get_Ar11(T, rhoL, xL))
deltabeta = dpdTV-dpdTL
drhovecdTL = (np.dot(deltas,rhovecV)-deltabeta)/np.dot(Hliq@deltarho, xL)*xL
drhovecdTV = np.linalg.solve(Hvap, Hliq@drhovecdTL-deltas)
return np.array(drhovecdTL.tolist()+drhovecdTV.tolist())
def robust_VLE(model, T, Tc, rhoc):
DeltaT = 1e-2*Tc
Tg = Tc-DeltaT
if Tg < 0:
raise ValueError("T > Tc")
rholiq, rhovap = teqp.extrapolate_from_critical(model, Tc, rhoc, Tg)
if rholiq < 0 or rhovap < 0:
raise ValueError("Negative density obtained from critical extrapolation")
rholiq, rhovap = teqp.pure_VLE_T(model, Tg, rholiq, rhovap, 100)
z = np.array([1.0])
R = model.get_R(z)
def dpdTsat(model, T, rhoL, rhoV):
""" dp/dT = ((hV-hL)/T)/(vV-vL)
The ideal parts of enthalpy cancel, leaving just the residual parts
"""
numV = R*(model.get_Ar01(T, rhoV, z) + model.get_Ar10(T,rhoV,z))
numL = R*(model.get_Ar01(T, rhoL, z) + model.get_Ar10(T,rhoL,z))
num = numV-numL
den = 1/rhoV-1/rhoL
return num/den
def get_drhodT_sat(model, T, rhoL, rhoV, Q):
dpdT = dpdTsat(model, T, rhoL, rhoV)
rho = rhoL if Q == 0 else rhoV
Ar01 = model.get_Ar01(T, rho, z)
Ar02 = model.get_Ar02n(T, rho, z)[2]
Ar11 = model.get_Ar11(T, rho, z)
pL = [teqp.get_pr(model, T, soly[0:2,j]) + R*T*soly[0:2,j].sum() for j in range(len(t))]
pV = [teqp.get_pr(model, T, soly[2:4,j]) + R*T*soly[2:4,j].sum() for j in range(len(t))]
plt.plot(x, pL, 'r', ms=0.1, **kwargs)
plt.plot(y, pV, 'g', ms=0.1, **kwargs)
# plt.plot(y, p, 'g', **kwargs)
dpdrho__T = R*T*(1 + 2*Ar01 + Ar02)
dpdT__rho = R*rho*(1 + Ar01 - Ar11)
drhodT__p = -dpdT__rho/dpdrho__T
drhodp__T = 1/dpdrho__T
return drhodT__p + drhodp__T*dpdT
def rhoprime(T, rhoLrhoV):
rhoL, rhoV = rhoLrhoV
drhoLdT = get_drhodT_sat(model, T, rhoL, rhoV, 0)
drhoVdT = get_drhodT_sat(model, T, rhoL, rhoV, 1)
return np.array([drhoLdT, drhoVdT])
sol = solve_ivp(rhoprime, [Tg, T], y0=[rholiq, rhovap], method='RK45', dense_output=False, rtol=1e-8)
rholiq, rhovap = sol.y[:,-1]
rholiq, rhovap = teqp.pure_VLE_T(model, T, rholiq, rhovap, 100)
return rholiq, rhovap
def prettyPX(model, puremodels, ipure, Tvec):
def prettyPX(model, puremodels, names, ipure, Tvec, *, plot_critical=True, show=True, ofname=None):
"""
ipure: int
index (0-based) of the pure fluid model to use
......@@ -237,7 +169,7 @@ def prettyPX(model, puremodels, ipure, Tvec):
tic = timeit.default_timer()
Tcvec = model.get_Tcvec()
rhocvec = 1/model.get_vcvec()
k = 1 # Have to start at pure second component for now... In theory either are possible.
k = 1
T0 = Tcvec[k]
rho0 = rhocvec
rho0[1-k] = 0
......@@ -247,6 +179,7 @@ def prettyPX(model, puremodels, ipure, Tvec):
df = pandas.DataFrame(curveJSON)
df['z0 / mole frac.'] = df['rho0 / mol/m^3']/(df['rho0 / mol/m^3']+df['rho1 / mol/m^3'])
plt.plot(df['z0 / mole frac.'], df['p / Pa'], lw=2)
#print(toc-tic, 'complete critical locus')
def plot_isotherms():
for T in Tvec:
......@@ -265,149 +198,19 @@ def prettyPX(model, puremodels, ipure, Tvec):
except BaseException as BE:
print(BE)
# raise
plot_critical_curve()
plot_isotherms()
plt.gca().set(xlabel='$x_1$ / mole frac.', ylabel='$p$ / Pa')
plt.tight_layout(pad=0.2)
plt.savefig('pxdiagram.pdf')
plt.show()
def main(names):
T = 300
backend = 'HEOS'
tracer = vle.VLEIsolineTracer(IMPOSED_T, T, backend, names)
# Set flags
tracer.set_allowable_error(1e-6)
tracer.polishing(True)
tracer.set_debug_polishing(False)
tracer.set_forwards_integration(True)
tracer.set_unstable_termination(False)
tracer.set_stepping_variable(vle.VLEIsolineTracer.stepping_variable.STEP_IN_RHO1)
tracer.trace()
data = tracer.get_tracer_data()
print(tracer.get_termination_reason())
print(tracer.get_tracing_time())
fig, ax4 = plt.subplots(1,1,figsize=(6, 4))
x = np.array(data.x).T[0]
x1 = np.array(data.x).T[1]
y = np.array(data.y).T[0]
pL = np.array(data.pL)
rhoLmat = np.array(data.rhoL)
rhoVmat = np.array(data.rhoV)
# ax4.plot(x, rhoLmat[:,0],color='b')
# ax4.plot(x, rhoLmat[:,1],color='r')
# ax4.plot(x, rhoVmat[:,0],dashes = [2,2], color='b')
# ax4.plot(x, rhoVmat[:,1],dashes = [2,2], color='r')
dx0dpL = np.diff(x)/np.diff(pL)
dx1dpL = np.diff(x1)/np.diff(pL)
ax4.plot(x, pL, dashes=[2,2])
ax4.plot(y, pL, dashes=[2,2])
model = teqp.build_multifluid_model(names, '../mycp', '../mycp/dev/mixtures/mixture_binary_pairs.json')
# For now, using the EOS directly to get started... (ultimately either superancillary or extrapolation from critical)
AS = CP.AbstractState('HEOS', names[1])
AS.update(CP.QT_INPUTS, 0, T)
rhoL = AS.saturated_liquid_keyed_output(CP.iDmolar)
rhoV = AS.saturated_vapor_keyed_output(CP.iDmolar)
rhovec = np.array([0, rhoL, 0, rhoV])
# Initial slopes of dpdx1 and dpdy1
for m in [model]:
drhodp_liq, drhodp_vap = get_drhovecdp_sat(m, T, rhovec[0:2], rhovec[2::])
dpdxL = 1.0/get_dxdp(rhovec=rhovec[0:2], drhovecdp=drhodp_liq)
dpdxV = 1.0/get_dxdp(rhovec=rhovec[2::], drhovecdp=drhodp_vap)
xx = np.linspace(0, 0.4)
plt.plot(xx, xx*dpdxL[0]+AS.p(), dashes=[2, 2])
plt.plot(xx, xx*dpdxV[0]+AS.p(), dashes=[2, 2])
traceT(model, T, rhovec.copy())
# for f in np.linspace(0.9, 1.1, 100):
# modeltweak = teqp.build_BIPmodified(model, {'betaT': 1.0, 'gammaT': 1.001633952*f, 'betaV': 1.0, 'gammaV': 1.006268954})
# try:
# traceT(modeltweak, T, rhovec.copy(), kwargs=dict(lw=0.5, dashes=[3,1,1,1]))
# except:
# pass
for increment in np.linspace(-120,60,30):
T2 = T+increment
try:
AS = CP.AbstractState('HEOS', names[1])
AS.update(CP.QT_INPUTS, 0, T2)
rhoL = AS.saturated_liquid_keyed_output(CP.iDmolar)
rhoV = AS.saturated_vapor_keyed_output(CP.iDmolar)
rhovec = np.array([0, rhoL, 0, rhoV])
traceT_arclength(model, T2, rhovec.copy(), kwargs={'lw':0.5})
except BaseException as BE:
print(T2, BE)
tic = timeit.default_timer()
Tcvec = model.get_Tcvec()
rhocvec = 1/model.get_vcvec()
T0 = Tcvec[1]
rho0 = rhocvec*np.array([1,0])
curveJSON = teqp.trace_critical_arclength_binary(model, T0, rho0, "")
toc = timeit.default_timer()
print(toc-tic, 'to trace critical locus')
df = pandas.DataFrame(curveJSON)
df['z0 / mole frac.'] = df['rho0 / mol/m^3']/(df['rho0 / mol/m^3']+df['rho1 / mol/m^3'])
plt.plot(df['z0 / mole frac.'], df['p / Pa'], lw=2)
# plt.ylim(0.00019e6, 0.026e6)
if plot_critical:
plot_critical_curve()
plot_isotherms()
plt.gca().set(xlabel='$x_1$ / mole frac.', ylabel='$p$ / Pa')
plt.tight_layout(pad=0.2)
plt.show()
def trace_isopleth(names):
AS = CP.AbstractState('HEOS', '&'.join(names))
AS.set_mole_fractions([0.25, 0.75])
AS.build_phase_envelope("")
PE = AS.get_phase_envelope_data()
model = teqp.build_multifluid_model(names, '../mycp', '../mycp/dev/mixtures/mixture_binary_pairs.json')
i = 0
rhovec0 = np.array([
PE.rhomolar_liq[i]*PE.x[0][i], PE.rhomolar_liq[i]*PE.x[1][i],
PE.rhomolar_vap[i]*PE.y[0][i], PE.rhomolar_vap[i]*PE.y[1][i]])
T = PE.T[i]
dT = 1e-1
R = model.get_R(np.array([1.0]))
rhovec = rhovec0.copy()
def rhovecprime(T, rhovec):
return drhovecdT_isopleth(model=model, T=T, rhovec=rhovec)
tic = timeit.default_timer()
Tmax = 430
sol = solve_ivp(rhovecprime, [T, Tmax], y0=rhovec0.copy(), method='RK45',dense_output=True, t_eval = np.linspace(T,Tmax,1000))
p = [teqp.get_pr(model, T, sol.y[0:2,j]) + R*T*sol.y[0:2,j].sum() for T,j in zip(sol.t,range(len(sol.t)))]
x = sol.y[0,:]/(sol.y[0,:]+sol.y[1,:])
y = sol.y[2,:]/(sol.y[2,:]+sol.y[3,:])
toc = timeit.default_timer()
plt.plot(sol.t, p)
print(toc-tic)
# plt.plot(PE.T, PE.rhomolar_vap)
# plt.plot(PE.T, PE.rhomolar_liq)
plt.plot(PE.T, PE.p)
plt.show()
if ofname:
plt.savefig(ofname)
if show:
plt.show()
if __name__ == '__main__':
# main(names)
# trace_isopleth(names)
# names = ['Methane', 'n-Propane']
# model = teqp.build_multifluid_model(names, '../mycp', '../mycp/dev/mixtures/mixture_binary_pairs.json')
# puremodel = teqp.build_multifluid_model(['n-Propane'], '../mycp', '../mycp/dev/mixtures/mixture_binary_pairs.json')
# prettyPX(model, puremodel, np.linspace(180, 365, 100))
names = ['n-Butane', 'HydrogenSulfide']
model = teqp.build_multifluid_model(names, '../mycp', '../mycp/dev/mixtures/mixture_binary_pairs.json',{'estimate':'linear'})
names = ['Methane', 'n-Butane']
model = teqp.build_multifluid_model(names, '../mycp', '../mycp/dev/mixtures/mixture_binary_pairs.json')#,{'estimate':'Lorentz-Berthelot'})
puremodels = [teqp.build_multifluid_model([name], '../mycp', '../mycp/dev/mixtures/mixture_binary_pairs.json') for name in names]
prettyPX(model, puremodels, 0, np.linspace(230, 420, 100))
\ No newline at end of file
prettyPX(model, puremodels, names, 1, np.linspace(190, 423, 20), show=False, plot_critical=True, ofname='pxdiagram.pdf')
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment