fix issue with EOS-activity coeff vs pressure

Author: yoelcortes

thermosteam/equilibrium/activity_coefficients.py

@@ -131,7 +131,7 @@ def fill_group_psis(group_psis, psis, group_mask):
 @njit(cache=True)
 def gamma_UNIFAC(x, T, interactions, 
                  group_psis, group_mask, qs, rs, Qs,
-                 chemgroups, chem_Qfractions, index):
+                 chemgroups, chem_Qfractions, index, *, P=None):
     N_chemicals = index.size
     gamma = np.ones(x.size)
     if N_chemicals > 1:
@@ -157,7 +157,7 @@ def gamma_UNIFAC(x, T, interactions,
 @njit(cache=True)
 def gamma_modified_UNIFAC(x, T, interactions, 
                    group_psis, group_mask, qs, rs, Qs,
-                   chemgroups, chem_Qfractions, index):
+                   chemgroups, chem_Qfractions, index, *, P=None):
     N_chemicals = index.size
     gamma = np.ones(x.size)
     if N_chemicals > 1:
@@ -225,11 +225,12 @@ class IdealActivityCoefficients(ActivityCoefficients):
     
     """
     __slots__ = ()
+    P_dependent = False
 
     def __init__(self, chemicals):
         self._chemicals = tuple(chemicals)
 
-    def __call__(self, xs, T):
+    def __call__(self, xs, T, P=101325):
         return np.ones(len(xs))
 
 
@@ -244,6 +245,7 @@ class GroupActivityCoefficients(ActivityCoefficients):
     chemicals : Iterable[:class:`~thermosteam.Chemical`]
     
     """
+    P_dependent = False
     __slots__ = ('_rs', '_qs', '_Qs','_chemgroups',
                  '_group_psis',  '_chem_Qfractions',
                  '_group_mask', '_interactions',
@@ -304,7 +306,7 @@ def args(self):
                 self._chemgroups, self._chem_Qfractions, 
                 self._index)
 
-    def activity_coefficients(self, x, T):
+    def activity_coefficients(self, x, T, P=101325):
         """
         Return activity coefficients of chemicals with defined functional groups.
         
@@ -324,7 +326,7 @@ def activity_coefficients(self, x, T):
                                            self._chem_Qfractions,
                                            self._group_psis)
 
-    def __call__(self, x, T):
+    def __call__(self, x, T, P=101325):
         """
         Return activity coefficients.
         
@@ -485,6 +487,7 @@ class GCEOSActivityCoefficients(ActivityCoefficients):
     
     """
     __slots__ = ('_chemicals', '_eos')
+    P_dependent = True
     EOS = None # type[GCEOSMIX] Subclasses must implement this attribute.
     chemsep_db = None # Optional[str] Name of chemsep data base for interaction parameters.
     cache = None # [dict] Subclasses must implement this attribute.
@@ -553,6 +556,7 @@ def __call__(self, x, T, P=101325):
     activity_coefficient_classes.append(cls)
     dct[clsname] = cls
 
-dct['PRActivityCoefficients'].chemsep_db = 'ChemSep PR'
+for i in ('PR', 'PR78'):
+    dct[i + 'ActivityCoefficients'].chemsep_db = 'ChemSep PR'
 __all__ = (*__all__, *clsnames, 'activity_coefficient_classes')
 del dct, clsnames

thermosteam/equilibrium/bubble_point.py

@@ -146,7 +146,7 @@ def _T_error(self, T, P, z_over_P, z_norm, y):
         Psats = np.array([i(T) for i in self.Psats], dtype=float)
         y_phi =  (z_over_P
                   * Psats
-                  * self.gamma(z_norm, T) 
+                  * self.gamma(z_norm, T, P) 
                   * self.pcf(T, P, Psats))
         y[:] = solve_y(y_phi, self.phi, T, P, y)
         return 1. - y.sum()
@@ -157,6 +157,12 @@ 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 _P_error_dep(self, P, T, z, Psats, z_Psats, y):
+        if P <= 0: raise InfeasibleRegion('negative pressure')
+        y_phi = z_Psats * self.gamma(z, T, P) * self.pcf(T, P, Psats) / P
+        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_conversion):
         if T <= 0: raise InfeasibleRegion('negative temperature')
         dz[:] = liquid_conversion(z, T, P, 'l')
@@ -165,7 +171,7 @@ def _T_error_reactive(self, T, P, z, dz, y, x, liquid_conversion):
         Psats = np.array([i(T) for i in self.Psats], dtype=float)
         y_phi =  (x / P
                   * Psats
-                  * self.gamma(x, T) 
+                  * self.gamma(x, T, P) 
                   * self.pcf(T, P, Psats))
         y[:] = solve_y(y_phi, self.phi, T, P, y)
         return 1. - y.sum()
@@ -175,7 +181,7 @@ def _P_error_reactive(self, P, T, Psats, z, dz, y, x, liquid_conversion):
         dz[:] = liquid_conversion(z, T, P, 'l')
         x[:] = z + dz
         x /= x.sum()
-        z_Psat_gamma = x * Psats * self.gamma(x, T)
+        z_Psat_gamma = x * Psats * self.gamma(x, T, P)
         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()
@@ -341,10 +347,16 @@ def solve_Py(self, z, T, liquid_conversion=None):
             elif T < self.Tmin: T = self.Tmin
             Psats = np.array([i(T) for i in self.Psats])
             z_norm = z / z.sum()
-            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, Psats, y)
+            if self.gamma.P_dependent:
+                f = self._P_error_dep
+                z_Psats = z_norm * Psats
+                P_guess, y = self._Py_ideal(z_Psats)
+                args = (T, z_norm, Psats, z_Psats, y)
+            else:
+                z_Psat_gamma = z_norm * Psats * self.gamma(z_norm, T, 101325)
+                P_guess, y = self._Py_ideal(z_Psat_gamma)
+                f = self._P_error
+                args = (T, z_Psat_gamma, Psats, y)
             try:
                 P = flx.aitken_secant(f, P_guess, P_guess-1, self.P_tol, 1e-9,
                                       args, checkiter=False, maxiter=self.maxiter)
@@ -362,7 +374,7 @@ def solve_Py(self, z, T, liquid_conversion=None):
             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)
+            z_Psat_gamma = z * Psats * self.gamma(z_norm, T, 101325)
             P_guess, y = self._Py_ideal(z_Psat_gamma)
             args = (T, Psats, z_norm, dz, y, x, liquid_conversion)
             try:

thermosteam/equilibrium/dew_point.py

@@ -30,13 +30,13 @@ def gamma_iter(gamma, x_gamma, T, P, f_gamma, gamma_args):
     mask = x < 1e-32
     x[mask] = 1e-32
     x = fn.normalize(x)
-    return f_gamma(x, T, *gamma_args)
+    return f_gamma(x, T, *gamma_args, P=P)
 
 def solve_x(x_guess, x_gamma, T, P, f_gamma, gamma_args):
     mask = x_guess < 1e-32
     x_guess[mask] = 1e-32
     x_guess = fn.normalize(x_guess)
-    gamma = f_gamma(x_guess, T, *gamma_args)
+    gamma = f_gamma(x_guess, T, *gamma_args, P=P)
     args = (x_gamma, T, P, f_gamma, gamma_args)
     gamma = flx.wegstein(
         gamma_iter, gamma, 1e-12, args=args, checkiter=False,

thermosteam/equilibrium/fugacities.py

@@ -53,7 +53,7 @@ def __init__(self, chemicals, thermo=None):
 
     def unweighted(self, x, T, P=101325.):
         Psats = np.array([i.Psat(T) for i in self.chemicals], dtype=float)
-        return self.gamma(x, T) * self.pcf(T, P, Psats) * Psats
+        return self.gamma(x, T, P) * self.pcf(T, P, Psats) * Psats
 
     def __call__(self, x, T, P=101325., reduce=False):
         f_reduced = x * self.gamma(x, T)

thermosteam/equilibrium/vle.py

@@ -87,7 +87,7 @@ def xVlogK_iter(
     x = xVlogK[:n]
     Ks = np.exp(xVlogK[n+1:])
     x, y = xy(x, Ks)
-    Ks[:] = pcf_Psat_over_P * f_gamma(x, T, *gamma_args) / f_phi(y, T, P)
+    Ks[:] = pcf_Psat_over_P * f_gamma(x, T, *gamma_args, P=P) / f_phi(y, T, P)
     V = xVlogK[n]
     if V < 0.: V = 0.
     elif V > 1.: V = 1.