diff --git a/thermosteam/equilibrium/activity_coefficients.py b/thermosteam/equilibrium/activity_coefficients.py index 05790458..66cce4ff 100644 --- a/thermosteam/equilibrium/activity_coefficients.py +++ b/thermosteam/equilibrium/activity_coefficients.py @@ -16,6 +16,7 @@ """ """ import numpy as np +from warnings import warn from .unifac import DOUFSG, DOUFIP2016, UFIP, UFSG from flexsolve import njitable @@ -65,12 +66,17 @@ def get_interaction(all_interactions, i, j, no_interaction): def get_chemgroups(chemicals, field): getfield = getattr chemgroups = [] - for chemical in chemicals: + index = [] + for i, chemical in enumerate(chemicals): group = getfield(chemical, field) - if not group: - raise RuntimeError(f"{chemical} has no defined {field} UNIFAC groups") - chemgroups.append(group) - return chemgroups + if group: + chemgroups.append(group) + index.append(i) + else: + warn(f"{chemical} has no defined {field} groups; " + "functional group interactions are ignored", + RuntimeWarning, stacklevel=3) + return index, chemgroups @njitable(cache=True) def loggammacs_UNIFAC(qs, rs, x): @@ -157,7 +163,7 @@ class GroupActivityCoefficients(ActivityCoefficients): __slots__ = ('_rs', '_qs', '_Qs','_chemgroups', '_group_psis', '_chem_Qfractions', '_group_mask', '_interactions', - '_chemicals') + '_chemicals', '_index') def __new__(cls, chemicals): chemicals = tuple(chemicals) @@ -165,7 +171,8 @@ def __new__(cls, chemicals): return cls._cached[chemicals] else: self = super().__new__(cls) - chemgroups = get_chemgroups(chemicals, self.group_name) + index, chemgroups = get_chemgroups(chemicals, self.group_name) + self._index = None if len(index) == len(chemicals) else index all_groups = set() for groups in chemgroups: all_groups.update(groups) index = {group:i for i,group in enumerate(all_groups)} @@ -204,31 +211,60 @@ def __new__(cls, chemicals): def __reduce__(self): return type(self), (self.chemicals,) - def __call__(self, x, T): - """Return UNIFAC coefficients. + def activity_coefficients(self, x, T): + """ + Return activity coefficients of chemicals with defined functional groups. Parameters ---------- x : array_like Molar fractions T : float - Temperature (K) + Temperature [K] """ - x = np.asarray(x) psis = self.psi(T, self._interactions.copy()) self._group_psis[self._group_mask] = psis[self._group_mask] - gamma = group_activity_coefficients(x, self._chemgroups, - self.loggammacs(self._qs, self._rs, x), - self._Qs, psis, - self._chem_Qfractions, - self._group_psis) + return group_activity_coefficients(x, self._chemgroups, + self.loggammacs(self._qs, self._rs, x), + self._Qs, psis, + self._chem_Qfractions, + self._group_psis) + + def __call__(self, x, T): + """ + Return activity coefficients. + + Parameters + ---------- + x : array_like + Molar fractions + T : float + Temperature [K] + + Notes + ----- + Activity coefficients of chemicals with missing groups are default to 1. + + """ + x = np.asarray(x, float) + if self._index: + N_chemicals = x.size + x = x[self._index] + xsum = x.sum() + gamma = np.ones(N_chemicals) + if xsum: gamma[self._index] = self.activity_coefficients(x / xsum, T) + else: + gamma = self.activity_coefficients(x, T) gamma[np.isnan(gamma)] = 1 return gamma class UNIFACActivityCoefficients(GroupActivityCoefficients): - """Create a UNIFACActivityCoefficients that estimates activity coefficients using the UNIFAC group contribution method when called with a composition and a temperature (K). + """ + Create a UNIFACActivityCoefficients that estimates activity coefficients + using the UNIFAC group contribution method when called with a composition + and a temperature (K). Parameters ---------- @@ -251,13 +287,41 @@ def psi(T, a): class DortmundActivityCoefficients(GroupActivityCoefficients): - """Create a DortmundActivityCoefficients that estimates activity coefficients using the Dortmund UNIFAC group contribution method when called with a composition and a temperature (K). + """ + Create a DortmundActivityCoefficients that estimates activity coefficients + using the Dortmund UNIFAC group contribution method when called with a + composition and a temperature (K). Parameters ---------- chemicals : Iterable[Chemical] + Examples + -------- + >>> import thermosteam as tmo + >>> chemicals = tmo.Chemicals(['Water', 'Ethanol'], cache=True) + >>> Gamma = tmo.equilibrium.DortmundActivityCoefficients(chemicals) + >>> composition = [0.5, 0.5] + >>> T = 350. + >>> Gamma(composition, T) + array([1.475, 1.242]) + + >>> chemicals = tmo.Chemicals(['Dodecane', 'Tridecane'], cache=True) + >>> Gamma = tmo.equilibrium.DortmundActivityCoefficients(chemicals) + >>> # Note how both hydrocarbons have similar lenghts and structure, + >>> # so activities should be very close + >>> Gamma([0.5, 0.5], 350.) + array([1., 1.]) + + >>> chemicals = tmo.Chemicals(['Water', 'O2'], cache=True) + >>> Gamma = tmo.equilibrium.DortmundActivityCoefficients(chemicals) + >>> # The following warning is issued because O2 does not have Dortmund groups + >>> # RuntimeWarning: O2 has no defined Dortmund groups; + >>> # functional group interactions are ignored + >>> Gamma([0.5, 0.5], 350.) + array([1., 1.]) + """ __slots__ = () all_subgroups = DOUFSG diff --git a/thermosteam/equilibrium/bubble_point.py b/thermosteam/equilibrium/bubble_point.py index 1714ec3e..f11a4b98 100644 --- a/thermosteam/equilibrium/bubble_point.py +++ b/thermosteam/equilibrium/bubble_point.py @@ -66,7 +66,7 @@ class BubblePoint: """ __slots__ = ('chemicals', 'IDs', 'gamma', 'phi', 'pcf', 'P', 'T', 'y', 'Psats', 'Tmin', 'Tmax', 'Pmin', 'Pmax') - Tmin_default = 200. + Tmin_default = 150. _cached = {} def __init__(self, chemicals=(), thermo=None): thermo = settings.get_default_thermo(thermo) diff --git a/thermosteam/equilibrium/dew_point.py b/thermosteam/equilibrium/dew_point.py index 0f43feda..da95115d 100644 --- a/thermosteam/equilibrium/dew_point.py +++ b/thermosteam/equilibrium/dew_point.py @@ -68,7 +68,7 @@ class DewPoint: __slots__ = ('chemicals', 'phi', 'gamma', 'IDs', 'pcf', 'Psats', 'P', 'T', 'x', 'Tmin', 'Tmax', 'Pmin', 'Pmax') - Tmin_default = 200. + Tmin_default = 150. _cached = {} def __init__(self, chemicals=(), thermo=None): thermo = settings.get_default_thermo(thermo)