diff --git a/tests/test_reaction.py b/tests/test_reaction.py index 52d8f39a..c8280f3c 100644 --- a/tests/test_reaction.py +++ b/tests/test_reaction.py @@ -9,6 +9,7 @@ """ import pytest import thermosteam as tmo +import numpy as np from thermosteam.reaction import ( Reaction, ParallelReaction, SeriesReaction, ReactionSystem, ReactionItem ) @@ -213,8 +214,10 @@ def test_reaction_enthalpy_with_phases(): reactant='Glucose', X=1, correct_atomic_balance=True) assert_allclose(combustion.dH, -2787650.3239119546, rtol=1e-3) - -def test_reactive_phase_equilibrium(): + +def test_reactive_phase_equilibrium_no_kinetics(): + import thermosteam as tmo + from numpy.testing import assert_allclose tmo.settings.set_thermo(['EthylLactate', 'LacticAcid', 'H2O', 'Ethanol'], cache=True) rxn = tmo.Reaction('LacticAcid + Ethanol -> H2O + EthylLactate', reactant='LacticAcid', X=0.2) stream = tmo.Stream( @@ -236,6 +239,153 @@ def test_reactive_phase_equilibrium(): 1.8234109156732896] ) + stream.vle(T=340, P=101325, liquid_reaction=rxn) + assert_allclose( + stream.imol['l'], + stream.mol + ) + assert_allclose( + stream.imol['g'], + 0 + ) + stream.vle(T=450, P=101325, liquid_reaction=rxn) + assert_allclose( + stream.imol['g'], + stream.mol + ) + assert_allclose( + stream.imol['l'], + 0 + ) + +def test_reactive_phase_equilibrium_with_kinetics(): + import thermosteam as tmo + from math import exp + from numpy.testing import assert_allclose + tmo.settings.set_thermo(['EthylLactate', 'LacticAcid', 'H2O', 'Ethanol'], cache=True) + + class Esterification(tmo.Reaction): + __slots__ = ('mcat',) + + def __init__(self): + super().__init__( + 'LacticAcid + Ethanol -> H2O + EthylLactate', + reactant='LacticAcid', X=0.2 + ) + self.mcat = 0.01 + + def _rate(self, stream): + R = tmo.constants.R + T = stream.T + kf = 6.52e3 * exp(-4.8e4 / (R * T)) + kr = 2.72e3 * exp(-4.8e4 / (R * T)) + LaEt, La, H2O, EtOH = stream.mol / stream.F_mol + return self.mcat * (kf * La * EtOH - kr * LaEt * H2O) + + rxn = Esterification() + stream = tmo.Stream( + H2O=2, Ethanol=5, LacticAcid=1, T=355, + ) + rxn(stream) + assert_allclose( + stream.mol, + [0.7998037291462218, 0.20019627085377822, + 2.799803729146222, 4.200196270853779] + ) + stream = tmo.Stream( + H2O=2, Ethanol=5, LacticAcid=1, T=355, + ) + T = 360 + P = 101325 + stream.vle(T=T, P=P, liquid_reaction=rxn) + assert_allclose( + stream.imol['l'], + [0.35781246995914223, + 0.28312078337665264, + 1.291316286965713, + 1.3287108564655643], + rtol=1e-6, + atol=1e-6, + ) + assert_allclose( + stream.imol['g'], + [0.0, + 0.35906674666420507, + 1.0664961829934292, + 3.313476673575294], + rtol=1e-6, + atol=1e-6, + ) + V = stream.vapor_fraction + H = stream.H + stream.Hf + stream = tmo.Stream( + H2O=2, Ethanol=5, LacticAcid=1, T=355, + ) + stream.vle(V=V, P=P, liquid_reaction=rxn) + assert_allclose( + stream.imol['l'], + [0.35781246995914223, + 0.28312078337665264, + 1.291316286965713, + 1.3287108564655643], + rtol=1e-6, + atol=1e-6, + ) + assert_allclose( + stream.imol['g'], + [0.0, + 0.35906674666420507, + 1.0664961829934292, + 3.313476673575294], + rtol=1e-6, + atol=1e-6, + ) + stream = tmo.Stream( + H2O=2, Ethanol=5, LacticAcid=1, T=355, + ) + stream.vle(V=V, T=T, liquid_reaction=rxn) + assert_allclose( + stream.imol['l'], + [0.35781246995914223, + 0.28312078337665264, + 1.291316286965713, + 1.3287108564655643], + rtol=1e-6, + atol=1e-6, + ) + assert_allclose( + stream.imol['g'], + [0.0, + 0.35906674666420507, + 1.0664961829934292, + 3.313476673575294], + rtol=1e-6, + atol=1e-6, + ) + stream = tmo.Stream( + H2O=2, Ethanol=5, LacticAcid=1, T=355, + ) + stream.vle(H=H, P=P, liquid_reaction=rxn) + assert_allclose( + stream.imol['l'], + [0.35781246995914223, + 0.28312078337665264, + 1.291316286965713, + 1.3287108564655643], + rtol=1e-6, + atol=1e-6, + ) + assert_allclose( + stream.imol['g'], + [0.0, + 0.35906674666420507, + 1.0664961829934292, + 3.313476673575294], + rtol=1e-6, + atol=1e-6, + ) + + def test_repr(): cal2joule = 4.184 Glucan = tmo.Chemical('Glucan', search_db=False, formula='C6H10O5', Hf=-233200*cal2joule, phase='s', default=True) @@ -264,5 +414,6 @@ def test_repr(): test_reaction() test_reaction_enthalpy_balance() test_reaction_enthalpy_with_phases() - test_reactive_phase_equilibrium() + test_reactive_phase_equilibrium_no_kinetics() + test_reactive_phase_equilibrium_with_kinetics() test_repr() \ No newline at end of file diff --git a/thermosteam/_multi_stream.py b/thermosteam/_multi_stream.py index 9dbf76c8..8263a5fe 100644 --- a/thermosteam/_multi_stream.py +++ b/thermosteam/_multi_stream.py @@ -361,9 +361,7 @@ def reset_cache(self): self._streams = {} self._vle_cache = eq.VLECache(self._imol, self._thermal_condition, - self._thermo, - self._bubble_point_cache, - self._dew_point_cache) + self._thermo) self._lle_cache = eq.LLECache(self._imol, self._thermal_condition, self._thermo) @@ -381,8 +379,6 @@ def __getitem__(self, phase): stream._imol = self._imol.get_phase(phase) stream._thermal_condition = self._thermal_condition stream._thermo = self._thermo - stream._bubble_point_cache = self._bubble_point_cache - stream._dew_point_cache = self._dew_point_cache stream._property_cache = {} stream.characterization_factors = {} stream._property_cache_key = None, None diff --git a/thermosteam/_stream.py b/thermosteam/_stream.py index 9d2173c0..17c93b3a 100644 --- a/thermosteam/_stream.py +++ b/thermosteam/_stream.py @@ -27,7 +27,6 @@ from .base import SparseVector, SparseArray from numpy.typing import NDArray from typing import Optional, Sequence, Callable - import biosteam as bst # from .constants import g MaterialIndexer = tmo.indexer.MaterialIndexer @@ -274,7 +273,6 @@ class Stream(AbstractStream): """ __slots__ = ( '_imol', '_thermal_condition', '_streams', - '_bubble_point_cache', '_dew_point_cache', '_vle_cache', '_lle_cache', '_sle_cache', '_price', '_property_cache_key', '_property_cache', 'characterization_factors', 'equations', @@ -767,8 +765,6 @@ def _init_indexer(self, flow, phase, chemicals, chemical_flows): def reset_cache(self): """Reset cache regarding equilibrium methods.""" - self._bubble_point_cache = eq.BubblePointCache() - self._dew_point_cache = eq.DewPointCache() self._property_cache_key = None, None self._property_cache = {} @@ -1965,8 +1961,6 @@ def proxy(self, ID=None): new._thermal_condition = self._thermal_condition new._property_cache = self._property_cache new._property_cache_key = self._property_cache_key - new._bubble_point_cache = self._bubble_point_cache - new._dew_point_cache = self._dew_point_cache new.equations = self.equations new.characterization_factors = self.characterization_factors return new @@ -2022,9 +2016,7 @@ def vlle(self, T, P): imol = self.imol vle = eq.VLE(imol, self._thermal_condition, - self._thermo, - self._bubble_point_cache, - self._dew_point_cache) + self._thermo) lle = eq.LLE(imol, self._thermal_condition, self._thermo) @@ -2094,9 +2086,7 @@ def get_bubble_point(self, IDs: Optional[Sequence[str]]=None): BubblePoint([Water, Ethanol]) """ - chemicals = self.chemicals[IDs] if IDs else self.vle_chemicals - bp = self._bubble_point_cache(chemicals, self._thermo) - return bp + return eq.BubblePoint(self.chemicals[IDs] if IDs else self.vle_chemicals, self._thermo) def get_dew_point(self, IDs: Optional[Sequence[str]]=None): """ @@ -2116,9 +2106,7 @@ def get_dew_point(self, IDs: Optional[Sequence[str]]=None): DewPoint([Water, Ethanol]) """ - chemicals = self.chemicals[IDs] if IDs else self.vle_chemicals - dp = self._dew_point_cache(chemicals, self._thermo) - return dp + return eq.DewPoint(self.chemicals[IDs] if IDs else self.vle_chemicals, self._thermo) def bubble_point_at_T(self, T: Optional[float]=None, IDs: Optional[Sequence[str]]=None): """ @@ -2502,9 +2490,7 @@ def phases(self, phases): self._streams = {} self._vle_cache = eq.VLECache(self._imol, self._thermal_condition, - self._thermo, - self._bubble_point_cache, - self._dew_point_cache) + self._thermo) self._lle_cache = eq.LLECache(self._imol, self._thermal_condition, self._thermo) diff --git a/thermosteam/equilibrium/bubble_point.py b/thermosteam/equilibrium/bubble_point.py index 174a2456..374b065c 100644 --- a/thermosteam/equilibrium/bubble_point.py +++ b/thermosteam/equilibrium/bubble_point.py @@ -13,11 +13,10 @@ from .domain import vle_domain from ..exceptions import InfeasibleRegion from .. import functional as fn -from ..utils import fill_like, Cache from .._settings import settings __all__ = ( - 'BubblePoint', 'BubblePointValues', 'BubblePointCache', + 'BubblePoint', 'BubblePointValues', 'BubblePointBeta' ) @@ -54,6 +53,22 @@ def __repr__(self): return f"{type(self).__name__}(T={self.T:.2f}, P={self.P:.0f}, IDs={self.IDs}, z={self.z}, y={self.y})" +class ReactiveBubblePointValues: + __slots__ = ('T', 'P', 'IDs', 'z0', 'dz', 'y', 'x') + + def __init__(self, T, P, IDs, z0, dz, y, x): + self.T = T + self.P = P + self.IDs = IDs + self.z0 = z0 + self.dz = dz + self.y = y + self.x = x + + def __repr__(self): + return f"{type(self).__name__}(T={self.T:.2f}, P={self.P:.0f}, IDs={self.IDs}, z0={self.z0}, dz={self.dz}, y={self.y}, x={self.x})" + + # %% Bubble point calculation class BubblePoint: @@ -90,15 +105,15 @@ class BubblePoint: _cached = {} T_tol = 1e-9 P_tol = 1e-3 - def __init__(self, chemicals=(), thermo=None): + def __new__(cls, chemicals=(), thermo=None): thermo = settings.get_default_thermo(thermo) chemicals = tuple(chemicals) key = (chemicals, thermo.Gamma, thermo.Phi, thermo.PCF) - cached = self._cached + cached = cls._cached if key in cached: - other = cached[key] - fill_like(self, other, self.__slots__) + return cached[key] else: + self = super().__new__(cls) self.IDs = tuple([i.ID for i in chemicals]) self.gamma = thermo.Gamma(chemicals) self.phi = thermo.Phi(chemicals) @@ -111,6 +126,7 @@ def __init__(self, chemicals=(), thermo=None): self.Pmax = max([i(Tmax) for i in Psats]) self.chemicals = chemicals cached[key] = self + return self def _T_error(self, T, P, z_over_P, z_norm, y): if T <= 0: raise InfeasibleRegion('negative temperature') @@ -128,6 +144,29 @@ def _P_error(self, P, T, z_Psat_gamma, Psats, y): y[:] = solve_y(y_phi, self.phi, T, P, y) return 1. - y.sum() + def _T_error_reactive(self, T, P, z, dz, y, x, liquid_reaction): + if T <= 0: raise InfeasibleRegion('negative temperature') + dz[:] = liquid_reaction.conversion(z, T, P, 'l') + x[:] = z + dz + x /= x.sum() + Psats = np.array([i(T) for i in self.Psats], dtype=float) + y_phi = (x / P + * Psats + * self.gamma(x, T) + * self.pcf(T, P, Psats)) + y[:] = solve_y(y_phi, self.phi, T, P, y) + return 1. - y.sum() + + def _P_error_reactive(self, P, T, Psats, z, dz, y, x, liquid_reaction): + if P <= 0: raise InfeasibleRegion('negative pressure') + dz[:] = liquid_reaction.conversion(z, T, P, 'l') + x[:] = z + dz + x /= x.sum() + z_Psat_gamma = x * Psats * self.gamma(x, T) + y_phi = z_Psat_gamma * self.pcf(T, P, Psats) / P + y[:] = solve_y(y_phi, self.phi, T, P, y) + return 1. - y.sum() + def _T_error_ideal(self, T, z_over_P, y): y[:] = z_over_P * np.array([i(T) for i in self.Psats], dtype=float) return 1 - y.sum() @@ -153,18 +192,21 @@ def _Py_ideal(self, z_Psat_gamma_pcf): y = z_Psat_gamma_pcf / P return P, y - def __call__(self, z, *, T=None, P=None): + def __call__(self, z, *, T=None, P=None, liquid_reaction=None): z = np.asarray(z, float) if T: if P: raise ValueError("may specify either T or P, not both") - P, y = self.solve_Py(z, T) + P, *args = self.solve_Py(z, T, liquid_reaction) elif P: - T, y = self.solve_Ty(z, P) + T, *args = self.solve_Ty(z, P, liquid_reaction) else: raise ValueError("must specify either T or P") - return BubblePointValues(T, P, self.IDs, z, y) + if liquid_reaction: + return ReactiveBubblePointValues(T, P, self.IDs, z, *args) + else: + return BubblePointValues(T, P, self.IDs, z, *args) - def solve_Ty(self, z, P): + def solve_Ty(self, z, P, liquid_reaction=None): """ Bubble point at given composition and pressure. @@ -201,7 +243,8 @@ def solve_Ty(self, z, P): chemical = self.chemicals[fn.first_true_index(positives)] T = chemical.Tsat(P, check_validity=False) if P <= chemical.Pc else chemical.Tc y = z.copy() - else: + return T, fn.normalize(y) + elif liquid_reaction is None: f = self._T_error z_norm = z / z.sum() z_over_P = z/P @@ -217,9 +260,28 @@ def solve_Ty(self, z, P): f(Tmin, *args), f(Tmax, *args), T_guess, self.T_tol, 5e-12, args, checkiter=False, checkbounds=False) - return T, fn.normalize(y) + return T, fn.normalize(y) + else: + f = self._T_error_reactive + z_norm = z / z.sum() + x = z_norm.copy() + dz = z_norm.copy() + z_over_P = z / P + T_guess, y = self._Ty_ideal(z_over_P) + args = (P, z_norm, dz, y, x, liquid_reaction) + try: + T = flx.aitken_secant(f, T_guess, T_guess + 1e-3, + self.T_tol, 5e-12, args, + checkiter=False) + except RuntimeError: + Tmin = self.Tmin; Tmax = self.Tmax + T = flx.IQ_interpolation(f, Tmin, Tmax, + f(Tmin, *args), f(Tmax, *args), + T_guess, self.T_tol, 5e-12, args, + checkiter=False, checkbounds=False) + return T, dz, fn.normalize(y), x - def solve_Py(self, z, T): + def solve_Py(self, z, T, liquid_reaction=None): """ Bubble point at given composition and temperature. @@ -256,15 +318,17 @@ def solve_Py(self, z, T): chemical = self.chemicals[fn.first_true_index(positives)] P = chemical.Psat(T) if T <= chemical.Tc else chemical.Pc y = z.copy() - else: + return P, fn.normalize(y) + elif liquid_reaction is None: if T > self.Tmax: T = self.Tmax elif T < self.Tmin: T = self.Tmin - Psat = np.array([i(T) for i in self.Psats]) + Psats = np.array([i(T) for i in self.Psats]) z_norm = z / z.sum() - z_Psat_gamma = z * Psat * self.gamma(z_norm, T) + z_Psat_gamma = z_norm * Psats * self.gamma(z_norm, T) f = self._P_error P_guess, y = self._Py_ideal(z_Psat_gamma) - args = (T, z_Psat_gamma, Psat, y) + args = (T, z_Psat_gamma, Psats, y) + breakpoint() try: P = flx.aitken_secant(f, P_guess, P_guess-1, self.P_tol, 1e-9, args, checkiter=False) @@ -274,7 +338,26 @@ def solve_Py(self, z, T): f(Pmin, *args), f(Pmax, *args), P_guess, self.P_tol, 5e-12, args, checkiter=False, checkbounds=False) - return P, fn.normalize(y) + return P, fn.normalize(y) + else: + f = self._P_error_reactive + z_norm = z / z.sum() + Psats = np.array([i(T) for i in self.Psats]) + x = z_norm.copy() + dz = z_norm.copy() + z_Psat_gamma = z * Psats * self.gamma(z_norm, T) + P_guess, y = self._Py_ideal(z_Psat_gamma) + args = (T, Psats, z_norm, dz, y, x, liquid_reaction) + try: + P = flx.aitken_secant(f, P_guess, P_guess-1, self.P_tol, 1e-9, + args, checkiter=False) + except RuntimeError: + Pmin = self.Pmin; Pmax = self.Pmax + P = flx.IQ_interpolation(f, Pmin, Pmax, + f(Pmin, *args), f(Pmax, *args), + P_guess, self.P_tol, 5e-12, args, + checkiter=False, checkbounds=False) + return P, dz, fn.normalize(y), x def __repr__(self): chemicals = ", ".join([i.ID for i in self.chemicals]) @@ -318,7 +401,7 @@ def __init__(self, chemicals=(), flasher=None): __call__ = BubblePoint.__call__ - def solve_Ty(self, z, P): + def solve_Ty(self, z, P, liquid_reaction=None): """ Bubble point at given composition and pressure. @@ -360,7 +443,7 @@ def solve_Ty(self, z, P): T = results.T return T, fn.normalize(y) - def solve_Py(self, z, T): + def solve_Py(self, z, T, liquid_reaction=None): """ Bubble point at given composition and temperature. @@ -406,6 +489,3 @@ def __repr__(self): chemicals = ", ".join([i.ID for i in self.chemicals]) return f"{type(self).__name__}([{chemicals}])" -class BubblePointCache(Cache): load = BubblePoint -del Cache - diff --git a/thermosteam/equilibrium/dew_point.py b/thermosteam/equilibrium/dew_point.py index ca5960ac..5b9a20c5 100644 --- a/thermosteam/equilibrium/dew_point.py +++ b/thermosteam/equilibrium/dew_point.py @@ -12,11 +12,10 @@ from numba import njit from .. import functional as fn from ..exceptions import InfeasibleRegion -from ..utils import fill_like, Cache from .._settings import settings from .domain import vle_domain -__all__ = ('DewPoint', 'DewPointCache') +__all__ = ('DewPoint',) # %% Solvers @@ -62,6 +61,22 @@ def __repr__(self): return f"{type(self).__name__}(T={self.T:.2f}, P={self.P:.0f}, IDs={self.IDs}, z={self.z}, x={self.x})" +class ReactiveDewPointValues: + __slots__ = ('T', 'P', 'IDs', 'z0', 'dz', 'y', 'x') + + def __init__(self, T, P, IDs, z0, dz, y, x): + self.T = T + self.P = P + self.IDs = IDs + self.z0 = z0 + self.dz = dz + self.y = y + self.x = x + + def __repr__(self): + return f"{type(self).__name__}(T={self.T:.2f}, P={self.P:.0f}, IDs={self.IDs}, z0={self.z0}, dz={self.dz}, y={self.y}, x={self.x})" + + # %% Dew point calculation class DewPoint: @@ -99,15 +114,15 @@ class DewPoint: _cached = {} T_tol = 1e-9 P_tol = 1e-3 - def __init__(self, chemicals=(), thermo=None): + def __new__(cls, chemicals=(), thermo=None): thermo = settings.get_default_thermo(thermo) chemicals = tuple(chemicals) key = (chemicals, thermo.Gamma, thermo.Phi, thermo.PCF) - cached = self._cached + cached = cls._cached if key in cached: - other = cached[key] - fill_like(self, other, self.__slots__) + return cached[key] else: + self = super().__new__(cls) self.IDs = tuple([i.ID for i in chemicals]) self.gamma = thermo.Gamma(chemicals) self.phi = thermo.Phi(chemicals) @@ -120,6 +135,7 @@ def __init__(self, chemicals=(), thermo=None): self.Pmax = max([i(Tmax) for i in Psats]) self.chemicals = chemicals cached[key] = self + return self def _solve_x(self, x_gamma, T, P, x): gamma = self.gamma @@ -135,6 +151,19 @@ def _T_error(self, T, P, z_norm, zP, x): x[:] = self._solve_x(x_gamma, T, P, x) return 1 - x.sum() + def _T_error_reactive(self, T, P, z, dz, y, x, gas_reaction): + if T <= 0: raise InfeasibleRegion('negative temperature') + dz[:] = gas_reaction.conversion(z, T, P, 'g') + y[:] = z + dz + y /= y.sum() + Psats = np.array([i(T) for i in self.Psats]) + Psats[Psats < 1e-16] = 1e-16 # Prevent floating point error + phi = self.phi(y, T, P) + pcf = self.pcf(T, P, Psats) + x_gamma = phi * y * P / Psats / pcf + x[:] = self._solve_x(x_gamma, T, P, x) + return 1 - x.sum() + def _T_error_ideal(self, T, zP, x): Psats = np.array([i(T) for i in self.Psats]) Psats[Psats < 1e-16] = 1e-16 # Prevent floating point error @@ -147,6 +176,15 @@ def _P_error(self, P, T, z_norm, z_over_Psats, Psats, x): x[:] = self._solve_x(x_gamma, T, P, x) return 1 - x.sum() + def _P_error_reactive(self, P, T, Psats, z, dz, y, x, gas_reaction): + if P <= 0: raise InfeasibleRegion('negative pressure') + dz[:] = gas_reaction.conversion(z, T, P, 'g') + y[:] = z + dz + y /= y.sum() + x_gamma = y / Psats * P * self.phi(y, T, P) / self.pcf(T, P, Psats) + x[:] = self._solve_x(x_gamma, T, P, x) + return 1 - x.sum() + def _Tx_ideal(self, zP): f = self._T_error_ideal Tmin = self.Tmin + 10. @@ -167,18 +205,21 @@ def _Px_ideal(self, z_over_Psats): x = z_over_Psats * P return P, x - def __call__(self, z, *, T=None, P=None): + def __call__(self, z, *, T=None, P=None, gas_reaction=None): z = np.asarray(z, float) if T: if P: raise ValueError("may specify either T or P, not both") - P, x = self.solve_Px(z, T) + P, *args = self.solve_Px(z, T, gas_reaction) elif P: - T, x = self.solve_Tx(z, P) + T, *args = self.solve_Tx(z, P, gas_reaction) else: raise ValueError("must specify either T or P") - return DewPointValues(T, P, self.IDs, z, x) + if gas_reaction: + return ReactiveDewPointValues(T, P, self.IDs, z, *args) + else: + return DewPointValues(T, P, self.IDs, z, *args) - def solve_Tx(self, z, P): + def solve_Tx(self, z, P, gas_reaction=None): """ Dew point given composition and pressure. @@ -215,7 +256,8 @@ def solve_Tx(self, z, P): chemical = self.chemicals[fn.first_true_index(positives)] T = chemical.Tsat(P, check_validity=False) if P <= chemical.Pc else chemical.Tc x = z.copy() - else: + return T, fn.normalize(x) + elif gas_reaction is None: f = self._T_error z_norm = z/z.sum() zP = z * P @@ -230,11 +272,30 @@ def solve_Tx(self, z, P): Tmax = self.Tmax T = flx.IQ_interpolation(f, Tmin, Tmax, f(Tmin, *args), f(Tmax, *args), - T_guess, 1e-9, 5e-12, args, + T_guess, self.T_tol, 5e-12, args, + checkiter=False, checkbounds=False) + return T, fn.normalize(x) + else: + f = self._T_error_reactive + z_norm = z / z.sum() + x = z_norm.copy() + dz = z_norm.copy() + zP = z * P + T_guess, y = self._Tx_ideal(zP) + args = (P, z_norm, dz, y, x, gas_reaction) + try: + T = flx.aitken_secant(f, T_guess, T_guess + 1e-3, + self.T_tol, 5e-12, args, + checkiter=False) + except RuntimeError: + Tmin = self.Tmin; Tmax = self.Tmax + T = flx.IQ_interpolation(f, Tmin, Tmax, + f(Tmin, *args), f(Tmax, *args), + T_guess, self.T_tol, 5e-12, args, checkiter=False, checkbounds=False) - return T, fn.normalize(x) + return T, dz, fn.normalize(y), x - def solve_Px(self, z, T): + def solve_Px(self, z, T, gas_reaction=None): """ Dew point given composition and temperature. @@ -271,10 +332,11 @@ def solve_Px(self, z, T): chemical = self.chemicals[fn.first_true_index(z)] P = chemical.Psat(T) if T <= chemical.Tc else chemical.Pc x = z.copy() - else: - z_norm = z/z.sum() + return P, fn.normalize(x) + elif gas_reaction is None: + z_norm = z / z.sum() Psats = np.array([i(T) for i in self.Psats], dtype=float) - z_over_Psats = z/Psats + z_over_Psats = z / Psats P_guess, x = self._Px_ideal(z_over_Psats) args = (T, z_norm, z_over_Psats, Psats, x) f = self._P_error @@ -288,11 +350,28 @@ def solve_Px(self, z, T): f(Pmin, *args), f(Pmax, *args), P_guess, self.P_tol, 5e-12, args, checkiter=False, checkbounds=False) - return P, fn.normalize(x) + return P, fn.normalize(x) + else: + f = self._P_error_reactive + z_norm = z / z.sum() + y = z_norm.copy() + dz = z_norm.copy() + Psats = np.array([i(T) for i in self.Psats], dtype=float) + z_over_Psats = z / Psats + P_guess, x = self._Px_ideal(z_over_Psats) + args = (T, Psats, z_norm, dz, y, x, gas_reaction) + try: + P = flx.aitken_secant(f, P_guess, P_guess-1, self.P_tol, 1e-9, + args, checkiter=False) + except RuntimeError: + Pmin = self.Pmin; Pmax = self.Pmax + P = flx.IQ_interpolation(f, Pmin, Pmax, + f(Pmin, *args), f(Pmax, *args), + P_guess, self.P_tol, 5e-12, args, + checkiter=False, checkbounds=False) + return P, dz, fn.normalize(y), x def __repr__(self): chemicals = ", ".join([i.ID for i in self.chemicals]) return f"{type(self).__name__}([{chemicals}])" -class DewPointCache(Cache): load = DewPoint -del Cache diff --git a/thermosteam/equilibrium/vle.py b/thermosteam/equilibrium/vle.py index fce734a7..83c86688 100644 --- a/thermosteam/equilibrium/vle.py +++ b/thermosteam/equilibrium/vle.py @@ -15,8 +15,8 @@ from ..exceptions import InfeasibleRegion, NoEquilibrium from . import binary_phase_fraction as binary from .equilibrium import Equilibrium -from .dew_point import DewPointCache -from .bubble_point import BubblePointCache +from .dew_point import DewPoint +from .bubble_point import BubblePoint from .fugacity_coefficients import IdealFugacityCoefficients from .poyinting_correction_factors import MockPoyintingCorrectionFactors from . import activity_coefficients as ac @@ -40,34 +40,23 @@ def update_xV(xV, V, Ks, z): xV[-1] = V xV[:-1] = z/(1. + V * (Ks - 1.)) -def xV_iter(xVlogK, pcf_Psat_over_P, T, P, - z, z_light, z_heavy, f_gamma, gamma_args, - f_phi, n): - xVlogK = xVlogK.copy() - xV = xVlogK[:-n] - Ks = np.exp(xVlogK[-n:]) - try: - x, y = xy(xV, Ks) - except: - breakpoint() - Ks[:] = pcf_Psat_over_P * f_gamma(x, T, *gamma_args) / f_phi(y, T, P) - V = binary.solve_phase_fraction_Rashford_Rice(z, Ks, xV[-1], z_light, z_heavy) - update_xV(xV, V, Ks, z) - xVlogK[-n:] = np.log(Ks) - return xVlogK - -def xV_iter_2n(xVlogK, pcf_Psat_over_P, T, P, z, f_gamma, gamma_args, f_phi, n): +def xV_iter_2n(xVlogK, pcf_Psat_over_P, T, P, z, f_gamma, gamma_args, f_phi, + n, gas_reaction, liquid_reaction, index): xVlogK = xVlogK.copy() xV = xVlogK[:-n] Ks = np.exp(xVlogK[-n:]) x, y = xy(xV, Ks) Ks[:] = pcf_Psat_over_P * f_gamma(x, T, *gamma_args) / f_phi(y, T, P) + if gas_reaction: + z = z + gas_reaction.conversion(material=y * xV[n], T=T, P=P, phase='g') + if liquid_reaction: + z = z + liquid_reaction.conversion(material=x * (1 - xV[n]), T=T, P=P, phase='l') V = binary.compute_phase_fraction_2N(z, Ks) update_xV(xV, V, Ks, z) xVlogK[-n:] = np.log(Ks) return xVlogK -def xV_reactive_iter( +def xV_iter( xVlogK, pcf_Psat_over_P, T, P, z, z_light, z_heavy, f_gamma, gamma_args, f_phi, n, gas_reaction, liquid_reaction, index @@ -78,14 +67,9 @@ def xV_reactive_iter( Ks = np.exp(xVlogK[-n:]) x, y = xy(xV, Ks) if gas_reaction: - material = SparseVector.from_size(gas_reaction.chemicals.size) - material[index] = y * V - z = z + gas_reaction.conversion(material=material, T=T, P=P, phase='g')[index] + z = z + gas_reaction.conversion(material=y * V, T=T, P=P, phase='g') if liquid_reaction: - material = SparseVector.from_size(liquid_reaction.chemicals.size) - material[index] = x * (1 - V) - # breakpoint() - z = z + liquid_reaction.conversion(material=material, T=T, P=P, phase='l')[index] + z = z + liquid_reaction.conversion(material=x * (1 - V), T=T, P=P, phase='l') Ks[:] = pcf_Psat_over_P * f_gamma(x, T, *gamma_args) / f_phi(y, T, P) V = binary.solve_phase_fraction_Rashford_Rice(z, Ks, xV[-1], z_light, z_heavy) update_xV(xV, V, Ks, z) @@ -109,10 +93,6 @@ class VLE(Equilibrium, phases='lg'): thermo=None : :class:`~thermosteam.Thermo`, optional Themodynamic property package for equilibrium calculations. Defaults to `thermosteam.settings.get_thermo()`. - bubble_point_cache=None : :class:`~thermosteam.utils.Cache`, optional - Cache to retrieve bubble point object. - dew_point_cache=None : :class:`~thermosteam.utils.Cache`, optional - Cache to retrieve dew point object Examples -------- @@ -299,8 +279,6 @@ class VLE(Equilibrium, phases='lg'): '_F_mol_vle', # [float] Total moles in equilibrium. '_F_mol_light', # [float] Total moles of gas chemicals not included in equilibrium calculation. '_F_mol_heavy', # [float] Total moles of heavy chemicals not included in equilibrium calculation. - '_dew_point_cache', # [Cache] Retrieves the DewPoint object if arguments are the same. - '_bubble_point_cache', # [Cache] Retrieves the BubblePoint object if arguments are the same. '_dmol_vle', '_dF_mol', ) @@ -315,62 +293,15 @@ class VLE(Equilibrium, phases='lg'): default_method = 'fixed-point' def __init__(self, imol=None, thermal_condition=None, - thermo=None, bubble_point_cache=None, dew_point_cache=None): + thermo=None): self.method = self.default_method self._T = self._P = self._H_hat = self._V = 0 - self._dew_point_cache = dew_point_cache or DewPointCache() - self._bubble_point_cache = bubble_point_cache or BubblePointCache() super().__init__(imol, thermal_condition, thermo) self._x = None self._z_last = None self._nonzero = None self._index = () - def reactive_vle(self, *, T=None, P=None, V=None, H=None, S=None, - gas_reaction=None, liquid_reaction=None): - ### Decide what kind of equilibrium to run ### - T_spec = T is not None - P_spec = P is not None - V_spec = V is not None - H_spec = H is not None - S_spec = S is not None - N_specs = (T_spec + P_spec + V_spec + H_spec + S_spec) - assert N_specs == 2, ("must pass two and only two of the following " - "specifications: T, P, V, H, S") - - # Run equilibrium - imol = self._imol - if T_spec: - if P_spec: - self.set_thermal_condition(T, P) - self.set_thermal_condition_reactive(T, P, gas_reaction, liquid_reaction) - elif V_spec: - self.set_TV(T, V) - raise NotImplementedError - elif H_spec: - self.set_TH(T, H) - raise NotImplementedError - else: - self.set_TS(T, S) - raise NotImplementedError - elif P_spec: - if V_spec: - self.set_PV(P, V) - raise NotImplementedError - elif H_spec: - self.set_PH(P, H, stacklevel=1) - raise NotImplementedError - else: - self.set_PS(P, S, stacklevel=1) - raise NotImplementedError - elif S_spec: # pragma: no cover - if H_spec: - raise NotImplementedError('specification H and S is invalid') - else: # V_spec - raise NotImplementedError('specification V and S not implemented') - elif H_spec: # pragma: no cover - raise NotImplementedError('specification V and H not implemented') - def __call__(self, *, T=None, P=None, V=None, H=None, S=None, x=None, y=None, gas_reaction=None, liquid_reaction=None): """ @@ -401,19 +332,9 @@ def __call__(self, *, T=None, P=None, V=None, H=None, S=None, x=None, y=None, """ if gas_reaction or liquid_reaction: - reactants = self.imol[ - [i.reactant for i in (gas_reaction, liquid_reaction) - if i is not None] - ] - if reactants.any(): - if x is not None or y is not None: - raise ValueError( - "can not pass either 'x' or 'y' arguments with reactions preset" - ) - return self.reactive_vle( - liquid_reaction=liquid_reaction, - gas_reaction=gas_reaction, - T=T, P=P, V=V, H=H, S=S, + if x is not None or y is not None: + raise ValueError( + "can not pass either 'x' or 'y' arguments with reactions preset" ) ### Decide what kind of equilibrium to run ### @@ -432,21 +353,21 @@ def __call__(self, *, T=None, P=None, V=None, H=None, S=None, x=None, y=None, if T_spec: if P_spec: try: - self.set_thermal_condition(T, P) + self.set_thermal_condition(T, P, gas_reaction, liquid_reaction) except NoEquilibrium: thermal_condition = self._thermal_condition thermal_condition.T = T thermal_condition.P = P elif V_spec: try: - self.set_TV(T, V) + self.set_TV(T, V, gas_reaction, liquid_reaction) except NoEquilibrium: thermal_condition = self._thermal_condition thermal_condition.T = T elif H_spec: - self.set_TH(T, H) + self.set_TH(T, H, gas_reaction, liquid_reaction) elif S_spec: - self.set_TS(T, S) + self.set_TS(T, S, gas_reaction, liquid_reaction) elif x_spec: self.set_Tx(T, np.asarray(x)) else: # y_spec @@ -454,22 +375,22 @@ def __call__(self, *, T=None, P=None, V=None, H=None, S=None, x=None, y=None, elif P_spec: if V_spec: try: - self.set_PV(P, V) + self.set_PV(P, V, gas_reaction, liquid_reaction) except NoEquilibrium: thermal_condition = self._thermal_condition thermal_condition.P = P elif H_spec: try: - self.set_PH(P, H, stacklevel=1) + self.set_PH(P, H, gas_reaction, liquid_reaction) except NoEquilibrium: thermal_condition = self._thermal_condition thermal_condition.P = P elif S_spec: try: - self.set_PS(P, S, stacklevel=1) + self.set_PS(P, S, gas_reaction, liquid_reaction) except: try: - self.set_PS(P, S, stacklevel=1) + self.set_PS(P, S, gas_reaction, liquid_reaction) except NoEquilibrium: thermal_condition = self._thermal_condition thermal_condition.P = P @@ -500,68 +421,7 @@ def __call__(self, *, T=None, P=None, V=None, H=None, S=None, x=None, y=None, raise ValueError("specification V and x is invalid") else: # pragma: no cover raise ValueError("can only pass either 'x' or 'y' arguments, not both") - - def _setup(self): - # Get flow rates - imol = self._imol - self._phase_data = tuple(imol) - self._liquid_mol = liquid_mol = imol['l'] - self._vapor_mol = vapor_mol = imol['g'] - mol = liquid_mol + vapor_mol - nonzero = mol.nonzero_keys() - chemicals = self.chemicals - if self._nonzero == nonzero: - index = self._index - reset = False - else: - # Set up indices for both equilibrium and non-equilibrium species - index = chemicals.get_vle_indices(nonzero) - eq_chems = chemicals.tuple - eq_chems = [eq_chems[i] for i in index] - reset = True - self._nonzero = set(nonzero) - self._index = index - - # Get overall composition - if not mol.any(): raise NoEquilibrium('no chemicals to perform equilibrium') - self._F_mass = (chemicals.MW * mol).sum() - self._mol_vle = mol_vle = mol[index] - - # Set light and heavy keys - LNK_index = chemicals._light_indices - HNK_index = chemicals._heavy_indices - vapor_mol[HNK_index] = 0 - vapor_mol[LNK_index] = light_mol = mol[LNK_index] - liquid_mol[LNK_index] = 0 - liquid_mol[HNK_index] = heavy_mol = mol[HNK_index] - self._F_mol_light = F_mol_light = light_mol.sum() - self._F_mol_heavy = F_mol_heavy = (heavy_mol * chemicals._heavy_solutes).sum() - self._F_mol_vle = F_mol_vle = mol_vle.sum() - self._F_mol = F_mol = F_mol_vle + F_mol_light + F_mol_heavy - if F_mol == 0.: raise NoEquilibrium('no chemicals to perform equilibrium') - self._z = mol_vle / F_mol - self._z_light = z_light = F_mol_light / F_mol - self._z_heavy = z_heavy = F_mol_heavy / F_mol - self._z_norm = mol_vle / F_mol_vle - N = len(index) - if N: - N += z_light > 0. - N += z_heavy > 0. - self._N = N - if reset: - if N == 0: - self._phi = self._gamma = self._pcf = self._dew_point = self._bubble_point = None - elif N == 1: - self._chemical, = eq_chems - else: - # Set equilibrium objects - thermo = self._thermo - self._bubble_point = bp = self._bubble_point_cache(eq_chems, thermo) - self._dew_point = self._dew_point_cache(eq_chems, thermo) - self._pcf = bp.pcf - self._gamma = bp.gamma - self._phi = bp.phi - + @property def imol(self): return self._imol @@ -770,35 +630,51 @@ def set_Py(self, P, y): self._thermal_condition.T, x = self._dew_point.solve_Tx(y, P) self._lever_rule(x, y) - def set_thermal_condition(self, T, P): - self._setup() + def set_thermal_condition(self, T, P, gas_reaction=None, liquid_reaction=None): + self._setup(gas_reaction, liquid_reaction) thermal_condition = self._thermal_condition self._T = thermal_condition.T = T self._P = thermal_condition.P = P if self._N == 0: return if self._N == 1: return self._set_thermal_condition_chemical(T, P) # Check if there is equilibrium - P_dew, x_dew = self._dew_point.solve_Px(self._z, T) - if P <= P_dew and not self._F_mol_heavy: - self._vapor_mol[self._index] = self._mol_vle - self._liquid_mol[self._index] = 0 - return - P_bubble, y_bubble = self._bubble_point.solve_Py(self._z, T) - if P >= P_bubble and not self._F_mol_light: - self._vapor_mol[self._index] = 0 - self._liquid_mol[self._index] = self._mol_vle - return + if gas_reaction: + P_dew, dz_dew, x_dew, _ = self._dew_point.solve_Px(self._z, T, gas_reaction) + if P <= P_dew and not self._F_mol_heavy: + self._vapor_mol[self._index] = self._mol_vle + dz_dew * self._F_mol_vle + self._liquid_mol[self._index] = 0 + return + else: + P_dew, x_dew = self._dew_point.solve_Px(self._z, T) + if P <= P_dew and not self._F_mol_heavy: + self._vapor_mol[self._index] = self._mol_vle + self._liquid_mol[self._index] = 0 + return + dz_dew = None + if liquid_reaction: + P_bubble, dz_bubble, y_bubble, _ = self._bubble_point.solve_Py(self._z, T, liquid_reaction) + if P >= P_bubble and not self._F_mol_light: + self._vapor_mol[self._index] = 0 + self._liquid_mol[self._index] = self._mol_vle + dz_bubble * self._F_mol_vle + return + else: + P_bubble, y_bubble = self._bubble_point.solve_Py(self._z, T) + if P >= P_bubble and not self._F_mol_light: + self._vapor_mol[self._index] = 0 + self._liquid_mol[self._index] = self._mol_vle + return + dz_bubble = None # Guess composition in the vapor is a # weighted average of bubble/dew points dP = (P_bubble - P_dew) V = (P - P_dew) / dP if dP > 1. else 0.5 - self._refresh_v(V, y_bubble) - set_flows(self._vapor_mol, self._liquid_mol, self._index, self._solve_v(T, P), self._mol_vle) - + self._refresh_v(V, y_bubble, x_dew, dz_bubble, dz_dew) + v = self._solve_v(T, P, gas_reaction, liquid_reaction) + mol_vle = self._mol_vle + self._dmol_vle if (gas_reaction or liquid_reaction) else self._mol_vle + set_flows(self._vapor_mol, self._liquid_mol, self._index, v, mol_vle) - def set_TV(self, T, V): - self._setup() - mol = self._mol_vle + def set_TV(self, T, V, gas_reaction=None, liquid_reaction=None): + self._setup(gas_reaction, liquid_reaction) thermal_condition = self._thermal_condition thermal_condition.T = self._T = T if self._N == 0: raise RuntimeError('no chemicals present to perform VLE') @@ -806,33 +682,59 @@ def set_TV(self, T, V): if self._F_mol_heavy and V == 1.: V = 1. - 1e-3 if self._F_mol_light and V == 0.: V = 1e-3 if V == 1: - P_dew, x_dew = self._dew_point.solve_Px(self._z, T) - self._vapor_mol[self._index] = self._mol_vle + if gas_reaction: + P_dew, dz, y, x_dew = self._dew_point.solve_Px(self._z, T, gas_reaction) + self._vapor_mol[self._index] = self._mol_vle + dz * self._F_mol_vle + else: + P_dew, x_dew = self._dew_point.solve_Px(self._z, T) + self._vapor_mol[self._index] = self._mol_vle self._liquid_mol[self._index] = 0 thermal_condition.P = P_dew elif V == 0: - P_bubble, y_bubble = self._bubble_point.solve_Py(self._z, T) + if liquid_reaction: + P_bubble, dz, y_bubble, x = self._bubble_point.solve_Py(self._z, T) + self._liquid_mol[self._index] = self._mol_vle + dz * self._F_mol_vle + else: + P_bubble, y_bubble = self._bubble_point.solve_Py(self._z, T) + self._liquid_mol[self._index] = self._mol_vle self._vapor_mol[self._index] = 0 - self._liquid_mol[self._index] = self._mol_vle thermal_condition.P = P_bubble else: - P_bubble, y_bubble = self._bubble_point.solve_Py(self._z, T) - P_dew, x_dew = self._dew_point.solve_Px(self._z, T) - self._refresh_v(V, y_bubble) - if self._F_mol_light: P_bubble = self._bubble_point.Pmax - if self._F_mol_heavy: P_dew = self._bubble_point.Pmin + if liquid_reaction: + P_bubble, dz_bubble, y_bubble, _ = self._bubble_point.solve_Py(self._z, T, liquid_reaction) + else: + P_bubble, y_bubble = self._bubble_point.solve_Py(self._z, T) + dz_bubble = None + if gas_reaction: + P_dew, dz_dew, y, x_dew = self._dew_point.solve_Px(self._z, T, gas_reaction) + else: + P_dew, x_dew = self._dew_point.solve_Px(self._z, T, gas_reaction) + dz_dew = None + self._refresh_v(V, y_bubble, x_dew, dz_bubble, dz_dew) + if self._F_mol_light: P_bubble = 0.1 * self._bubble_point.Pmax + 0.9 * P_bubble + if self._F_mol_heavy: P_dew = 0.1 * self._bubble_point.Pmin + 0.9 * P_dew - V_bubble = self._V_err_at_P(P_bubble, 0.) + V_bubble = self._V_err_at_P(P_bubble, 0., gas_reaction, liquid_reaction) if V_bubble > V: - F_mol_vapor = self._F_mol * V + F_mol = self._F_mol + mol = self._mol_vle + if liquid_reaction: + mol = mol + self._dmol_vle + F_mol = F_mol + self._dF_mol + F_mol_vapor = F_mol * V v = y_bubble * F_mol_vapor mask = v > mol v[mask] = mol[mask] P = P_bubble else: - V_dew = self._V_err_at_P(P_dew, 0.) + V_dew = self._V_err_at_P(P_dew, 0., gas_reaction, liquid_reaction) if V_dew < V: - l = x_dew * self._F_mol * (1. - V) + F_mol = self._F_mol + mol = self._mol_vle + if gas_reaction: + mol = mol + self._dmol_vle + F_mol = F_mol + self._dF_mol + l = x_dew * F_mol * (1. - V) mask = l > mol l[mask] = mol[mask] v = mol - l @@ -842,18 +744,28 @@ def set_TV(self, T, V): self._V_err_at_P, P_bubble, P_dew, V_bubble - V, V_dew - V, self._P, self.P_tol, self.V_tol, - (V,), checkiter=False, checkbounds=False, + (V, gas_reaction, liquid_reaction), + checkiter=False, checkbounds=False, maxiter=self.maxiter, ) v = self._v + mol = self._mol_vle + if gas_reaction or liquid_reaction: + mol = mol + self._dmol_vle self._P = thermal_condition.P = P set_flows(self._vapor_mol, self._liquid_mol, self._index, v, mol) - try: self._H_hat = self.mixture.xH(self._phase_data, T, P) / self._F_mass - except: pass + if liquid_reaction or gas_reaction: + try: self._H_hat = ( + self.mixture.xH(self._phase_data, T, P) + (self.chemicals.Hf * mol).sum() + ) / self._F_mass + except: pass + else: + try: self._H_hat = self.mixture.xH(self._phase_data, T, P) / self._F_mass + except: pass - def set_TH(self, T, H): - self._setup() + def set_TH(self, T, H, gas_reaction=None, liquid_reaction=None): + self._setup(gas_reaction, liquid_reaction) if self._N == 0: raise RuntimeError('no chemicals present to perform VLE') if self._N == 1: return self._set_TH_chemical(T, H) self._T = T @@ -887,7 +799,7 @@ def set_TH(self, T, H): V = dH_bubble/(H_dew - H_bubble) # Guess composition in the vapor is a weighted average of boiling points - self._refresh_v(V, y_bubble) + self._refresh_v(V, y_bubble, x_dew) F_mass = self._F_mass H_hat = H/F_mass P = flx.IQ_interpolation( @@ -901,8 +813,8 @@ def set_TH(self, T, H): self._P = self._thermal_condition.P = P self._thermal_condition.T = T - def set_TS(self, T, S): - self._setup() + def set_TS(self, T, S, gas_reaction=None, liquid_reaction=None): + self._setup(gas_reaction, liquid_reaction) if self._N == 0: raise RuntimeError('no chemicals present to perform VLE') if self._N == 1: return self._set_TS_chemical(T, S) self._T = T @@ -936,7 +848,7 @@ def set_TS(self, T, S): V = dS_bubble/(S_dew - S_bubble) # Guess composition in the vapor is a weighted average of boiling points - self._refresh_v(V, y_bubble) + self._refresh_v(V, y_bubble, x_dew) F_mass = self._F_mass S_hat = S/F_mass P = flx.IQ_interpolation( @@ -950,8 +862,8 @@ def set_TS(self, T, S): self._P = self._thermal_condition.P = P self._thermal_condition.T = T - def set_PV(self, P, V): - self._setup() + def set_PV(self, P, V, gas_reaction=None, liquid_reaction=None): + self._setup(gas_reaction, liquid_reaction) self._thermal_condition.P = self._P = P if self._N == 0: raise RuntimeError('no chemicals present to perform VLE') if self._N == 1: return self._set_PV_chemical(P, V) @@ -959,38 +871,63 @@ def set_PV(self, P, V): # Setup bounderies thermal_condition = self._thermal_condition index = self._index - mol = self._mol_vle vapor_mol = self._vapor_mol liquid_mol = self._liquid_mol if self._F_mol_heavy and V == 1.: V = 1. - 1e-3 if self._F_mol_light and V == 0.: V = 1e-3 - if V == 1: - T_dew, x_dew = self._dew_point.solve_Tx(self._z, P) - vapor_mol[index] = mol - liquid_mol[index] = 0 + if V == 1 and not self._F_mol_heavy: + if gas_reaction: + T_dew, dz, y, x_dew = self._dew_point.solve_Tx(self._z, P, gas_reaction) + self._vapor_mol[self._index] = self._mol_vle + dz * self._F_mol_vle + else: + T_dew, x_dew = self._dew_point.solve_Tx(self._z, P) + self._vapor_mol[self._index] = self._mol_vle + self._liquid_mol[self._index] = 0 thermal_condition.T = T_dew elif V == 0 and not self._F_mol_light: - T_bubble, y_bubble = self._bubble_point.solve_Ty(self._z, P) - vapor_mol[index] = 0 - liquid_mol[index] = mol + if liquid_reaction: + T_bubble, dz, y_bubble, x = self._bubble_point.solve_Ty(self._z, P, liquid_reaction) + self._liquid_mol[self._index] = self._mol_vle + dz * self._F_mol_vle + else: + T_bubble, y_bubble = self._bubble_point.solve_Ty(self._z, P) + self._liquid_mol[self._index] = self._mol_vle + self._vapor_mol[self._index] = 0 thermal_condition.T = T_bubble else: - T_dew, x_dew = self._dew_point.solve_Tx(self._z, P) - T_bubble, y_bubble = self._bubble_point.solve_Ty(self._z, P) - self._refresh_v(V, y_bubble) - if self._F_mol_heavy: T_dew = 0.9 * T_dew + 0.1 * self._dew_point.Tmax - if self._F_mol_light: T_bubble = 0.9 * T_bubble + 0.1 * self._bubble_point.Tmin - V_bubble = self._V_err_at_T(T_bubble, 0.) + if liquid_reaction: + T_bubble, dz_bubble, y_bubble, _ = self._bubble_point.solve_Ty(self._z, P, liquid_reaction) + else: + T_bubble, y_bubble = self._bubble_point.solve_Ty(self._z, P) + dz_bubble = None + if gas_reaction: + T_dew, dz_dew, y, x_dew = self._dew_point.solve_Tx(self._z, P, gas_reaction) + else: + T_dew, x_dew = self._dew_point.solve_Tx(self._z, P) + dz_dew = None + self._refresh_v(V, y_bubble, x_dew, dz_bubble, dz_dew) + if self._F_mol_light: T_bubble = 0.1 * self._bubble_point.Tmin + 0.9 * T_bubble + if self._F_mol_heavy: T_dew = 0.1 * self._bubble_point.Tmax + 0.9 * T_dew + V_bubble = self._V_err_at_T(T_bubble, 0., gas_reaction, liquid_reaction) if V_bubble > V: - F_mol_vapor = self._F_mol * V + F_mol = self._F_mol + mol = self._mol_vle + if liquid_reaction: + mol = mol + self._dmol_vle + F_mol = F_mol + self._dF_mol + F_mol_vapor = F_mol * V v = y_bubble * F_mol_vapor mask = v > mol v[mask] = mol[mask] T = T_bubble else: - V_dew = self._V_err_at_T(T_dew, 0.) + V_dew = self._V_err_at_T(T_dew, 0., gas_reaction, liquid_reaction) if V_dew < V: - l = x_dew * self._F_mol * (1. - V) + F_mol = self._F_mol + mol = self._mol_vle + if gas_reaction: + mol = mol + self._dmol_vle + F_mol = F_mol + self._dF_mol + l = x_dew * F_mol * (1. - V) mask = l > mol l[mask] = mol[mask] v = mol - l @@ -1000,18 +937,27 @@ def set_PV(self, P, V): self._V_err_at_T, T_bubble, T_dew, V_bubble - V, V_dew - V, self._T, self.T_tol, self.V_tol, - (V,), checkiter=False, checkbounds=False, + (V, gas_reaction, liquid_reaction), + checkiter=False, checkbounds=False, maxiter=self.maxiter, ) - v = self._v + mol = self._mol_vle + if gas_reaction or liquid_reaction: + mol = mol + self._dmol_vle self._T = thermal_condition.T = T set_flows(vapor_mol, liquid_mol, index, v, mol) - try: self._H_hat = self.mixture.xH(self._phase_data, T, P)/self._F_mass - except: pass + if liquid_reaction or gas_reaction: + try: self._H_hat = ( + self.mixture.xH(self._phase_data, T, P) + (self.chemicals.Hf * mol).sum(0) + ) / self._F_mass + except: pass + else: + try: self._H_hat = self.mixture.xH(self._phase_data, T, P) / self._F_mass + except: pass - def set_PS(self, P, S, stacklevel=0): - self._setup() + def set_PS(self, P, S, gas_reaction=None, liquid_reaction=None): + self._setup(gas_reaction, liquid_reaction) thermal_condition = self._thermal_condition thermal_condition.P = self._P = P if self._N == 0: @@ -1055,7 +1001,7 @@ def set_PS(self, P, S, stacklevel=0): # Guess T, overall vapor fraction, and vapor flow rates V = dS_bubble/(S_dew - S_bubble) - self._refresh_v(V, y_bubble) + self._refresh_v(V, y_bubble, x_dew) F_mass = self._F_mass S_hat = S/F_mass @@ -1124,8 +1070,8 @@ def set_PS(self, P, S, stacklevel=0): ) self._S_hat = S_hat - def set_PH(self, P, H, stacklevel=0): - self._setup() + def set_PH(self, P, H, gas_reaction=None, liquid_reaction=None): + self._setup(gas_reaction, liquid_reaction) thermal_condition = self._thermal_condition thermal_condition.P = self._P = P if self._N == 0: @@ -1137,23 +1083,38 @@ def set_PH(self, P, H, stacklevel=0): # Setup bounderies index = self._index - mol = self._mol_vle vapor_mol = self._vapor_mol liquid_mol = self._liquid_mol # Check if subcooled liquid - T_bubble, y_bubble = self._bubble_point.solve_Ty(self._z, P) + if liquid_reaction: + T_bubble, dz_bubble, y_bubble, x = self._bubble_point.solve_Ty(self._z, P, liquid_reaction) + mol = self._mol_vle + dz_bubble * self._F_mol_vle + else: + T_bubble, y_bubble = self._bubble_point.solve_Ty(self._z, P) + mol = self._mol_vle + dz_bubble = None if self._F_mol_light: T_bubble = 0.9 * T_bubble + 0.1 * self._bubble_point.Tmin vapor_mol[index] = 0 liquid_mol[index] = mol H_bubble = self.mixture.xH(self._phase_data, T_bubble, P) + if liquid_reaction: + Hf = (self.chemicals.Hf * mol).sum() + H_bubble += Hf dH_bubble = H - H_bubble if dH_bubble <= 0: + if liquid_reaction: H -= Hf thermal_condition.T = self.mixture.xsolve_T_at_HP(self._phase_data, H, T_bubble, P) return # Check if super heated vapor - T_dew, x_dew = self._dew_point.solve_Tx(self._z, P) + if gas_reaction: + T_dew, dz_dew, y, x_dew = self.dew_point.solve_Tx(self._z, P, gas_reaction) + mol = self._mol_vle + dz_dew + else: + T_dew, x_dew = self._dew_point.solve_Tx(self._z, P) + mol = self._mol_vle + dz_dew = None if T_dew <= T_bubble: T_dew, T_bubble = T_bubble, T_dew T_dew += 0.5 @@ -1162,31 +1123,39 @@ def set_PH(self, P, H, stacklevel=0): vapor_mol[index] = mol liquid_mol[index] = 0 H_dew = self.mixture.xH(self._phase_data, T_dew, P) + if gas_reaction: + Hf = (self.chemicals.Hf * mol).sum() + H_dew += Hf dH_dew = H - H_dew if dH_dew >= 0: + if gas_reaction: H -= Hf thermal_condition.T = self.mixture.xsolve_T_at_HP(self._phase_data, H, T_dew, P) return # Guess T, overall vapor fraction, and vapor flow rates - V = dH_bubble/(H_dew - H_bubble) - self._refresh_v(V, y_bubble) + V = abs(dH_bubble / (H_dew - H_bubble)) + self._refresh_v(V, y_bubble, x_dew, dz_bubble, dz_dew) F_mass = self._F_mass H_hat = H/F_mass - if (H_hat_bubble:=self._H_hat_err_at_T(T_bubble, 0.)) > H_hat: + if (H_hat_bubble:=self._H_hat_err_at_T(T_bubble, 0., gas_reaction, liquid_reaction)) > H_hat: T = T_bubble - elif (H_hat_dew:=self._H_hat_err_at_T(T_dew, 0.)) < H_hat: + elif (H_hat_dew:=self._H_hat_err_at_T(T_dew, 0., gas_reaction, liquid_reaction)) < H_hat: T = T_dew else: T = flx.IQ_interpolation( self._H_hat_err_at_T, T_bubble, T_dew, H_hat_bubble - H_hat, H_hat_dew - H_hat, self._T, self.T_tol, self.H_hat_tol, - (H_hat,), checkiter=False, checkbounds=False, + (H_hat, gas_reaction, liquid_reaction), + checkiter=False, checkbounds=False, maxiter=self.maxiter ) - + + if gas_reaction or liquid_reaction: + H -= (self.chemicals.Hf * (self._mol_vle + self._dmol_vle)).sum() + # Make sure energy balance is correct by vaporizing a fraction # of the liquid or condensing a fraction of the vapor self._T = thermal_condition.T = T @@ -1239,16 +1208,25 @@ def set_PH(self, P, H, stacklevel=0): ) self._H_hat = H_hat - def _estimate_v(self, V, y_bubble): - return (V*self._z_norm + (1-V)*y_bubble) * V * self._F_mol_vle + def _estimate_v(self, V, y_bubble, x_dew, dz_bubble, dz_dew): + F_mol_vle = self._F_mol_vle + z = self._z_norm + L = (1 - V) + if dz_bubble is not None: + F_mol_vle = F_mol_vle - dz_bubble * L * F_mol_vle + if dz_dew is not None: + F_mol_vle = F_mol_vle - dz_dew * V * F_mol_vle + v_bubble = (V * z + L * y_bubble) * V * F_mol_vle + v_dew = (z - L * (L * z + V * x_dew)) * F_mol_vle + v = L * v_bubble + V * v_dew + return v - def _refresh_v(self, V, y_bubble): - self._v = v = self._estimate_v(V, y_bubble) + def _refresh_v(self, V, y_bubble, x_dew, dz_bubble=None, dz_dew=None): # TODO: use reaction data here for better estimate + self._v = v = self._estimate_v(V, y_bubble, x_dew, dz_bubble, dz_dew) self._V = V self._y = fn.normalize(v, v.sum() + self._F_mol_light) - z_last = self._z_last try: - reload_cache = self._x is None or np.abs(z_last - self._z).sum() > 0.001 + reload_cache = self._x is None or np.abs(self._z_last - self._z).sum() > 0.001 except: reload_cache = True if reload_cache: @@ -1256,16 +1234,30 @@ def _refresh_v(self, V, y_bubble): l[l < 0] = 1e-12 self._x = fn.normalize(l, l.sum() + self._F_mol_heavy) - def _H_hat_err_at_T(self, T, H_hat): - set_flows(self._vapor_mol, self._liquid_mol, self._index, - self._solve_v(T, self._P), self._mol_vle) - self._H_hat = self.mixture.xH(self._phase_data, T, self._P)/self._F_mass + def _H_hat_err_at_T(self, T, H_hat, gas_reaction, liquid_reaction): + v = self._solve_v(T, self._P, gas_reaction, liquid_reaction) + if gas_reaction or liquid_reaction: + mol_vle = self._mol_vle + self._dmol_vle + set_flows(self._vapor_mol, self._liquid_mol, self._index, v, mol_vle) + H = self.mixture.xH(self._phase_data, T, self._P) + (self.chemicals.Hf * mol_vle).sum() + else: + mol_vle = self._mol_vle + set_flows(self._vapor_mol, self._liquid_mol, self._index, v, mol_vle) + H = self.mixture.xH(self._phase_data, T, self._P) + self._H_hat = H / self._F_mass return self._H_hat - H_hat - def _H_hat_err_at_P(self, P, H_hat): - set_flows(self._vapor_mol, self._liquid_mol, self._index, - self._solve_v(self._T, P), self._mol_vle) - self._H_hat = self.mixture.xH(self._phase_data, self._T, P)/self._F_mass + def _H_hat_err_at_P(self, P, H_hat, gas_reaction=None, liquid_reaction=None): + v = self._solve_v(self._T, P, gas_reaction, liquid_reaction) + if gas_reaction or liquid_reaction: + mol_vle = self._mol_vle + self._dmol_vle + set_flows(self._vapor_mol, self._liquid_mol, self._index, v, mol_vle) + H = self.mixture.xH(self._phase_data, self._T, P) + (self.chemicals.Hf * mol_vle).sum() + else: + mol_vle = self._mol_vle + set_flows(self._vapor_mol, self._liquid_mol, self._index, v, mol_vle) + H = self.mixture.xH(self._phase_data, self._T, P) + self._H_hat = H / self._F_mass return self._H_hat - H_hat def _S_hat_err_at_T(self, T, S_hat): @@ -1280,50 +1272,28 @@ def _S_hat_err_at_P(self, P, S_hat): self._S_hat = self.mixture.xS(self._phase_data, self._T, P)/self._F_mass return self._S_hat - S_hat - def _V_err_at_P(self, P, V): - return self._solve_v(self._T , P).sum()/self._F_mol_vle - V - - def _V_err_at_T(self, T, V): - return self._solve_v(T, self._P).sum()/self._F_mol_vle - V - - def _solve_y(self, y, pcf_Psat_over_P, T, P): - gamma = self._gamma - x = self._x - pcf_Psat_over_P = pcf_Psat_over_P - N = self._N - z = self._z - Ks = pcf_Psat_over_P.copy() - n = z.size - if N > 2 or self._z_light or self._z_heavy: - f = xV_iter - args = (pcf_Psat_over_P, T, P, z, self._z_light, - self._z_heavy, gamma.f, gamma.args, self._phi, n) - elif N == 2: - f = xV_iter_2n - args = (pcf_Psat_over_P, T, P, z, gamma.f, gamma.args, self._phi, n) - xVlogK = np.zeros(2*x.size + 1) - xVlogK[:n] = x - xVlogK[n] = self._V - xVlogK[-n:] = np.log(Ks) - xVlogK = flx.aitken(f, xVlogK, self.x_tol, args, checkiter=False, - checkconvergence=False, convergenceiter=5, - maxiter=self.maxiter, subset=n) - x = xVlogK[:n] - self._V = V = xVlogK[n] - x[x < 1e-32] = 1e-32 - self._x = x = fn.normalize(x) - if V == 0: - Ks = 0 + def _V_err_at_P(self, P, V, gas_reaction, liquid_reaction): + v = self._solve_v(self._T , P, gas_reaction, liquid_reaction).sum() + if gas_reaction or liquid_reaction: + F_mol_vle = self._F_mol_vle + self._dF_mol else: - Ks = (z / x - 1) / V + 1. - self._z_last = z - v = self._F_mol * V * x * Ks - return fn.normalize(v, v.sum() + self._F_mol_light) + F_mol_vle = self._F_mol_vle + return v / F_mol_vle - V + + def _V_err_at_T(self, T, V, gas_reaction, liquid_reaction): + v = self._solve_v(T, self._P, gas_reaction, liquid_reaction).sum() + if gas_reaction or liquid_reaction: + F_mol_vle = self._F_mol_vle + self._dF_mol + else: + F_mol_vle = self._F_mol_vle + return v / F_mol_vle - V - def _solve_v(self, T, P): + def _solve_v(self, T, P, gas_reaction=None, liquid_reaction=None): """Solve for vapor mol""" method = self.method if method == 'shgo': + if gas_reaction or liquid_reaction: + raise RuntimeError('shgo is not a valid method (yet) when reactions are present') gamma = self._gamma phi = self._phi Psats = np.array([i(T) for i in self._bubble_point.Psats]) @@ -1341,28 +1311,13 @@ def _solve_v(self, T, P): self._bubble_point.Psats]) pcf_Psats_over_P = self._pcf(T, P, Psats) * Psats / P self._T = T - y = self._solve_y(self._y, pcf_Psats_over_P, T, P) - self._v = v = self._F_mol * self._V * y - mask = v > self._mol_vle - v[mask] = self._mol_vle[mask] - v[v < 0.] = 0. - else: - raise RuntimeError(f"invalid method '{method}'") - return v - - def _solve_v_reactive(self, T, P, gas_reaction, liquid_reaction): - """Solve for vapor mol""" - method = self.method - if method == 'shgo': - raise NotImplementedError("'shgo' method not implemented yet") - elif method == 'fixed-point': - Psats = np.array([i(T) for i in - self._bubble_point.Psats]) - pcf_Psats_over_P = self._pcf(T, P, Psats) * Psats / P - self._T = T - y = self._solve_y_reactive(self._y, pcf_Psats_over_P, T, P, gas_reaction, liquid_reaction) - self._v = v = (self._F_mol + self._dF_mol) * self._V * y - mol_vle = self._mol_vle + self._dmol_vle + y = self._solve_y(self._y, pcf_Psats_over_P, T, P, gas_reaction, liquid_reaction) + if gas_reaction or liquid_reaction: + self._v = v = (self._F_mol + self._dF_mol) * self._V * y + mol_vle = self._mol_vle + self._dmol_vle + else: + self._v = v = self._F_mol * self._V * y + mol_vle = self._mol_vle mask = v > mol_vle v[mask] = mol_vle[mask] v[v < 0.] = 0. @@ -1370,18 +1325,25 @@ def _solve_v_reactive(self, T, P, gas_reaction, liquid_reaction): raise RuntimeError(f"invalid method '{method}'") return v - def _solve_y_reactive(self, y, pcf_Psat_over_P, T, P, gas_reaction, liquid_reaction): + def _solve_y(self, y, pcf_Psat_over_P, T, P, gas_reaction, liquid_reaction): gamma = self._gamma x = self._x pcf_Psat_over_P = pcf_Psat_over_P z = self._z Ks = pcf_Psat_over_P.copy() - f = xV_reactive_iter index = self._index n = z.size - args = (pcf_Psat_over_P, T, P, z, self._z_light, self._z_heavy, - gamma.f, gamma.args, self._phi, n, - gas_reaction, liquid_reaction, index) + if n > 2 or self._z_light or self._z_heavy: + f = xV_iter + args = (pcf_Psat_over_P, T, P, z, + self._z_light, self._z_heavy, + gamma.f, gamma.args, self._phi, n, + gas_reaction, liquid_reaction, index) + else: + f = xV_iter_2n + args = (pcf_Psat_over_P, T, P, z, + gamma.f, gamma.args, self._phi, n, + gas_reaction, liquid_reaction, index) xVlogK = np.zeros(2*x.size + 1) xVlogK[:n] = x xVlogK[n] = self._V @@ -1400,36 +1362,32 @@ def _solve_y_reactive(self, y, pcf_Psat_over_P, T, P, gas_reaction, liquid_react else: Ks = (z / x - 1) / V + 1. self._z_last = z - dmol_vle = 0 - if gas_reaction: - material = SparseVector.from_size(gas_reaction.chemicals.size) - material[index] = y * V - dmol_vle += gas_reaction.conversion(material=material, T=T, P=P, phase='g')[index] - if liquid_reaction: - material = SparseVector.from_size(liquid_reaction.chemicals.size) - material[index] = x * (1 - V) - dmol_vle += liquid_reaction.conversion(material=material, T=T, P=P, phase='l')[index] - dmol_vle *= self._F_mol - self._dmol_vle = dmol_vle - self._dF_mol = dmol_vle.sum() - v = (self._F_mol - self._dF_mol) * V * x * Ks + if gas_reaction or liquid_reaction: + dmol_vle = 0 + if gas_reaction: + dmol_vle += gas_reaction.conversion(material=y * V, T=T, P=P, phase='g') + if liquid_reaction: + dmol_vle += liquid_reaction.conversion(material=x * (1 - V), T=T, P=P, phase='l') + dmol_vle *= self._F_mol + self._dmol_vle = dmol_vle + self._dF_mol = dmol_vle.sum() + v = (self._F_mol - self._dF_mol) * V * x * Ks + else: + v = self._F_mol * V * x * Ks return fn.normalize(v, v.sum() + self._F_mol_light) - - def _setup_reactive(self, gas_reaction, liquid_reaction): + def _setup(self, gas_reaction=None, liquid_reaction=None): imol = self._imol self._phase_data = tuple(imol) self._liquid_mol = liquid_mol = imol['l'] self._vapor_mol = vapor_mol = imol['g'] mol = liquid_mol + vapor_mol nonzero = set(mol.nonzero_keys()) - has_gas_reaction = gas_reaction and imol['g', gas_reaction.reactant] - has_liquid_reaction = liquid_reaction and imol['l', liquid_reaction.reactant] - if has_gas_reaction: + if gas_reaction: nonzero.update( gas_reaction.stoichiometry.nonzero_keys() ) - if has_liquid_reaction: + if liquid_reaction: nonzero.update( liquid_reaction.stoichiometry.nonzero_keys() ) @@ -1445,8 +1403,13 @@ def _setup_reactive(self, gas_reaction, liquid_reaction): reset = True self._nonzero = set(nonzero) self._index = index + if gas_reaction and [i.ID for i in eq_chems] != [i.ID for i in gas_reaction.chemicals]: + gas_reaction.reset_thermo(gas_reaction.thermo.subset(eq_chems)) + if liquid_reaction and [i.ID for i in eq_chems] != [i.ID for i in liquid_reaction.chemicals]: + liquid_reaction.reset_thermo(liquid_reaction.thermo.subset(eq_chems)) # Get overall composition + if not mol.any(): raise NoEquilibrium('no chemicals to perform equilibrium') self._F_mass = (chemicals.MW * mol).sum() self._mol_vle = mol_vle = mol[index] @@ -1462,7 +1425,9 @@ def _setup_reactive(self, gas_reaction, liquid_reaction): self._F_mol_vle = F_mol_vle = mol_vle.sum() self._F_mol = F_mol = F_mol_vle + F_mol_light + F_mol_heavy self._z = mol_vle / F_mol - self._z_light = z_light = F_mol_light / F_mol + try: self._z_light = z_light = F_mol_light / F_mol + except: + breakpoint() self._z_heavy = z_heavy = F_mol_heavy / F_mol self._z_norm = mol_vle / F_mol_vle N = len(index) @@ -1478,28 +1443,12 @@ def _setup_reactive(self, gas_reaction, liquid_reaction): else: # Set equilibrium objects thermo = self._thermo - self._bubble_point = bp = self._bubble_point_cache(eq_chems, thermo) - self._dew_point = self._dew_point_cache(eq_chems, thermo) + self._bubble_point = bp = BubblePoint(eq_chems, thermo) + self._dew_point = DewPoint(eq_chems, thermo) self._pcf = bp.pcf self._gamma = bp.gamma self._phi = bp.phi - - - def set_thermal_condition_reactive(self, T, P, gas_reaction, liquid_reaction): - self._setup_reactive(gas_reaction, liquid_reaction) - thermal_condition = self._thermal_condition - self._T = thermal_condition.T = T - self._P = thermal_condition.P = P - - # Guess composition in the vapor is a - # weighted average of bubble/dew points - v = self._vapor_mol[self._index] - self._refresh_v(self._V, v / v.sum()) - set_flows( - self._vapor_mol, self._liquid_mol, - self._index, self._solve_v_reactive(T, P, gas_reaction, liquid_reaction), - self._mol_vle + self._dmol_vle - ) + class VLECache(Cache): load = VLE del Cache, Equilibrium diff --git a/thermosteam/reaction/_reaction.py b/thermosteam/reaction/_reaction.py index e8d45257..a567583b 100644 --- a/thermosteam/reaction/_reaction.py +++ b/thermosteam/reaction/_reaction.py @@ -17,7 +17,7 @@ _parse as prs, _xparse as xprs, ) -from ..utils import chemicals_user, AbstractMethod +from ..utils import thermo_user, AbstractMethod from .._phase import NoPhase, phase_tuple from ..indexer import ChemicalIndexer, MaterialIndexer from ..exceptions import InfeasibleRegion @@ -99,7 +99,7 @@ def as_material_array(material, basis, phases, chemicals): # %% Stoichiometric reactions -@chemicals_user +@thermo_user class Reaction: """ Create a Reaction object which defines a stoichiometric reaction and @@ -247,19 +247,20 @@ class Reaction: Glucose + O2 -> Ethanol + CO2 Glucose 35.00 """ - _kinetics = AbstractMethod + _rate = AbstractMethod phases = MaterialIndexer.phases __slots__ = ( + '_stream', '_basis', '_phases', - '_chemicals', + '_thermo', '_X_index', '_stoichiometry', '_X' ) def __init__(self, reaction, reactant=None, X=1., - chemicals=None, basis='mol', *, + chemicals=None, basis='mol', thermo=None, *, phases=None, check_mass_balance=False, check_atomic_balance=False, @@ -270,7 +271,8 @@ def __init__(self, reaction, reactant=None, X=1., else: raise ValueError("basis must be either by 'wt' or by 'mol'") self.X = X - chemicals = self._load_chemicals(chemicals) + thermo = self._load_thermo(thermo or chemicals) + chemicals = self.chemicals if reaction: self._phases = phases = phase_tuple(phases) if phases else xprs.get_phases(reaction) if phases: @@ -283,7 +285,7 @@ def __init__(self, reaction, reactant=None, X=1., else: raise ValueError('must pass reactant when multiple reactants are involved') else: - reactant_index = self._chemicals.index(reactant) + reactant_index = self.chemicals.index(reactant) for phase_index, x in enumerate(stoichiometry[:, reactant_index]): if x: break self._X_index = (phase_index, reactant_index) @@ -297,7 +299,7 @@ def __init__(self, reaction, reactant=None, X=1., else: raise ValueError('must pass reactant when multiple reactants are involved') else: - self._X_index = self._chemicals.index(reactant) + self._X_index = self.chemicals.index(reactant) self._rescale() if correct_atomic_balance: self.correct_atomic_balance() @@ -311,7 +313,7 @@ def __init__(self, reaction, reactant=None, X=1., else: self._phases = () self._stoichiometry = SparseVector.from_size(chemicals.size) - self._X_index = self._chemicals.index(reactant) + self._X_index = self.chemicals.index(reactant) def backwards(self, reactant=None, X=None): new = self.copy() @@ -331,7 +333,7 @@ def backwards(self, reactant=None, X=None): else: raise ValueError('must pass reactant when multiple reactants are involved') else: - new._X_index = new._chemicals.index(reactant) + new._X_index = new.chemicals.index(reactant) if X is not None: new.X = X new._rescale() return new @@ -340,31 +342,36 @@ def backwards(self, reactant=None, X=None): def reaction_chemicals(self): """Return all chemicals involved in the reaction.""" keys = sorted(self._stoichiometry.nonzero_keys()) - chemicals = self._chemicals + chemicals = self.chemicals return [chemicals[i] for i in keys] - def reset_chemicals(self, chemicals): - if self._chemicals is chemicals: return + def reset_thermo(self, thermo): + if self._thermo is thermo: return + chemicals = thermo.chemicals phases = self.phases stoichiometry = self._stoichiometry reactant = self.reactant if phases: M, N = stoichiometry.shape new_stoichiometry = SparseArray.from_shape([M, chemicals.size]) - IDs = self._chemicals.IDs + IDs = self.chemicals.IDs for (i, j), k in stoichiometry.nonzero_items(): new_stoichiometry[i, chemicals.index(IDs[j])] = k phase, reactant = reactant X_index = phases.index(phase), chemicals.index(reactant) else: new_stoichiometry = SparseVector.from_size(chemicals.size) - IDs = self._chemicals.IDs + IDs = self.chemicals.IDs for i, j in stoichiometry.nonzero_items(): new_stoichiometry[chemicals.index(IDs[i])] = j X_index = chemicals.index(reactant) - self._chemicals = chemicals + self._thermo = thermo self._stoichiometry = new_stoichiometry self._X_index = X_index + if hasattr(self, '_stream'): del self._stream + + def reset_chemicals(self, chemicals): + raise Exception('`reset_chemicals` is deprecated; use `reset_thermo` instead') def copy(self, basis=None): """Return copy of Reaction object.""" @@ -373,7 +380,7 @@ def copy(self, basis=None): copy._phases = self._phases copy._stoichiometry = self._stoichiometry.copy() copy._X_index = self._X_index - copy._chemicals = self._chemicals + copy._thermo = self._thermo copy._X = self._X if basis: set_reaction_basis(copy, basis) return copy @@ -385,7 +392,7 @@ def has_reaction(self): def _math_compatible_reaction(self, rxn, copy=True): basis = self.basis if copy or basis != rxn._basis: rxn = rxn.copy(basis) - if self._chemicals is not rxn._chemicals: + if self.chemicals is not rxn.chemicals: raise ValueError('chemicals must be the same to add/subtract reactions') if self._phases != rxn._phases: raise ValueError('phases must be the same to add/subtract reactions') @@ -451,13 +458,15 @@ def __isub__(self, rxn): self.X = self.X - rxn.X return self - def kinetics(self, material, time=None, T=None, P=None, phase=None): + def _get_stream(self, material, T=None, P=None, phase=None): if isinstance(material, tmo.Stream): try: stream = self._stream except: - if self._chemicals is material.chemicals: + if self._thermo is material._thermo: self._stream = stream = material.copy() + elif self.chemicals is material.chemicals: + self._stream = stream = tmo.Stream(thermo=self._thermo) else: raise ValueError('incompatible chemicals for kinetics') stream.copy_like(material) @@ -468,31 +477,59 @@ def kinetics(self, material, time=None, T=None, P=None, phase=None): try: stream = self._stream except: - thermo = tmo.settings.get_thermo() - if self._chemicals is thermo.chemicals: - self._stream = stream = tmo.Stream() - else: - raise ValueError('incompatible chemicals for kinetics') + self._stream = stream = tmo.Stream(thermo=self._thermo) stream.imol.data[:] = material stream.T = T stream.P = P stream.phase = phase - self._kinetics(stream, time) + return stream - def conversion(self, material, time=None, T=None, P=None, phase=None): + def rate(self, material, T=None, P=None, phase=None): + return self._rate(self._get_stream(material, T, P, phase)) + + def conversion_bounds(self, material): + material, *_ = as_material_array(material, self._basis, self._phases, self.chemicals) + X = -material / (material[self._X_index] * self._stoichiometry) + stoichiometry = self.stoichiometry + return X[stoichiometry > 0].max(), X[stoichiometry < 0].min() + + def _equilibrium_objective(self, X, data): + stream = self._stream + self.X = X + self._reaction(stream.mol) + rate = self._rate(stream) + stream.set_data(data) + return rate + + def equilibrium(self, material, T, P, phase): + stream = self._get_stream(material, T, P, phase) + initial_condition = stream.get_data() + f = self._equilibrium_objective + args = (initial_condition,) + self.X = flx.IQ_interpolation(f, *self.conversion_bounds(stream), args=args) + + def conversion(self, material, T=None, P=None, phase=None, time=None, kinetics=None): values, config, original = as_material_array( - material, self._basis, self._phases, self._chemicals + material, self._basis, self._phases, self.chemicals ) - if self._kinetics: self.kinetics(material, time, T, P, phase) + if kinetics or (kinetics is None and self._rate): + if time is None or time == np.inf: + self.equilibrium(material, T, P, phase) + else: + raise NotImplementedError('kinetic integration not yet implemented') conversion = self._conversion(values) if config: material._imol.reset_chemicals(*config) return conversion - def __call__(self, material, time=None, T=None, P=None, phase=None): + def __call__(self, material, T=None, P=None, phase=None, time=None, kinetics=None): values, config, original = as_material_array( - material, self._basis, self._phases, self._chemicals + material, self._basis, self._phases, self.chemicals ) - if self._kinetics: self.kinetics(material, time, T, P, phase) + if kinetics or (kinetics is None and self._rate): + if time is None or time == np.inf: + self.equilibrium(material, T, P, phase) + else: + raise NotImplementedError('kinetic integration not yet implemented') self._reaction(values) if tmo.reaction.CHECK_FEASIBILITY: has_negatives = values.has_negatives() @@ -518,7 +555,7 @@ def __call__(self, material, time=None, T=None, P=None, phase=None): def force_reaction(self, material): """React material ignoring feasibility checks.""" values, config, original = as_material_array( - material, self._basis, self._phases, self._chemicals + material, self._basis, self._phases, self.chemicals ) self._reaction(values) fn.remove_negligible_negative_values(values) @@ -576,12 +613,12 @@ def product_yield(self, product, basis=None, product_yield=None): reactant_index = self._X_index[1] else: reactant_index = self._X_index - product_index = self._chemicals.index(product) + product_index = self.chemicals.index(product) product_coefficient = stoichiometry[product_index] if product_yield is None: product_yield = product_coefficient * self.X if basis and self.basis != basis: - chemicals_tuple = self._chemicals.tuple + chemicals_tuple = self.chemicals.tuple MW_reactant = chemicals_tuple[reactant_index].MW MW_product = chemicals_tuple[product_index].MW if basis == 'wt': @@ -596,7 +633,7 @@ def product_yield(self, product, basis=None, product_yield=None): else: X = product_yield / product_coefficient if basis and self.basis != basis: - chemicals_tuple = self._chemicals.tuple + chemicals_tuple = self.chemicals.tuple MW_reactant = chemicals_tuple[reactant_index].MW MW_product = chemicals_tuple[product_index].MW if basis == 'wt': @@ -763,7 +800,7 @@ def dH(self): * The basis of the reaction (by mol or by wt) """ - chemicals = self._chemicals + chemicals = self.chemicals stoichiometry = self._stoichiometry.to_array() phases = self.phases Hfs = chemicals.Hf @@ -818,12 +855,12 @@ def istoichiometry(self): stoichiometry = self._stoichiometry if stoichiometry.ndim == 1: return tmo.indexer.ChemicalIndexer.from_data(self._stoichiometry, - chemicals=self._chemicals, + chemicals=self.chemicals, check_data=False) else: return tmo.indexer.MaterialIndexer.from_data(self._stoichiometry, phases=self._phases, - chemicals=self._chemicals, + chemicals=self.chemicals, check_data=False) @property @@ -831,14 +868,14 @@ def reactant(self): """[str] Reactant associated to conversion.""" if self._phases: phase_index, chemical_index = self._X_index - return self._phases[phase_index], self._chemicals.IDs[chemical_index] + return self._phases[phase_index], self.chemicals.IDs[chemical_index] else: - return self._chemicals.IDs[self._X_index] + return self.chemicals.IDs[self._X_index] @property def MWs(self): """[1d array] Molecular weights of all chemicals [g/mol].""" - return self._chemicals.MW + return self.chemicals.MW @property def basis(self): @@ -1092,28 +1129,29 @@ class ReactionItem(Reaction): """ __slots__ = ('_index', '_parent') - kinetics = Reaction.kinetics - _kinetics = Reaction._kinetics - phases = MaterialIndexer.phases def __init__(self, rxnset, index): self._stoichiometry = rxnset._stoichiometry[index] self._phases = rxnset._phases self._basis = rxnset._basis self._X = rxnset._X - self._chemicals = rxnset._chemicals + self._thermo = rxnset._thermo self._X_index = rxnset._X_index[index] self._index = index self._parent = rxnset - def reset_chemicals(self, chemicals): - if self._chemicals is chemicals: return + def reset_thermo(self, thermo): + if self._thermo is thermo: return parent = self._parent - parent.reset_chemicals(chemicals) + parent.reset_thermo(thermo) index = self._index self._stoichiometry = parent._stoichiometry[index] self._X_index = parent._X_index[index] - self._chemicals = chemicals + self._thermo = thermo + if hasattr(self, '_stream'): del self._stream + + def reset_chemicals(self, chemicals): + raise Exception('`reset_chemicals` is deprecated; use `reset_thermo` instead') @property def basis(self): @@ -1130,7 +1168,8 @@ def copy(self, basis=None): copy._phases = self._phases copy._stoichiometry = self._stoichiometry.copy() copy._X_index = self._X_index - copy._chemicals = self._chemicals + copy._thermo = self._thermo + copy.chemicals = self.chemicals copy._X = self.X if basis: set_reaction_basis(copy, basis) return copy @@ -1144,6 +1183,7 @@ def X(self, X): self._X[self._index] = X +@thermo_user class ReactionSet: """ Create a ReactionSet that contains all reactions and conversions as an array. @@ -1156,8 +1196,9 @@ class ReactionSet: __slots__ = (*Reaction.__slots__, '_parent_index') copy = Reaction.copy phases = MaterialIndexer.phases - kinetics = Reaction.kinetics - _kinetics = Reaction._kinetics + rate = Reaction.rate + _rate = Reaction._rate + _equilibrium_objective = Reaction._equilibrium_objective _get_stoichiometry_by_mol = Reaction._get_stoichiometry_by_mol _get_stoichiometry_by_wt = Reaction._get_stoichiometry_by_wt force_reaction = Reaction.force_reaction @@ -1171,9 +1212,9 @@ def __init__(self, reactions): if len(phases_set) > 1: raise ValueError('all reactions must implement the same phases') self._phases, = phases_set - chemicals = {i.chemicals for i in reactions} - try: self._chemicals, = chemicals + try: chemicals, = {i.chemicals for i in reactions} except: raise ValueError('all reactions must have the same chemicals') + self._thermo = reactions[0]._thermo basis = {i.basis for i in reactions} try: self._basis, = basis except: raise ValueError('all reactions must have the same basis') @@ -1181,6 +1222,9 @@ def __init__(self, reactions): self._X = np.array([i.X for i in reactions]) X_index = [i._X_index for i in reactions] self._X_index = tuple(X_index) if self._phases else np.array(X_index) + + def equilibrium(self, material, T, P, phase): + raise NotImplementedError('equilibrium of reaction sets not implemented in BioSTEAM (yet)') def __iter__(self): for i in range(self._X.size): yield ReactionItem(self, i) @@ -1194,7 +1238,7 @@ def __getitem__(self, index): rxnset._stoichiometry = stoichiometry rxnset._X = self._X[index] rxnset._X_index = self._X_index[index] - rxnset._chemicals = self._chemicals + rxnset._thermo = self._thermo rxnset._parent_index = (self, index) return rxnset else: @@ -1207,25 +1251,27 @@ def _rescale(self): row = stoichiometry[i] row /= -row[index] - def reset_chemicals(self, chemicals): - if chemicals is self._chemicals: return + def reset_thermo(self, thermo): + if thermo is self._thermo: return if hasattr(self, '_parent_index'): parent, index = self._parent_index - parent.reset_chemicals(chemicals) + parent.reset_thermo(thermo) self._stoichiometry = parent._stoichiometry[index] self._X_index = parent._X_index[index] - self._chemicals = parent._chemicals + self._thermo = thermo + if hasattr(self, '_stream'): del self._stream return phases = self.phases stoichiometry = self._stoichiometry reactants = self.reactants + chemicals = thermo.chemicals if phases: new_stoichiometry = [] for stoic in stoichiometry: A, B = stoic.shape new_stoic = SparseArray.from_shape([A, chemicals.size]) new_stoichiometry.append(new_stoic) - IDs = self._chemicals.IDs + IDs = self.chemicals.IDs for i in range(A): for j in range(B): value = stoic[i, j] @@ -1234,15 +1280,19 @@ def reset_chemicals(self, chemicals): else: A, B = np.array(stoichiometry).shape new_stoichiometry = SparseArray.from_shape([A, chemicals.size]) - IDs = self._chemicals.IDs + IDs = self.chemicals.IDs for i in range(A): for j in range(B): value = stoichiometry[i, j] if value: new_stoichiometry[i, chemicals.index(IDs[j])] = value X_index = [chemicals.index(i) for i in reactants] - self._chemicals = chemicals + self._thermo = thermo self._stoichiometry = new_stoichiometry self._X_index = X_index + if hasattr(self, '_stream'): del self._stream + + def reset_chemicals(self, chemicals): + raise Exception('`reset_chemicals` is deprecated; use `reset_thermo` instead') @property def basis(self): @@ -1262,10 +1312,6 @@ def X(self, X): if X is not self._X: self._X[:] = X @property - def chemicals(self): - """[Chemicals] Chemicals corresponing to each entry in the stoichiometry array.""" - return self._chemicals - @property def stoichiometry(self): """[2d array] Stoichiometry coefficients.""" return self._stoichiometry @@ -1273,7 +1319,7 @@ def stoichiometry(self): @property def reactants(self): """tuple[str] Reactants associated to conversion.""" - IDs = self._chemicals.IDs + IDs = self.chemicals.IDs phases = self._phases X_index = self._X_index if phases: @@ -1284,11 +1330,11 @@ def reactants(self): @property def MWs(self): """[2d array] Molecular weights of all chemicals.""" - return self._chemicals.MW[np.newaxis, :] + return self.chemicals.MW[np.newaxis, :] def to_df(self, index=None): columns = [f'Stoichiometry (by {self.basis})', 'Reactant', 'Conversion [%]'] - chemicals = self._chemicals + chemicals = self.chemicals phases = self._phases rxns = [get_stoichiometric_string(i, phases, chemicals) for i in self._stoichiometry] cmps = [ID + ',' + phase for phase, ID in self.reactants] if phases else self.reactants @@ -1306,7 +1352,7 @@ def __repr__(self): def _info(self, index_name='index'): info = f"{type(self).__name__} (by {self.basis}):" - chemicals = self._chemicals + chemicals = self.chemicals phases = self._phases length = len string = str @@ -1532,6 +1578,8 @@ def X_net(self, indexer=False): else: return X_net + +@thermo_user class ReactionSystem: """ Create a ReactionSystem object that can react a stream across a series of @@ -1616,12 +1664,13 @@ class ReactionSystem: Biomass 0.0437 """ - kinetics = Reaction.kinetics - _kinetics = Reaction._kinetics + rate = Reaction.rate + _rate = Reaction._rate + _equilibrium_objective = Reaction._equilibrium_objective __slots__ = ('_reactions', '_basis', - '_chemicals', + '_thermo', '_phases') def __init__(self, *reactions, basis=None): @@ -1629,17 +1678,20 @@ def __init__(self, *reactions, basis=None): self._reactions = reactions try: self._phases, = set([i._phases for i in reactions]) except: raise ValueError('all reactions must have the same phases') - try: self._chemicals, = set([i._chemicals for i in reactions]) + try: chemicals, = set([i.chemicals for i in reactions]) except: raise ValueError('all reactions must have the same chemicals') + self._thermo = reactions[0]._thermo try: self._basis, = set([i._basis for i in reactions]) except: raise ValueError('all reactions must have the same basis') force_reaction = Reaction.force_reaction adiabatic_reaction = Reaction.adiabatic_reaction - chemicals = Reaction.chemicals __call__ = Reaction.__call__ show = Reaction.show + def equilibrium(self, material, T, P, phase): + raise NotImplementedError('equilibrium of reaction systems not implemented in BioSTEAM (yet)') + def __iter__(self): return iter(self.reactions) @@ -1682,7 +1734,7 @@ def X_net(self, indexer=False): else: X_net[i] = j if indexer: - chemicals = self._chemicals + chemicals = self.chemicals data = chemicals.kwarray(X_net) return ChemicalIndexer.from_data(data, NoPhase, chemicals, False) else: @@ -1703,7 +1755,7 @@ def reactant_flux(self, material, index, subindex=None): """ values, config, original = as_material_array( - material, self._basis, self._phases, self._chemicals + material, self._basis, self._phases, self.chemicals ) preconverted_material = values if original else values.copy() reactions = self.reactions diff --git a/thermosteam/utils/decorators/thermo_user.py b/thermosteam/utils/decorators/thermo_user.py index 35b4a2ad..5b00107e 100644 --- a/thermosteam/utils/decorators/thermo_user.py +++ b/thermosteam/utils/decorators/thermo_user.py @@ -19,7 +19,12 @@ def thermo_user(cls): return cls def _load_thermo(self, thermo): - self._thermo = thermo = tmo.settings.get_default_thermo(thermo) + try: + self._thermo = thermo = tmo.settings.get_default_thermo(thermo) + except ValueError: + self._thermo = thermo = tmo.settings.thermo.subset( + tmo.settings.get_default_chemicals(thermo) + ) return thermo @property