better procedure for saturation lines, and python-API

c_interface/yaeos_c.f90

@@ -44,6 +44,7 @@ module yaeos_c
    ! Phase equilibria
    public :: flash, flash_grid
    public :: saturation_pressure, saturation_temperature
+   public :: pure_saturation_line
    public :: pt2_phase_envelope, px2_phase_envelope
    public :: critical_point, critical_line
    public :: stability_zpt, tm
@@ -707,6 +708,40 @@ subroutine saturation_temperature(id, z, P, kind, T0, T, x, y, Vx, Vy, beta)
       call equilibria_state_to_arrays(sat, x, y, aux, T, Vx, Vy, beta)
    end subroutine saturation_temperature
 
+   subroutine pure_saturation_line(id, comp_id, stop_P, stop_T, P, T, Vx, Vy)
+      use yaeos, only: fsat => pure_saturation_line, PurePsat, pr
+      integer(c_int), intent(in) :: id
+      integer(c_int), intent(in) :: comp_id
+      real(c_double), intent(in) :: stop_P
+      real(c_double), intent(in) :: stop_T
+      real(c_double), intent(out) :: P(800)
+      real(c_double), intent(out) :: T(800)
+      real(c_double), intent(out) :: Vx(800)
+      real(c_double), intent(out) :: Vy(800)
+
+      integer :: npoints
+      type(PurePsat) :: sat
+
+      real(8) :: nan
+
+      nan = 0
+      nan = nan/nan
+
+      T = nan
+      P = nan
+      Vx = nan
+      Vy = nan
+
+      sat = fsat(ar_models(id)%model, comp_id, stop_P, stop_T)
+
+      npoints = minval([size(sat%T), 800])
+
+      T(:npoints) = sat%T(:npoints)
+      P(:npoints) = sat%P(:npoints)
+      Vx(:npoints) = sat%Vx(:npoints)
+      Vy(:npoints) = sat%Vy(:npoints)
+   end subroutine
+
    subroutine pt2_phase_envelope(id, z, kind, max_points, Ts, Ps, tcs, pcs, T0, P0)
       use yaeos, only: &
          saturation_pressure, saturation_temperature, pt_envelope_2ph, &

python/yaeos/__init__.py

@@ -4,8 +4,9 @@
 relevant constants, procedures and objects to have better access to them.
 """
 
-__version__ = "2.1.0"
+__version__ = "2.3.0"
 
+import yaeos.constants as constants
 from yaeos.lib import yaeos_c
 from yaeos.models.excess_gibbs import NRTL, UNIFACVLE, UNIQUAC
 from yaeos.models.residual_helmholtz.cubic_eos import (
@@ -22,6 +23,7 @@
 
 
 __all__ = [
+    "constants",
     "yaeos_c",
     "SoaveRedlichKwong",
     "PengRobinson76",

python/yaeos/core.py

@@ -978,6 +978,52 @@ def cp_residual_vt(
         """
         return yaeos_c.cp_residual_vt(self.id, moles, volume, temperature)
 
+    def pure_saturation_pressures(
+        self, component, stop_pressure=0.01, stop_temperature=100
+    ):
+        """Calculate pure component saturation pressures [bar].
+
+        Calculation starts from the critical point and goes down to the
+        stop pressure or stop temperature.
+
+        Parameters
+        ----------
+        component : int
+            Component index (starting from 1)
+        stop_pressure : float, optional
+            Stop pressure [bar], by default 0.01
+        stop_temperature : float, optional
+            Stop temperature [K], by default 100
+
+        Returns
+        -------
+        dict
+            Pure component saturation points dictionary with keys:
+                - T: Temperature [K]
+                - P: Pressure [bar]
+                - Vx: Liquid Phase Volume [L/mole]
+                - Vy: Vapor Phase Volume [L/mole]
+
+        Example
+        -------
+        .. code-block:: python
+
+            import numpy as np
+
+            from yaeos import PengRobinson76
+
+            tc = np.array([320.0, 375.0])
+            pc = np.array([45.0, 60.0])
+            w = np.array([0.0123, 0.045])
+
+            model = PengRobinson76(tc, pc, w)
+        """
+        P, T, Vx, Vy = yaeos_c.pure_saturation_line(
+            self.id, component, stop_pressure, stop_temperature
+        )
+
+        return {"T": T, "P": P, "Vx": Vx, "Vy": Vy}
+
     def flash_pt(
         self, z, pressure: float, temperature: float, k0=None
     ) -> dict:

src/equilibria/boundaries/pure_saturation.f90

@@ -22,7 +22,7 @@ module yaeos__equilibria_boundaries_pure_saturation
 
 contains
 
-   function pure_saturation_line(model, component) result(pt)
+   function pure_saturation_line(model, component, minP, minT) result(pt)
       !! # Pure saturation line
       !!
       !! Saturation pressures and temperatures for a pure component.
@@ -39,6 +39,8 @@ function pure_saturation_line(model, component) result(pt)
       use stdlib_optval, only: optval
       class(ArModel), intent(in) :: model !! Thermodyanmic model
       integer, intent(in) :: component !! Component index to calculate the line
+      real(pr), intent(in) :: minP !! Minimum pressure [bar]
+      real(pr), intent(in) :: minT !! Minimum temperature [K]
       type(PurePsat) :: pt
       ! ------------------------------------------------------------------------
 
@@ -64,49 +66,56 @@ function pure_saturation_line(model, component) result(pt)
       z(component) = 1
       call model%volume(z, P=Pc, T=Tc, V=Vc, root_type="vapor")
 
-      Vx = Vc*0.99
-      Vy = Vc*1.1
+      Vx = Vc*0.999
+      Vy = Vc*1.001
 
       X = [log(Vx), log(Vy), log(Pc), log(Tc)]
 
       ns = 1
-      S = log(0.95)
-      dS = -0.5
+      S = log(0.99)
+      dS = -0.15
       allocate(pt%T(0), pt%P(0), pt%Vx(0), pt%Vy(0))
 
       ! ========================================================================
       ! Trace the line using the continuation method.
       ! ------------------------------------------------------------------------
       T = Tc
+      P = Pc
       points = 0
-      do while(T > 50 .or. .not. isnan(T))
+      do while(T > minT .and. P > minP .and. .not. isnan(T))
          call solve_point(model, component, nc, X, ns, S, F, dF, dFdS, its)
          dXdS = solve_system(dF, -dFdS)
          ns = maxloc(abs(dXdS(3:4)), dim=1) + 2
          dS = dXdS(ns)*dS
          dXdS = dXdS/dXdS(ns)
-         X = X + dXdS*dS
-         S = X(ns)
+
+         do while (exp(X(4)) - exp(X(4) + dXdS(4)*dS) < 3 .and. ((Tc - T) > 10 .or. (Pc - P) > 2))
+            dS = dS*1.5
+         end do
 
          Vx = exp(X(1))
          Vy = exp(X(2))
-         P = exp(X(3))
-         T = exp(X(4))
+         P  = exp(X(3))
+         T  = exp(X(4))
 
-         if (.not. isnan(T)) then
+         if (isnan(T)) then
+            exit
+         else
             pt%T = [pt%T, T]
             pt%P = [pt%P, P]
             pt%Vx = [pt%Vx, Vx]
             pt%Vy = [pt%Vy, Vy]
             points = points + 1
          end if
+
+         X = X + dXdS*dS
+         S = X(ns)
       end do
-      
+
       ! Save interpolators to obtain particular values. The interpolator needs
       ! monothonic increasing values in x, so we need to reverse the arrays.
-
-      pt%P = pt%P(points:1:-1)
-      pt%T = pt%T(points:1:-1)
+      pt%P  = pt%P(points:1:-1)
+      pt%T  = pt%T(points:1:-1)
       pt%Vx = pt%Vx(points:1:-1)
       pt%Vy = pt%Vy(points:1:-1)
 
@@ -144,6 +153,7 @@ subroutine solve_point(model, ncomp, nc, X, ns, S, F, dF, dFdS, its)
       !!
       !! The vector of variables \(X\) is equal to
       !! \([ \ln V_z, \ln V_y, \ln P, \ln T ]\).
+!
       class(ArModel), intent(in) :: model
       !! Thermodynamic model
       integer, intent(in) :: ncomp
@@ -178,7 +188,8 @@ subroutine solve_point(model, ncomp, nc, X, ns, S, F, dF, dFdS, its)
       real(pr) :: dX(4), B
       real(pr) :: Xnew(4)
 
-      integer :: i, j
+      integer :: i
+!
       i = ncomp
 
       dX = 1
@@ -188,9 +199,10 @@ subroutine solve_point(model, ncomp, nc, X, ns, S, F, dF, dFdS, its)
       B = model%get_v0(z, 1._pr, 150._pr)
 
       its = 0
-      do while((maxval(abs(dX)) > 1e-7 .or. maxval(abs(F)) > 1e-7) .and. j < 100 )
+      do while((maxval(abs(dX)) > 1e-7 .and. maxval(abs(F)) > 1e-7))
          its = its+1
          call isofugacity(X, F, dF, dFdS)
+         if (any(isnan(F))) exit
          dX = solve_system(dF, -F)
          Xnew = X + dX
          X = Xnew

test/test_purePsat.f90

@@ -31,15 +31,15 @@ program main
    ! ==========================================================================
    ! Test the pure saturation line of propane
    ! -------------------------------------------------------------------------
-   pt = pure_saturation_line(model, 2)
+   pt = pure_saturation_line(model, 2, 0.001_pr, 100._pr)
    call assert(abs(pt%get_P(200._pr) - 0.2068) < 0.1_Pr, "Propane Psat at 140K")
    call assert(abs(pt%get_T(10._pr) - 300.08) < 0.1_Pr, "Propane Psat at 10 bar")
 
 
    ! ==========================================================================
    ! Test the pure saturation line of methane
    ! -------------------------------------------------------------------------
-   pt = pure_saturation_line(model, 1)
+   pt = pure_saturation_line(model, 1, 1._pr, 100._pr)
    call assert(abs(pt%get_P(140._pr) - 6.45) < 0.1_Pr, "Methane Psat at 140K")
-   call assert(abs(pt%get_T(10._pr) - 148.970) < 0.1_Pr, "Propane Psat at 10 bar")
+   call assert(abs(pt%get_T(10._pr) - 148.970) < 0.1_Pr, "Methane Psat at 10 bar")
 end program main

test/test_tx.f90

@@ -35,8 +35,8 @@ program main
       w=sus%critical%acentric_factor%value &
       )
 
-   vp1 = pure_saturation_line(model, 1)
-   vp2 = pure_saturation_line(model, 2)
+   vp1 = pure_saturation_line(model, 1, 1._pr, 100._pr)
+   vp2 = pure_saturation_line(model, 2, 1._pr, 100._pr)
 
    P = 10._pr
    T = vp1%get_T(P)