tests and PR76

src/models/residual_helmholtz/cubic/implementations/pr76.f90

@@ -0,0 +1,55 @@
+module yaeos_models_ar_cubic_pengrobinson76
+    use yaeos_constants, only: R, pr
+    use yaeos_substance, only: Substances
+    use yaeos_models_ar_genericcubic, only: CubicEoS, AlphaFunction, CubicMixRule
+    use yaeos_models_ar_genericcubic_quadratic_mixing, only: QMR
+    use yaeos_models_ar_cubic_alphas, only: AlphaSoave
+    implicit none
+
+    private
+
+    public :: PengRobinson76
+
+contains
+
+    type(CubicEoS) function PengRobinson76(tc, pc, w, kij, lij) result(model)
+        real(pr), intent(in) :: tc(:), pc(:), w(:)
+        real(pr), optional, intent(in) :: kij(:, :), lij(:, :)
+
+        type(Substances) :: composition
+        type(QMR) :: mixrule
+        type(AlphaSoave) :: alpha
+        integer :: nc
+        integer :: i
+
+        nc = size(tc)
+
+        composition%tc = tc
+        composition%pc = pc
+        composition%w = w
+
+        alpha%k = 0.37464_pr &
+                 + 1.54226_pr * composition%w &
+                 - 0.26993_pr * composition%w**2
+
+        if (present(kij)) then
+            mixrule%k = kij
+        else
+            mixrule%k = reshape([(i, i=1,nc**2)], [nc, nc])
+        endif
+
+        if (present(lij)) then
+            mixrule%l = lij
+        else
+            mixrule%l = reshape([(i, i=1,nc**2)], [nc, nc])
+        endif
+
+        model%components = composition
+        model%ac = 0.45723553_pr * R**2 * composition%tc**2 / composition%pc
+        model%b = 0.07779607_pr * R * composition%tc/composition%pc
+        model%del1 = [1 + sqrt(2.0_pr)]
+        model%del2 = [1 - sqrt(2.0_pr)]
+        model%alpha = alpha
+        model%mixrule = mixrule
+    end function
+end module

test/fixtures/fix_models.f90

@@ -0,0 +1,23 @@
+module fixtures_models
+    use yaeos, only: pr, R, CubicEoS
+
+contains
+
+    type(CubicEoS) function binary_PR76() result(eos)
+        use yaeos, only: PengRobinson76
+        integer, parameter :: n=2
+        real(pr) :: tc(n), pc(n), w(n)
+        real(pr) :: kij(n, n), lij(n, n)
+        ! z = [0.3_pr, 0.7_pr]
+        tc = [190._pr, 310._pr]
+        pc = [14._pr, 30._pr]
+        w = [0.001_pr, 0.03_pr]
+
+        kij = reshape([0., 0.1, 0.1, 0.], [n,n]) 
+        lij = kij / 2 
+
+
+        eos = PengRobinson76(tc, pc, w, kij, lij)
+    end function
+
+end module

test/test_thermoprops.f90

@@ -0,0 +1,105 @@
+module test_thermoprops
+   use testdrive, only: new_unittest, unittest_type, error_type, check
+   implicit none
+
+contains
+   subroutine collect_suite(testsuite)
+      !> Collection of tests
+      type(unittest_type), allocatable, intent(out) :: testsuite(:)
+
+      testsuite = [ &
+                  new_unittest("Fugacity VT", test_fugacity_VT), &
+                  new_unittest("Fugacity tp", test_fugacity_tp) &
+                  ]
+   end subroutine
+
+   subroutine test_fugacity_VT(error)
+      use fixtures_models, only: binary_PR76
+      use yaeos, only: pr, CubicEoS, fugacity_vt
+      type(error_type), allocatable, intent(out) :: error
+      type(CubicEoS) :: eos
+
+      real(pr) :: lnfug(2), dlnphidp(2), dlnphidt(2), dlnphidn(2,2)
+
+      real(pr), allocatable :: z(:)
+      real(pr) ::  v, t, p
+
+      real(pr) :: lnfug_val(2), dlnphidp_val(2), dlnphidt_val(2)
+
+      lnfug_val =   [2.0759140949373416, -2.2851989270402058]
+      dlnphidp_val = [  -0.99059224575177762, -0.99388122357848807]
+      dlnphidt_val = [ 3.0263769083149254E-002,  7.6204871541712640E-002]
+
+      eos = binary_PR76()
+
+      z = [0.3_pr, 0.7_pr]
+      v = 8.8780451065729321E-002_pr
+      t = 150
+
+
+      call fugacity_vt(eos, &
+         z, V, T, P, lnfug, dlnPhidP, dlnphidT, dlnPhidn &
+      )
+
+      call check(&
+         error, maxval(abs(lnfug - lnfug_val)) < 1e-5 &
+      )
+      call check(&
+         error, maxval(abs(dlnphidp - dlnphidp_val)) < 1e-5 &
+      )
+      call check(&
+         error, maxval(abs(dlnphidt - dlnphidt_val)) < 1e-5 &
+      )
+   end subroutine
+
+   subroutine test_fugacity_TP(error)
+      use fixtures_models, only: binary_PR76
+      use yaeos, only: pr, CubicEoS, fugacity_tp
+      type(error_type), allocatable, intent(out) :: error
+      type(CubicEoS) :: eos
+
+      real(pr) :: lnfug(2), dlnphidp(2), dlnphidt(2), dlnphidn(2,2)
+
+      real(pr), allocatable :: z(:)
+      real(pr) ::  v, t, p
+
+      real(pr) :: lnfug_val(2), dlnphidp_val(2), dlnphidt_val(2)
+
+      integer :: root_type
+
+      lnfug_val =   [2.0759140949373416, -2.2851989270402058]
+      dlnphidp_val = [-0.99059224575177762, -0.99388122357848807]
+      dlnphidt_val = [ 3.0263769083149254E-002,  7.6204871541712640E-002]
+
+      eos = binary_PR76()
+      z = [0.3, 0.7]
+
+      p = 1
+      t = 150
+
+      root_type = 1
+
+      call fugacity_tp(eos, &
+         z, T, P, V, root_type, lnfug, dlnPhidP, dlnphidT, dlnPhidn&
+      )
+
+      call check(&
+         error, maxval(rel_error(lnfug_val, lnfug)) < 1e-4 &
+      )
+      call check(&
+         error, maxval(rel_error(dlnphidp_val, dlnphidp)) < 1e-4 &
+      )
+      call check(&
+         error, maxval(rel_error(dlnphidt_val, dlnphidt)) < 1e-4 &
+      )
+   end subroutine
+
+   elemental function rel_error(x, y)
+      use yaeos, only: pr
+      real(pr), intent(in) :: x, y
+      real(pr) :: rel_error
+
+      rel_error = abs(x - y)/abs(x)
+
+   end function
+end module