MC alpha func

src/models/residual_helmholtz/ar_models.f90

@@ -291,23 +291,18 @@ subroutine fugacity_pt(eos, &
    end subroutine fugacity_pt
 
    subroutine fugacity_vt(eos, &
-      n, V, T, P, lnPhi, dlnPhidP, dlnPhidT, dlnPhidn, dPdV, dPdT, dPdn &
+      n, V, T, P, lnPhi, &
+      dlnPhidP, dlnPhidT, dlnPhidn, &
+      dPdV, dPdT, dPdn &
       )
       !! Calculate fugacity coefficent given volume and temperature.
-      !!
-      !!@note
-      !!While the natural output variable is \(ln \phi_i P\). The calculated
-      !!derivatives will be the derivatives of the fugacity coefficient
-      !!\(ln \phi_i\)
-      !!@endnote
-      !!
       class(ArModel) :: eos !! Model
       real(pr), intent(in) :: n(:) !! Mixture mole numbers
       real(pr), intent(in) :: V !! Volume [L]
       real(pr), intent(in) :: T !! Temperature [K]
 
       real(pr), optional, intent(out) :: P !! Pressure [bar]
-      real(pr), optional, intent(out) :: lnPhi(size(n)) !! \(\ln(\phi_i*P)\) vector
+      real(pr), optional, intent(out) :: lnPhi(size(n)) !! \(\ln(\phi_i)\) vector
       real(pr), optional, intent(out) :: dlnPhidT(size(n)) !! \(ln(phi_i)\) Temp derivative
       real(pr), optional, intent(out) :: dlnPhidP(size(n)) !! \(ln(phi_i)\) Presssure derivative
       real(pr), optional, intent(out) :: dlnPhidn(size(n), size(n)) !! \(ln(phi_i)\) compositional derivative

src/models/residual_helmholtz/cubic/alphas/alphas.f90

@@ -72,6 +72,8 @@ subroutine alpha_mc(self, Tr, a, dadt, dadt2)
 
       real(pr) :: sqrt_Tr(size(Tr))
 
+      real(pr) :: u(size(Tr)), dudt(size(Tr)), dudt2(size(Tr))
+
       sqrt_Tr = 1 - sqrt(Tr)
 
       ! The associate statement allows to abreviate the expresions
@@ -81,12 +83,13 @@ subroutine alpha_mc(self, Tr, a, dadt, dadt2)
             dadT = c1*(c1*(sqrt(Tr) - 1) - 1)/sqrt(Tr)
             dadT2 = (1.0_pr/2.0_pr)*c1*(c1 + 1)/Tr**(1.5_pr)
          elsewhere
-            a = (1 + c1 * (sqrt_Tr) + c2 * (sqrt_Tr) + c3 * (sqrt_Tr))**2
-            dadt = (c1 + c2 + c3) * (c1*(sqrt(Tr) - 1) &
-               + c2*(sqrt(Tr) - 1) + c3*(sqrt(Tr) - 1) - 1)/sqrt(Tr)
-            dadt2 = (1.0_pr/2.0_pr) * (&
-               c1**2 + 2*c1*c2 + 2*c1*c3 &
-               + c1 + c2**2 + 2*c2*c3 +  c2 + c3**2 + c3)/Tr**(3.0_pr/2.0_pr)
+            u = -c1*(sqrt(Tr) - 1) + c2*(sqrt(Tr) - 1)**2 - c3*(sqrt(Tr) - 1)**3 + 1
+            dudt = (1.0_pr/2.0_pr)*(-c1 + 2*c2*(sqrt(Tr) - 1) - 3*c3*(sqrt(Tr) - 1)**2) /sqrt(Tr)
+            dudt2 = (1.0_pr/4.0_pr)*(-3*Tr*c3 + c1 + 2*c2 + 3*c3)/Tr**(3.0_pr/2.0_pr)
+
+            a = u**2
+            dadt = 2*u*dudt
+            dadt2 = 2*(dudt**2 + u*dudt2)
          end where
       end associate
    end subroutine alpha_mc

test/test_psrk.f90

@@ -0,0 +1,65 @@
+program main
+   !! Running a PSRK cubic equation of state example/test
+   use yaeos
+
+   implicit none
+
+   ! ===========================================================================
+   ! Definition of variables that will be used
+   ! ---------------------------------------------------------------------------
+   type(CubicEoS) :: eos
+   type(Groups) :: molecules(2)
+   type(EquilibriumState) :: sat
+   type(PTEnvel2) :: env
+   real(pr) :: tc(2), pc(2), w(2), n(2)
+   real(pr) :: C(2, 3)
+
+   real(pr) :: v, lnphi(2)
+
+   real(pr) :: pressures(7) = [40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0]
+   real(pr) :: temperatures(7) = [450.0, 460.0, 470.0, 480.0, 490.0, 500.0, 510.0]
+
+   integer :: i
+
+   ! ===========================================================================
+   ! Definition of required parameters 
+   ! ---------------------------------------------------------------------------
+
+   ! Critical constants
+   tc = [304.21_pr, 553.8_pr]
+   pc = [7.383e6_pr, 4.080358e6_pr] / 1e5
+   w = [0.223621_pr, 0.213_pr]
+
+   ! Molecules groups
+   molecules(1)%groups_ids = [117]
+   molecules(1)%number_of_groups = [1]
+
+   molecules(2)%groups_ids = [2]
+   molecules(2)%number_of_groups = [6]
+
+   ! Mathias-Copeman parameters
+   C(1, :) = [0.8255_pr, 0.16755_pr, -1.7039_pr]
+   C(2, :) = [0.84082_pr, -0.39847_pr, 0.94148_pr]
+
+
+   ! Setting up the equation of state
+   eos = PSRK(&
+      tc=tc, &
+      pc=pc, &
+      w=w, &
+      molecules=molecules, &
+      c1=C(:, 1), c2=C(:, 2), c3=C(:, 3) &
+   )
+
+   ! Molar numbers of each component
+   n = [60, 40]
+
+   do i=1,7
+      call eos%lnphi_pt(n=n, P=pressures(i), T=temperatures(i), lnphi=lnphi, root_type="stable")
+      call eos%volume(n=n, P=pressures(i), T=temperatures(i), V=V, root_type="stable")
+      sat = saturation_pressure(eos, n, temperatures(i), kind="bubble", p0=130._pr)
+
+      print *, pressures(i), temperatures(i), exp(lnphi)
+      print *, pressures(i), temperatures(i), V
+   end do
+end program main