Implementation of Cubic Mixing Rules

src/extra/fluids.f90

@@ -21,7 +21,7 @@ module yaeos__extra_fluids
 
    type :: FixtureFluid
       integer :: nc
-         !! Number of components
+      !! Number of components
       character(len=:), allocatable :: reference
       !! Reference of where the fluid was taken from
       real(pr), allocatable :: z0(:)
@@ -69,12 +69,30 @@ type(FixtureFluid) function oil_gao()
 
       mixrule = QMR(k=kij, l=lij)
 
-
       call ar_model%set_mixrule(mixrule)
 
       oil_gao%nc = nc
       oil_gao%reference = "doi: 10.2118/180736-pa"
       oil_gao%ar_model = ar_model
       oil_gao%z0 = [0.6202_pr, 0.3702_pr, 0.0051_pr, 0.0023_pr, 0.0014_pr, 0.0008_pr]
    end function oil_gao
+
+   type(FixtureFluid) function co2_h2o_isop() result(fluid)
+      use yaeos__models, only: CubicEoS, SoaveRedlichKwong
+      integer, parameter :: nc = 3
+      real(pr) :: Tc(nc), Pc(nc), w(nc), kij(nc, nc)
+      Tc = [304.21, 647.13, 508.3]
+      Pc = [73.3, 220.55, 47.64]
+      w = [0.2236, 0.3449, 0.6669]
+
+      kij(1, :) = [0._pr, 0.19_pr, 0.1215_pr]
+      kij(2, :) = [0.19_pr, 0._pr, -0.1727_pr]
+      kij(3, :) = [0.1215_pr, -0.1727_pr, 0._pr]
+
+      fluid%nc = nc
+      fluid%ar_model = SoaveRedlichKwong(Tc, Pc, w, kij=kij)
+      fluid%z0 = [0.997493267163008, 0.000567120831996, 0.00193961200499]
+      fluid%zi = [0.122878702720943, 0.376116603168507, 0.50100469411055]
+   end function co2_h2o_isop
+
 end module yaeos__extra_fluids

src/models/models.f90

@@ -19,6 +19,10 @@ module yaeos__models
    !!       defaults to classic vdW mixing rules.
    !!    - `MHV` (Modified Huron-Vidal) type: Michelsens first order modified
    !!       Huron-Vidal mixing rule.
+   !!    - `HV` (Huron-Vidal) type: Huron-Vidal mixing rule.
+   !!    - `sDDLC` (segmented Density Dependent Local Composition) type:
+   !!    mixing rule proposed by Cismondi and Mollerup, including volume as 
+   !!    a dependance of the mixture's attractive parameter. 
    !! - **GERG2008 Equation of State**:
    !!    - GERG2008 multifluid equation of state
    !! - **SAFT Equations of State**:
@@ -47,6 +51,7 @@ module yaeos__models
 
    ! Mixing Rules
    use yaeos__models_ar_cubic_quadratic_mixing
+   use yaeos__models_ar_cubic_cubic_mixing
    use yaeos__models_cubic_mixing_rules_huron_vidal
    use yaeos__models_ar_cubic_mixing_sddlc
 

src/models/residual_helmholtz/cubic/generic_cubic.f90

@@ -102,6 +102,7 @@ module yaeos__models_ar_genericcubic
       procedure :: set_delta1 => set_delta1
       procedure :: set_mixrule => set_mixrule
       procedure :: set_alpha => set_alpha
+      procedure :: calculate_attractive_parameters => calculate_attractive_parameters
    end type CubicEoS
 
    abstract interface
@@ -131,6 +132,7 @@ subroutine abs_Bmix(self, n, bi, B, dBi, dBij)
          real(pr), intent(in) :: bi(:)
          real(pr), intent(out) :: B, dBi(:), dBij(:, :)
       end subroutine abs_Bmix
+
       subroutine abs_D1mix(self, n, d1i, D1, dD1i, dD1ij)
          import pr, CubicMixRule
          class(CubicMixRule), intent(in) :: self
@@ -270,6 +272,23 @@ subroutine set_alpha(self, alpha)
       self%alpha = alpha
    end subroutine set_alpha
 
+   subroutine calculate_attractive_parameters(self, T, a, dadt, dadt2)
+      class(CubicEoS), intent(in out) :: self !! Model
+      real(pr), intent(in) :: T !! Temperature [K]
+      real(pr), intent(out) :: a(:) !! Attractive parameter
+      real(pr), intent(out) :: dadt(:) !! Derivative wrt temperature
+      real(pr), intent(out) :: dadt2(:) !! Second derivative wrt temperature
+
+      real(pr) :: Tr(size(a)) !! Reduced temperatures for each pure component
+
+      Tr = T/self%components%Tc
+      call self%alpha%alpha(Tr, a, dadt, dadt2)
+
+      a = self%ac * a
+      dadt = self%ac * dadt / self%components%Tc
+      dadt2 = self%ac * dadt2 / self%components%Tc**2
+   end subroutine calculate_attractive_parameters
+
    function v0(self, n, p, t)
       !! Cubic EoS volume initializer.
       !! For a Cubic Equation of State, the covolume calculated with the mixing

src/models/residual_helmholtz/cubic/mixing_rules/base.f90

@@ -172,7 +172,7 @@
       )
       !! # `DmixHV`
       !! Attractive parameter calculation for the Huron-Vidal mixing rule.
-      !! 
+      !!
       !! # Description
       !! This subroutine calculates the attractive parameter \(D\) and its
       !! derivatives for a mixture using the Huron-Vidal mixing rule.
@@ -181,8 +181,8 @@
       !! The expression of the attractive parameter is:
       !!
       !! \[
-      !!   D(n, T) = 
-      !!     B\left(\sum_i n_i\frac{a_i}{b_i} 
+      !!   D(n, T) =
+      !!     B\left(\sum_i n_i\frac{a_i}{b_i}
       !!     - \frac{G^E}{\Lambda}\right)
       !! \]
       !!
@@ -243,4 +243,236 @@
 
    end subroutine DmixHV
 
+   ! ===========================================================================
+   ! Cubic Mixing Rules
+   ! ---------------------------------------------------------------------------
+
+   pure subroutine CMR_aijk(&
+      ai, daidt, daidt2, k, dkdt, dkdt2, &
+      a, dadt, dadt2 &
+      )
+      !! # `CMR_aijk`
+      !! Calculate the \(a_{ijk}\) tensor of the CMR.
+      real(pr), intent(in) :: ai(:) !! \(\a_i\)
+      real(pr), intent(in) :: daidt(:) !! \( \frac{d\a_i}{dT} \)
+      real(pr), intent(in) :: daidt2(:)!! \( \frac{d^2\a_i}{dT^2} \)
+      real(pr), intent(in) :: k(:, :, :) !! k_{ijk} matrix
+      real(pr), intent(in) :: dkdt(:, :, :) !! \(  \frac{d k_{ijk}}{dT} \)
+      real(pr), intent(in) :: dkdt2(:, :, :) !! \(  \frac{d^2k_{ijk}}{dT^2} \)
+      real(pr), intent(out) :: a(:, :, :) !! \(a_{ijk}\) matrix
+      real(pr), intent(out) :: dadt(:, :, :) !! \(\frac{da_{ijk}}{dT}\) matrix
+      real(pr), intent(out) :: dadt2(:, :, :) !! \(\frac{d^2a_{ijk}{dT^2}\) matrix
+
+      integer :: i, j, l, nc
+
+      real(pr) :: aijk_3
+      real(pr) :: aaa, daaadt
+
+      real(pr), parameter :: third = 1._pr/3._pr
+
+      a = 0
+      dadt = 0
+      dadt2 = 0
+
+      nc = size(ai)
+
+      do i = 1,nc
+         aaa = ai(i) * ai(i) * ai(i)
+         aijk_3 = aaa ** third
+
+         a(i, i, i) =  aijk_3
+
+         daaadt = (&
+            daidt(i) * ai(i) * ai(i) &
+            + ai(i) * daidt(i) * ai(i) &
+            + ai(i) * ai(i) * daidt(i) &
+            ) * third
+
+         dadt(i, i, i) = daidt(i)
+         dadt2(i, i, i) = daidt2(i)
+
+         do j = i,nc
+            do l = i,nc
+               aaa = ai(i) * ai(j) * ai(l)
+               aijk_3 = aaa ** third
+
+               a(i, j, l) =  aijk_3 * (1 - k(i, j, l))
+
+               a(j, i, l) = a(i, j, l)
+               a(j, l, i) = a(i, j, l)
+
+               daaadt = (&
+                  daidt(i) * ai(j) * ai(l) &
+                  + ai(i) * daidt(j) * ai(l) &
+                  + ai(i) * ai(j) * daidt(l) &
+                  ) * third
+
+               dadt(i, j, l) = - dkdt(i, j, l) * aijk_3 + daaadt / (aaa) * a(i, j, l)
+               dadt(j, i, l) = dadt(i, j, l)
+               dadt(j, l, i) = dadt(i, j, l)
+
+               dadt2(i, j, l) = &
+                  - 2 * dkdt(i, j, l) * daaadt * aijk_3 / aaa &
+                  - daidt(i) * daaadt * a(i, j, l) / (aaa * ai(i)) &
+                  - daidt(j) * daaadt * a(i, j, l) / (aaa * ai(j)) &
+                  - daidt(l) * daaadt * a(i, j, l) / (aaa * ai(l)) &
+                  - dkdt2(i, j, l) * aijk_3 &
+                  + ((&
+                  2 * daidt(i) * daidt(j) * ai(l) &
+                  + 2 * daidt(i) * daidt(l) * ai(j) &
+                  + 2 * daidt(j) * daidt(l) * ai(i) &
+                  + daidt2(i) * ai(j) * ai(l) &
+                  + daidt2(j) * ai(i) * ai(l) &
+                  + daidt2(l) * ai(j) * ai(i) &
+                  )/3._pr) * (a(i, j, l) / aaa) &
+                  + (a(i, j, l)*daaadt**2/aaa**2)
+               dadt2(j, i, l) = dadt2(i, j, l)
+               dadt2(j, l, i) = dadt2(i, j, l)
+
+            end do
+         end do
+      end do
+   end subroutine CMR_aijk
+
+   subroutine CMR_Dmix(n, V, T, &
+      a, dadt, dadt2, &
+      D, &
+      dDdV, dDdT, dDdV2, dDdT2, dDi, dDdTV, dDidV, dDidT, dDij &
+      )
+      real(pr), intent(in) :: V !! Volume [L] (unused)
+      real(pr), intent(in) :: T !! Temperature [K]
+      real(pr), intent(in) :: n(:) !! Moles vector [mol]
+      real(pr), intent(in) :: a(:, :, :) !!
+      real(pr), intent(in) :: dadt(:, :, :) !!
+      real(pr), intent(in) :: dadt2(:, :, :) !!
+
+      real(pr), intent(out) :: D !! Mixture attractive parameter \(n^2a_{mix}\)
+      real(pr), intent(out) :: dDdV !! \(\frac{dD}{dT}\)
+      real(pr), intent(out) :: dDdT !! \(\frac{dD}{dV}\)
+      real(pr), intent(out) :: dDdT2 !! \(\frac{d^2D}{dT^2}\)
+      real(pr), intent(out) :: dDdV2 !! \(\frac{d^2D}{dV^2}\)
+      real(pr), intent(out) :: dDdTV !! \(\frac{d^2D}{dTV\)
+      real(pr), intent(out) :: dDi(:) !! \(\frac{dD}{dn_i}\)
+      real(pr), intent(out) :: dDidV(:) !! \(\frac{d^2D}{dVn_i}\)
+      real(pr), intent(out) :: dDidT(:) !! \(\frac{d^2D}{dTn_i}\)
+      real(pr), intent(out) :: dDij(:, :)!! \(\frac{d^2D}{dn_{ij}}\)
+
+      real(pr) :: aux(size(n)),auxij(size(n),size(n))
+      real(pr) :: auxT(size(n)),auxTij(size(n),size(n)),auxT2(size(n)),auxT2ij(size(n),size(n))
+      real(pr) :: aijk(size(n),size(n),size(n)),daijkdT(size(n),size(n),size(n)),daijkdT2(size(n),size(n),size(n))
+
+      real(pr) :: totn
+      real(pr) :: cum
+
+      integer :: i, j, k, nc
+
+      nc = size(n)
+
+      TOTN = sum(n)
+      D = 0.0_pr
+      dDdT = 0.0_pr
+      dDdT2 = 0.0_pr
+      aux = 0.0_pr
+      auxij = 0.0_pr
+      auxT = 0.0_pr
+      auxT2 = 0.0_pr
+      auxTij = 0.0_pr
+
+      aijk = a
+      daijkdT = dadt
+      daijkdT2 = dadt2
+
+      do i=1,nc
+         do j=1,nc
+            do k=1,nc
+               auxij(i,j)=auxij(i,j) + n(k)*aijk(i,j,k)
+
+               auxTij(i,j)=auxTij(i,j) + n(k)*daijkdT(i,j,k)
+               auxT2ij(i,j)=auxT2ij(i,j) + n(k)*daijkdT2(i,j,k)
+            end do
+
+            aux(i) = aux(i)+n(j)*auxij(i,j)
+
+            auxT(i) = auxT(i) + n(j)*auxTij(i,j)
+            auxT2(i) = auxT2(i) + n(j)*auxT2ij(i,j)
+
+         end do
+
+         D = D + n(i)*aux(i)
+
+         dDdT = dDdT + n(i)*auxT(i)
+         dDdT2 = dDdT2 + n(i)*auxT2(i)
+      end do
+
+      D = D/TOTN
+      dDdT = dDdT/TOTN
+      dDdT2 = dDdT2/TOTN
+
+      do i=1,nc
+         dDi(i) = (3*aux(i)-D)/TOTN
+
+         do j=1,i
+            dDij(i,j) = (6*auxij(i,j)-dDi(i)-dDi(j)) / totn
+            dDij(j,i) = dDij(i,j)
+         end do
+
+      end do
+
+      ! D = 0
+      ! dDdT = 0
+      ! dDdT2 = 0
+      ! do i=1,nc
+      !    do j=1,nc
+      !       do k=1,nc
+      !          D = D + n(i) * n(j) * n(k) * a(i, j, k)
+      !          dDdT = dDdT + n(i) * n(j) * n(k) * dadt(i, j, k)
+      !          dDdT2 = dDdT2 + n(i) * n(j) * n(k) * dadt2(i, j, k)
+      !       end do
+      !    end do
+      ! end do
+
+      ! D = D / totn
+      ! dDdT = dDdT / totn
+      ! dDdT2 = dDdT2 / totn
+   end subroutine CMR_Dmix
+
+   pure subroutine CMR_Bmix(n, bijk, Bmix, dBi, dBij)
+      real(pr), intent(in) :: n(:)
+      real(pr), intent(in) :: bijk(:, :, :)
+      real(pr), intent(out) :: Bmix
+      real(pr), intent(out) :: dBi(:), dBij(:,:)
+      real(pr) :: aux(size(n)), auxij(size(n), size(n))
+
+      real(pr) :: totn !! Total number of moles
+      real(pr) :: sqn !! totn^2
+
+      integer :: i, j, k, nc
+
+      nc = size(n)
+      totn = sum(n)
+      sqn = totn*totn
+      Bmix = 0
+      aux = 0
+      auxij = 0
+
+      do i=1,nc
+         do j=1,nc
+            do k=1,nc
+               auxij(i,j) = auxij(i,j) + n(k)*bijk(i,j,k)
+            end do
+            aux(i) = aux(i) + n(j)*auxij(i,j)
+         end do
+         Bmix = Bmix + n(i)*aux(i)
+      end do
+
+      Bmix=Bmix/sqn
+
+      do i=1,nc
+         dBi(i)=(3*aux(i)-2*totn*Bmix)/sqn
+         do j=1,i
+            dBij(i,j)=(6*auxij(i,j)-2*(Bmix+totn*dBi(i)+totn*dBi(j)))/sqn
+            dBij(j,i)=dBij(i,j)
+         end do
+      end do
+   end subroutine CMR_Bmix
 end module yaeos__models_ar_cubic_mixing_base