From ebacd0eeaa186d26dc7e65b59bf75559c641ca86 Mon Sep 17 00:00:00 2001
From: Ian Bell <ian.bell@nist.gov>
Date: Fri, 21 May 2021 13:48:03 -0400
Subject: [PATCH] move the needed derivative into C++ and expose
---
include/teqp/algorithms/VLE.hpp | 77 +++++++++
include/teqp/derivs.hpp | 2 +-
interface/pybind11_wrapper.cpp | 6 +-
notebooks/isochoric.py | 277 ++++++++++++++++++++++++--------
4 files changed, 297 insertions(+), 65 deletions(-)
diff --git a/include/teqp/algorithms/VLE.hpp b/include/teqp/algorithms/VLE.hpp
index 5932ca7..f653ee3 100644
--- a/include/teqp/algorithms/VLE.hpp
+++ b/include/teqp/algorithms/VLE.hpp
@@ -130,4 +130,81 @@ auto extrapolate_from_critical(const Model& model, const Scalar Tc, const Scalar
auto rholiq = -drhohat/sqrt(1 - T/Tc) + rhoc;
auto rhovap = drhohat/sqrt(1 - T/Tc) + rhoc;
return (Eigen::ArrayXd(2) << rholiq, rhovap).finished();
+}
+
+template<class A, class B>
+auto linsolve(const A &a, const B& b) {
+ return a.matrix().colPivHouseholderQr().solve(b.matrix()).array().eval();
+}
+
+template<class Model, class Scalar, class VecType>
+auto get_drhovecdp_Tsat(const Model& model, const Scalar &T, const VecType& rhovecL, const VecType& rhovecV) {
+ //tic = timeit.default_timer();
+ using id = IsochoricDerivatives<Model, Scalar, VecType>;
+ auto Hliq = id::build_Psi_Hessian_autodiff(model, T, rhovecL).eval();
+ auto Hvap = id::build_Psi_Hessian_autodiff(model, T, rhovecV).eval();
+ //Hvap[~np.isfinite(Hvap)] = 1e20;
+ //Hliq[~np.isfinite(Hliq)] = 1e20;
+
+ auto N = rhovecL.size();
+ Eigen::MatrixXd A = decltype(Hliq)::Zero(N, N);
+ auto b = decltype(Hliq)::Ones(N, 1);
+ decltype(Hliq) drhodp_liq, drhodp_vap;
+ assert(rhovecL.size() == rhovecV.size());
+ if ((rhovecL != 0).all() && (rhovecV != 0).all()) {
+ // Normal treatment for all concentrations not equal to zero
+ A(0, 0) = Hliq.row(0).dot(rhovecV.matrix());
+ A(0, 1) = Hliq.row(1).dot(rhovecV.matrix());
+ A(1, 0) = Hliq.row(0).dot(rhovecL.matrix());
+ A(1, 1) = Hliq.row(1).dot(rhovecL.matrix());
+
+ drhodp_liq = linsolve(A, b);
+ drhodp_vap = linsolve(Hvap, Hliq*drhodp_liq);
+ }
+ else{
+ // Special treatment for infinite dilution
+ auto murL = id::build_Psir_gradient_autodiff(model, T, rhovecL);
+ auto murV = id::build_Psir_gradient_autodiff(model, T, rhovecV);
+ auto RL = model.R(rhovecL / rhovecL.sum());
+ auto RV = model.R(rhovecV / rhovecV.sum());
+
+ // First, for the liquid part
+ for (auto i = 0; i < N; ++i) {
+ for (auto j = 0; j < N; ++j) {
+ 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
+ auto Aij = (Hliq.row(j).array().cwiseProduct(((i == 0) ? rhovecV : rhovecL).array().transpose())).eval(); // coefficient - wise product
+ // A throwaway boolean for clarity
+ bool is_liq = (i == 1);
+ // Apply correction to the j term (RT if liquid, RT*phi for vapor)
+ Aij[j] = (is_liq) ? RL*T : RL*T*exp(-(murV[j] - murL[j])/(RL*T));
+ // Fill in entry
+ A(i, j) = Aij.sum();
+ }
+ else{
+ // Normal
+ A(i, j) = Hliq.row(j).dot(((i==0) ? rhovecV : rhovecL).matrix());
+ }
+ }
+ }
+ drhodp_liq = linsolve(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
+ auto diagrhovecL = rhovecL.matrix().asDiagonal();
+ auto PSIVstar = (diagrhovecL*Hvap).eval();
+ auto PSILstar = (diagrhovecL*Hliq).eval();
+ for (auto j = 0; j < N; ++j) {
+ if (rhovecL[j] == 0) {
+ PSILstar(j, j) = RL*T;
+ PSIVstar(j, j) = RV*T/exp(-(murV[j] - murL[j]) / (RV * T));
+ }
+ }
+ drhodp_vap = linsolve(PSIVstar, PSILstar*drhodp_liq);
+ }
+ return std::make_tuple(drhodp_liq, drhodp_vap);
}
\ No newline at end of file
diff --git a/include/teqp/derivs.hpp b/include/teqp/derivs.hpp
index 5568e0a..e1169bf 100644
--- a/include/teqp/derivs.hpp
+++ b/include/teqp/derivs.hpp
@@ -421,7 +421,7 @@ struct IsochoricDerivatives{
static auto build_Psi_Hessian_autodiff(const Model& model, const Scalar& T, const VectorType& rho) {
auto rhotot_ = rho.sum();
auto molefrac = (rho / rhotot_).eval();
- auto H = build_Psir_Hessian_autodiff(model, T, rho);
+ auto H = build_Psir_Hessian_autodiff(model, T, rho).eval();
for (auto i = 0; i < 2; ++i) {
H(i, i) += model.R(molefrac) * T / rho[i];
}
diff --git a/interface/pybind11_wrapper.cpp b/interface/pybind11_wrapper.cpp
index 1767bf1..117464b 100644
--- a/interface/pybind11_wrapper.cpp
+++ b/interface/pybind11_wrapper.cpp
@@ -23,6 +23,9 @@ namespace py = pybind11;
template<typename Model, typename Wrapper>
void add_derivatives(py::module &m, Wrapper &cls) {
+
+ using RAX = Eigen::Ref<Eigen::ArrayXd>;
+
using id = IsochoricDerivatives<Model, double, Eigen::Array<double,Eigen::Dynamic,1> >;
m.def("get_Ar00", &id::get_Ar00, py::arg("model"), py::arg("T"), py::arg("rho"));
m.def("get_Ar10", &id::get_Ar10, py::arg("model"), py::arg("T"), py::arg("rho"));
@@ -47,11 +50,12 @@ void add_derivatives(py::module &m, Wrapper &cls) {
m.def("extrapolate_from_critical", &extrapolate_from_critical<Model, double>);
m.def("pure_VLE_T", &pure_VLE_T<Model, double>);
+ m.def("get_drhovecdp_Tsat", &get_drhovecdp_Tsat<Model, double, RAX>, py::arg("model"), py::arg("T"), py::arg("rhovecL").noconvert(), py::arg("rhovecV").noconvert());
//cls.def("get_Ar01", [](const Model& m, const double T, const Eigen::ArrayXd& rhovec) { return id::get_Ar01(m, T, rhovec); });
//cls.def("get_Ar10", [](const Model& m, const double T, const Eigen::ArrayXd& rhovec) { return id::get_Ar10(m, T, rhovec); });
using tdx = TDXDerivatives<Model, double, Eigen::Array<double, Eigen::Dynamic, 1> >;
- using RAX = Eigen::Ref<Eigen::ArrayXd>;
+
cls.def("get_R", [](const Model& m, const RAX molefrac) { return m.R(molefrac); }, py::arg("molefrac").noconvert());
cls.def("get_Ar00", [](const Model& m, const double T, const double rho, const RAX molefrac) { return tdx::get_Ar00(m, T, rho, molefrac); }, py::arg("T"), py::arg("rho"), py::arg("molefrac").noconvert());
cls.def("get_Ar01", [](const Model& m, const double T, const double rho, const RAX molefrac) { return tdx::template get_Ar01<ADBackends::autodiff>(m, T, rho, molefrac); }, py::arg("T"), py::arg("rho"), py::arg("molefrac").noconvert());
diff --git a/notebooks/isochoric.py b/notebooks/isochoric.py
index e2542a6..3edea4c 100644
--- a/notebooks/isochoric.py
+++ b/notebooks/isochoric.py
@@ -14,6 +14,7 @@ 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
@@ -31,6 +32,7 @@ def get_drhovecdp_sat(model, T, rhovecL, rhovecV):
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
@@ -71,7 +73,8 @@ def get_drhovecdp_sat(model, T, rhovecL, rhovecV):
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):
@@ -96,36 +99,65 @@ def traceT(model, T, rhovec0, *, kwargs={}):
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)
- y = sol.y[2,:]/(sol.y[2,:]+sol.y[3,:])
+ y = sol.y[2,:]/(sol.y[2,:] + sol.y[3,:])
plt.plot(sol.t, p, 'r',**kwargs)
plt.plot(y, p, 'g',**kwargs)
-def drhovecdT_oldschool(model,T, rhovec, fluids):
- rhovecL = rhovec[0:2]
- rhovecV = rhovec[2::]
- tracer = vle.VLEIsolineTracer(vle.VLEIsolineTracer.imposed_variable.IMPOSED_T, 20000, 'HEOS', fluids)
- def get_derivs(T, rhovec):
- return tracer.get_derivs(T,rhovec)#vle.VLEIsolineTracer(vle.VLEIsolineTracer.imposed_variable.IMPOSED_T, 20000, 'HEOS', fluids).get_derivs(T, rhovec)
-
- derV = get_derivs(T, rhovecV)
- derL = get_derivs(T, rhovecL)
- # dT = 1e-3
- # print(derL.dpdT__constrhovec(), (get_derivs(T+dT, rhovecV).p() - get_derivs(T-dT, rhovecV).p())/(2*dT))
- AS = tracer.get_AbstractState_pointer()
- R = AS.gas_constant()
- # AS.set_mole_fractions(z)
- z = rhovecL/rhovecL.sum()
+def traceT_arclength(model, T, rhovec0, *, kwargs={}):
- rhoV = sum(rhovecV)
- rhoL = sum(rhovecL)
- DELTAs2 = np.array([derV.d2psir_dTdrhoi__constrhoj(i) - derL.d2psir_dTdrhoi__constrhoj(i) + R*np.log(rhovecV[i]/rhoV/(rhovecL[i]/rhoL)) for i in range(2)]) # Column vector
- DELTAs = np.array([derV.d2psir_dTdrhoi__constrhoj(i) - derL.d2psir_dTdrhoi__constrhoj(i) + R*np.log(rhovecV[i]/rhovecL[i]) for i in range(2)]) # Column vector
- DELTAbetarho = derV.dpdT__constrhovec() - derL.dpdT__constrhovec() # double
- PSIL = derL.get_Hessian()
- PSIV = derV.get_Hessian()
- drhovecL_dT = (np.dot(DELTAs,rhovecV)-DELTAbetarho)/np.dot(np.dot(PSIL,rhovecV-rhovecL),z)*np.array(z)
- drhovecV_dT = np.linalg.solve(PSIV, np.dot(PSIL,drhovecL_dT)-DELTAs)
- return np.array(drhovecL_dT.tolist()+drhovecV_dT.tolist())
+ R = model.get_R(np.array([0.0, 1.0]))
+ print(T)
+ c = 1.0
+ def rhovecprime(t, rhovec):
+ tic = timeit.default_timer()
+ drhovecdpL, drhovecdpV = teqp.get_drhovecdp_Tsat(model, T, rhovec[0:2], rhovec[2::])
+ 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])
+ return c*der
+
+ der_init = rhovecprime(0, rhovec0)
+ if np.any(rhovec0 + der_init*1e-6 < 0):
+ c *= -1
+ # print('init', der_init, T)
+
+ 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, 3500000], y0=rhovec0.copy(), method='RK45',dense_output=True, events=events)#, rtol=1e-8)
+ # print(sol.t_events, sol.y_events)
+ # print(sol.message)
+ 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)
+ 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)
@@ -137,7 +169,7 @@ def drhovecdT_isopleth(model, T, rhovec):
rhoV = sum(rhovecV)
xL = rhovecL/rhoL
xV = rhovecV/rhoV
- R = 8.31446161815324
+ 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))
@@ -147,6 +179,100 @@ def drhovecdT_isopleth(model, T, rhovec):
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)
+
+ 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):
+ """
+ ipure: int
+ index (0-based) of the pure fluid model to use
+ """
+ def plot_critical_curve():
+ 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.
+ T0 = Tcvec[k]
+ rho0 = rhocvec
+ rho0[1-k] = 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)
+
+ def plot_isotherms():
+ for T in Tvec:
+ puremodel = puremodels[ipure]
+ # rhoL, rhoV = robust_VLE(puremodel, T, puremodel.get_Tcvec()[0], 1/puremodel.get_vcvec()[0])
+ rhoL = CP.PropsSI('Dmolar','T',T,'Q',0,names[ipure])
+ rhoV = CP.PropsSI('Dmolar','T',T,'Q',1,names[ipure])
+ if ipure == 0:
+ # Both phases are pure of the first component with index 0
+ rhovec = np.array([rhoL, 0, rhoV, 0])
+ else:
+ # Both phases are pure of the second component with index 1
+ rhovec = np.array([0, rhoL, 0, rhoV])
+ try:
+ traceT_arclength(model, T, rhovec.copy(), kwargs={'lw': 0.5})
+ 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
@@ -154,8 +280,8 @@ def main(names):
tracer = vle.VLEIsolineTracer(IMPOSED_T, T, backend, names)
# Set flags
tracer.set_allowable_error(1e-6)
- tracer.polishing(False)
- tracer.set_debug_polishing(True)
+ tracer.polishing(True)
+ tracer.set_debug_polishing(False)
tracer.set_forwards_integration(True)
tracer.set_unstable_termination(False)
@@ -179,8 +305,8 @@ def main(names):
# 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)
- ax4.plot(y, pL, dashes=[2, 2])
+ 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')
@@ -202,14 +328,39 @@ def main(names):
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})
+ # 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:
- traceT(modeltweak, T, rhovec.copy(), kwargs=dict(lw=0.5, dashes=[3,1,1,1]))
- except:
- pass
+ 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)
- plt.ylim(0.00019e6, 0.026e6)
+ 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)
plt.gca().set(xlabel='$x_1$ / mole frac.', ylabel='$p$ / Pa')
plt.tight_layout(pad=0.2)
plt.show()
@@ -219,37 +370,28 @@ def trace_isopleth(names):
AS.set_mole_fractions([0.25, 0.75])
AS.build_phase_envelope("")
PE = AS.get_phase_envelope_data()
- print(dir(PE))
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[0]
- dT = 1e-2
+ T = PE.T[i]
+ dT = 1e-1
R = model.get_R(np.array([1.0]))
rhovec = rhovec0.copy()
- for i in range(10000):
- xL0 = rhovec[0]/rhovec[0:2].sum()
- molefrac = np.array([xL0, 1-xL0])
- rhovecL = rhovec[0:2]
- rhoL = rhovecL.sum()
- rhovecV = rhovec[2::]
- rhoV = rhovecV.sum()
-
- # plt.plot(T, rhoL, 'o')
- # plt.plot(T, rhoV, 'x')
-
- vec1 = drhovecdT_isopleth(model, T, rhovec)
- rhovec += vec1*dT
- T += dT
-
- rhovecL = rhovec[0:2]
- rhoL = rhovecL.sum()
- p = rhoL*R*T + teqp.get_pr(model, T, rhovecL)
- print(T, xL0, rhoL)
- plt.plot(T, p, 'x')
+
+ 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)
@@ -257,6 +399,15 @@ def trace_isopleth(names):
plt.show()
if __name__ == '__main__':
- names = ['n-Hexane', 'n-Octane']
# main(names)
- trace_isopleth(names)
\ No newline at end of file
+ # 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'})
+ 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
--
GitLab