Merge pull request #128 from ipqa-research/dev

Dev

c_interface/yaeos_c.f90

@@ -21,7 +21,7 @@ module yaeos_c
    ! CubicEoS
    public :: srk, pr76, pr78, rkpr, psrk
    ! Mixing rules
-   public :: set_mhv, set_qmr, set_hv
+   public :: set_mhv, set_qmr, set_qmrtd, set_hv
 
    ! __del__
    public :: make_available_ar_models_list
@@ -193,6 +193,27 @@ end subroutine make_available_ar_models_list
    ! ==========================================================================
    !  Cubic Mixing rules
    ! --------------------------------------------------------------------------
+   subroutine set_qmrtd(ar_id, kij_0, kij_inf, t_star, lij)
+      use yaeos, only: QMRTD, CubicEoS
+      integer(c_int), intent(in) :: ar_id
+      real(c_double), intent(in) :: kij_0(:, :)
+      real(c_double), intent(in) :: kij_inf(:, :)
+      real(c_double), intent(in) :: t_star(:, :)
+      real(c_double), intent(in) :: lij(:, :)
+
+      type(QMRTD) :: mixrule
+
+      mixrule = QMRTD(k=kij_inf, k0=kij_0, Tref=t_star, l=lij)
+
+      associate (ar_model => ar_models(ar_id)%model)
+         select type(ar_model)
+          class is(CubicEoS)
+            deallocate(ar_model%mixrule)
+            ar_model%mixrule = mixrule
+         end select
+      end associate
+   end subroutine set_qmrtd
+
    subroutine set_mhv(ar_id, ge_id, q)
       !! Michelsen's Modified Huron-Vidal 1 with constant `q_1` parameter
       use yaeos, only: MHV, CubicEoS
@@ -329,7 +350,7 @@ subroutine residual_helmholtz(id, n, v, t, ar, ArT, ArV, ArTV, ArV2, ArT2, Arn,
          n=n, V=V, T=T, &
          Ar=Ar,  ArV=ArV, ArT=ArT, ArTV=ArTV, &
          ArV2=ARV2, ArT2=ArT2, Arn=Arn, ArVn=ArVn, ArTn=ArTn, Arn2=Arn2)
-   end subroutine
+   end subroutine residual_helmholtz
 
    subroutine lnphi_vt(id, n, v, t, lnphi, dlnphidp, dlnphidt, dlnphidn)
       integer(c_int), intent(in) :: id
@@ -717,5 +738,5 @@ function makenan()
       real(c_double) :: makenan
       makenan = 0
       makenan = makenan/makenan
-   end function
+   end function makenan
 end module yaeos_c

python/yaeos/__init__.py

@@ -15,6 +15,7 @@
     PengRobinson76,
     PengRobinson78,
     QMR,
+    QMRTD,
     RKPR,
     SoaveRedlichKwong,
 )
@@ -28,6 +29,7 @@
     "RKPR",
     "PSRK",
     "QMR",
+    "QMRTD",
     "NRTL",
     "UNIFACVLE",
     "UNIQUAC",

python/yaeos/core.py

@@ -1106,7 +1106,7 @@ def critical_point(self, z, max_iters=100) -> dict:
             self.id, z0=z, zi=[0, 0], spec=1, max_iters=max_iters
         )
 
-        return {"Tc": t, "Pc": p, "Vc": v}
+        return {"x": x, "Tc": t, "Pc": p, "Vc": v}
 
     def critical_line(self, z0, zi, a0=1e-5, da0=1e-2, max_points=1000):
         """Critical Line calculation.

python/yaeos/models/residual_helmholtz/cubic_eos/__init__.py

@@ -20,7 +20,7 @@
     RKPR,
     SoaveRedlichKwong,
 )
-from .mixing_rules import CubicMixRule, HV, MHV, QMR
+from .mixing_rules import CubicMixRule, HV, MHV, QMR, QMRTD
 
 
 __all__ = [
@@ -32,6 +32,7 @@
     "PSRK",
     "CubicMixRule",
     "QMR",
+    "QMRTD",
     "MHV",
     "HV",
 ]

python/yaeos/models/residual_helmholtz/cubic_eos/mixing_rules.py

@@ -80,6 +80,72 @@ def set_mixrule(self, ar_model_id: int) -> None:
         yaeos_c.set_qmr(ar_model_id, self.kij, self.lij)
 
 
+class QMRTD(CubicMixRule):
+    r"""Quadratic mixing rule, with temperature dependence.
+
+    ..math::
+        k_{ij}(T) = k_{ij}^{\infty} + k_{ij}^0 \exp{\left(-T/T^{ref}\right)}
+
+    Parameters
+    ----------
+    kij_0 : array_like
+        kij_0 binary interaction parameters matrix
+    kij_inf : array_like
+        kij_inf binary interaction parameters matrix
+    t_ref: array_like
+        Reference temperature
+    lij : array_like
+        lij binary interaction parameters matrix
+
+    Attributes
+    ----------
+    kij_0 : array_like
+        kij_0 binary interaction parameters matrix
+    kij_inf : array_like
+        kij_inf binary interaction parameters matrix
+    t_ref: array_like
+        Reference temperature
+    lij : array_like
+        lij binary interaction parameters matrix
+
+    Example
+    -------
+    .. code-block:: python
+
+        from yaeos import QMRTD, SoaveRedlichKwong
+
+        kij_0 = [[0.0, 0.1], [0.1, 0.0]]
+        kij_inf = [[0.0, 0.1], [0.1, 0.0]]
+        Tref = [[0.0, 390], [390, 0.0]]
+        lij = [[0.0, 0.02], [0.02, 0.0]]
+
+        # Quadratic mixing rule instance
+        mixrule = QMRTD(kij_0, kij_inf, Tref, lij)
+
+        tc = [305.32, 469.7]        # critical temperature [K]
+        pc = [48.72, 33.7]          # critical pressure [bar]
+        w = [0.0995, 0.152]         # acentric factor
+
+        model = SoaveRedlichKwong(tc, pc, w, mixrule)
+    """
+
+    def __init__(self, kij_0, kij_inf, t_ref, lij) -> None:
+        self.kij_0 = np.array(kij_0, order="F")
+        self.kij_inf = np.array(kij_inf, order="F")
+        self.t_ref = np.array(t_ref, order="F")
+        self.lij = np.array(lij, order="F")
+
+    def set_mixrule(self, ar_model_id: int) -> None:
+        """Set mix rule method."""
+        yaeos_c.set_qmrtd(
+            ar_model_id,
+            kij_0=self.kij_0,
+            kij_inf=self.kij_inf,
+            t_star=self.t_ref,
+            lij=self.lij,
+        )
+
+
 class MHV(CubicMixRule):
     """Modified Huron-Vidal mixing rule.