huron_vidal.f90 Source File


Source Code

module yaeos__models_cubic_mixing_rules_huron_vidal
   !! # Huron-Vidal (like) mixing rules module
   !! This module contains the mixing rules that are based/similar to the
   !! mixing rules defined by Huron-Vidal
   !!
   !! # Description
   !! Huron-Vidal presented a way to link a \(G^E\) model with a Cubic EoS
   !! mixing rule. This makes it possible to make good predictions on
   !! polar compounds containing mixtures.
   !!
   !! # Examples
   !!
   !! ```fortran
   !!  A basic code example
   !! ```
   !!
   !! # References
   !!
   use yaeos__constants, only: pr, R, solving_volume
   use yaeos__models_ar_genericcubic, only: CubicMixRule
   use yaeos__models_ar_cubic_quadratic_mixing, only: QMR
   use yaeos__models_ar_cubic_mixing_base, only: bmix_qmr
   use yaeos__models_ge, only: GeModel
   use yaeos__models_ge_nrtlhv, only: NRTLHV
   implicit none

   private

   public :: HV
   public :: MHV
   public :: DmixMHV
   public :: HV_NRTL, init_hvnrtl

   type, extends(CubicMixRule) :: HV
      class(GeModel), allocatable :: ge
      real(pr), allocatable :: del1(:)
      real(pr), allocatable :: bi(:)
   contains
      procedure :: Bmix => BmixHV
      procedure :: D1Mix => D1MixHV
      procedure :: Dmix => DmixHV
   end type HV

   type, extends(CubicMixRule) :: MHV
      !! # Michelsen's modified Huron-Vidal mixing rule
      !! Mixing rule at zero-pressure which allows for the inclusion of an
      !! excess-gibbs model.
      !!
      !! # Description
      !! This mixing rule is based on the aproximate zero-pressure limit
      !!  of a cubic equation of state. At the aproximate zero-pressure limit the
      !! attractive parameter can be expressed as:
      !!
      !! \[
      !! \frac{D}{RTB}(n, T) = \sum_i n_i \frac{a_i(T)}{b_i} + \frac{1}{q}
      !!  \left(\frac{G^E(n, T)}{RT} + \sum_i n_i \ln \frac{B}{nb_i} \right)
      !! \]
      !! Where \(q\) is a weak function of temperature. In the case of `MHV`
      !! and simplicity it is considered that depends on the model used.
      !!
      !! # Examples
      !! To use the modified Huron-Vidal mixing rule it is necessary to define
      !! a `CubicEoS` and replace its original mixing rule with the one generated
      !! by the user.
      !! ```fortran
      !! type(MHV) :: mixrule
      !! type(NRTL) :: ge_model
      !! type(CubicEoS) :: model
      !!
      !! ! Define the Ge model to be used and the CubicEoS
      !! ge_model = NRTL(a, b, c)
      !! model = SoaveRedlichKwong(tc, pc, w)
      !!
      !! ! Use the initialization function to setup
      !! mixrule = MHV(ge=ge_model, q=-0.593_pr, bi=model%b)
      !!
      !! ! Replace the original mixrule on the previously defined model
      !! model%mixrule = mixrule
      !!
      !! ! Ready to do calculations
      !! call pressure(model, n, v, T)
      !! ```
      !!
      !! # References
      !!
      real(pr), allocatable :: l(:, :)
      real(pr), private, allocatable :: bi(:)
      real(pr), private, allocatable :: B, dBi(:), dBij(:, :)
      class(GeModel), allocatable :: ge
      real(pr) :: q
   contains
      procedure :: Bmix => BmixMHV
      procedure :: D1Mix => D1MixMHV
      procedure :: Dmix => DmixMHV
   end type MHV

   type, extends(CubicMixRule) :: HV_NRTL
      !! # HV_NRTL
      !! Huron-Vidal mixing rule including the NRTL model modified by Huron
      !! and Vidal.
      !!
      !! # Description
      !! This is the Huron-Vidal mixing rule that includes the NRTL model
      !! modified by Huron and Vidal. It is a mixing rule that allows to
      !! use the NRTL model as an excess Gibbs energy model and can. be
      !! simplified to the classic Quatratic mixing rules when the parameters
      !! are set to:
      !!
      !! \[
      !!  \alpha_{ji} = 0
      !! \]
      !!
      !! \[
      !!  g_{ii} = -\frac{a_i}{b_i} \lambda
      !! \]
      !!
      !! \[
      !!  g_{ji} = -2\frac{\sqrt{b_i b_j}}{b_i + b_j} \sqrt{g_{ii}g_{jj}}
      !!           \left(1 - k_{ij})\right)
      !! \]
      !!
      !! # Examples
      !!
      !! # References
      type(NRTLHV) :: ge
      real(pr), allocatable :: del1(:)
      real(pr), allocatable :: bi(:)
      logical, allocatable :: use_kij(:, :)
      real(pr), allocatable :: kij(:, :)
   contains
      procedure :: Bmix => BmixHVNRTL
      procedure :: D1Mix => D1MixHVNRTL
      procedure :: Dmix => DmixHVNRTL
   end type HV_NRTL

   interface MHV
      module procedure :: init_mhv
   end interface MHV

   interface HV_NRTL
      module procedure :: init_hvnrtl
   end interface HV_NRTL

contains

   ! ===========================================================================
   ! Huron-Vidal Mixing rule
   ! ---------------------------------------------------------------------------
   subroutine BmixHV(self, n, bi, B, dBi, dBij)
      !! # Repulsive parameter \(B\) mixing rule
      !! Quadratinc mixing rule for the repulsive parameter.
      !!
      !! # Description
      !! \[B = \sum_i n_i b_i\]
      use yaeos__models_ar_cubic_mixing_base, only: bmix_linear
      class(HV), intent(in) :: self
      real(pr), intent(in) :: n(:)
      real(pr), intent(in) :: bi(:)
      real(pr), intent(out) :: B, dBi(:), dBij(:, :)
      call bmix_linear(n, bi, b, dbi, dbij)
   end subroutine BmixHV

   subroutine D1MixHV(self, n, d1i, D1, dD1i, dD1ij)
      use yaeos__models_ar_cubic_mixing_base, only: d1mix_rkpr
      class(HV), intent(in) :: self
      real(pr), intent(in) :: n(:)
      real(pr), intent(in) :: d1i(:)
      real(pr), intent(out) :: D1
      real(pr), intent(out) :: dD1i(:)
      real(pr), intent(out) :: dD1ij(:, :)
      call d1mix_rkpr(n, d1i, D1, dD1i, dD1ij)
   end subroutine D1MixHV

   subroutine DmixHV(self, n, T, &
      ai, daidt, daidt2, &
      D, dDdT, dDdT2, dDi, dDidT, dDij &
      )
      use yaeos__models_ar_cubic_mixing_base, only: lamdba_hv, mix => DMixHV
      class(HV), intent(in) :: self
      real(pr), intent(in) :: T, n(:)
      real(pr), intent(in) :: ai(:), daidt(:), daidt2(:)
      real(pr), intent(out) :: D, dDdT, dDdT2, dDi(:), dDidT(:), dDij(:, :)

      real(pr) :: b, bi(size(n)), dbi(size(n)), dbij(size(n), size(n))
      real(pr) :: del1(size(n)), del2(size(n))
      real(pr) :: d1, d1i(size(n)), dd1i(size(n)), dd1ij(size(n), size(n))
      real(pr) :: Ge, GeT, GeT2, Gen(size(n)), GeTn(size(n)), Gen2(size(n), size(n))

      real(pr) :: totn !! Total number of moles

      integer :: i, j, nc
      real(pr) :: L, dL(size(n)), dL2(size(n), size(n))

      nc = size(n)
      totn = sum(n)

      del1 = self%del1
      bi = self%bi

      call self%Bmix(n, bi, B, dBi, dBij)
      call self%D1Mix(n, del1, D1, dD1i, dD1ij)
      call lamdba_hv(nc, D1, dD1i, dD1ij, L, dL, dL2)

      call self%ge%excess_gibbs( &
         n, T, Ge=Ge, GeT=GeT, GeT2=GeT2, Gen=Gen, GeTn=GeTn, Gen2=Gen2 &
         )
      call mix(n, T,&
         bi, B, dBi, dBij, &
         D1, dD1i, dD1ij, &
         ai, daidt, daidt2, &
         Ge, GeT, GeT2, Gen, GeTn, Gen2,&
         D, dDdT, dDdT2, dDi, dDidT, dDij)
   end subroutine DmixHV

   ! ===========================================================================
   ! Huron-Vidal Mixing rule with Huron-Vidal NRTL
   ! ---------------------------------------------------------------------------
   type(HV_NRTL) function init_hvnrtl(b, del1, alpha, gji, use_kij, kij) result(mixrule)
      !! # Huron-Vidal NRTL mixing rule
      !! This is the Huron-Vidal mixing rule that includes the NRTL model
      !! modified by Huron and Vidal.
      !!
      !! # Description
      !! This is the Huron-Vidal mixing rule that includes the NRTL model
      !! modified by Huron and Vidal. It is a mixing rule that allows to
      !! use the NRTL model as an excess Gibbs energy model and can. be
      !! simplified to the classic Quatratic mixing rules when the parameters
      !! are set to:
      !! \[
      !!  \alpha_{ji} = 0
      !! \]
      !!
      !! \[
      !!  g_{ii} = -\frac{a_i}{b_i} \lambda
      !! \]
      !!
      !! \[
      !!  g_{ji} = -2\frac{\sqrt{b_i b_j}}{b_i + b_j} \sqrt{g_{ii}g_{jj}}
      !!           \left(1 - k_{ij})\right)
      !! \]
      !!
      !! # Examples
      !!
      use yaeos__models_ge_nrtlhv, only: NRTLHV
      real(pr), intent(in) :: b(:)
      real(pr), intent(in) :: del1(:)
      real(pr), intent(in) :: alpha(:, :)
      real(pr), intent(in) :: gji(:, :)
      logical, intent(in) :: use_kij(:, :)
      real(pr), intent(in) :: kij(:, :)

      integer :: i, nc

      nc = size(b)

      mixrule%ge = NRTLHV(b=b, alpha=alpha, gij=gji)
      mixrule%bi = b
      mixrule%del1 = del1
      mixrule%use_kij = use_kij
      mixrule%kij = kij
   end function init_hvnrtl

   subroutine BmixHVNRTL(self, n, bi, B, dBi, dBij)
      !! # Repulsive parameter \(B\) mixing rule
      !! Quadratinc mixing rule for the repulsive parameter.
      !!
      !! # Description
      !! \[B = \sum_i n_i b_i\]
      use yaeos__models_ar_cubic_mixing_base, only: bmix_linear
      class(HV_NRTL), intent(in) :: self
      real(pr), intent(in) :: n(:)
      real(pr), intent(in) :: bi(:)
      real(pr), intent(out) :: B, dBi(:), dBij(:, :)
      call bmix_linear(n, bi, b, dbi, dbij)
   end subroutine BmixHVNRTL

   subroutine D1MixHVNRTL(self, n, d1i, D1, dD1i, dD1ij)
      use yaeos__models_ar_cubic_mixing_base, only: d1mix_rkpr
      class(HV_NRTL), intent(in) :: self
      real(pr), intent(in) :: n(:)
      real(pr), intent(in) :: d1i(:)
      real(pr), intent(out) :: D1
      real(pr), intent(out) :: dD1i(:)
      real(pr), intent(out) :: dD1ij(:, :)
      call d1mix_rkpr(n, d1i, D1, dD1i, dD1ij)
   end subroutine D1MixHVNRTL

   subroutine DmixHVNRTL(self, n, T, &
      ai, daidt, daidt2, &
      D, dDdT, dDdT2, dDi, dDidT, dDij &
      )
      use yaeos__models_ar_cubic_mixing_base, only: lamdba_hv, DmixHV
      use yaeos__models_ge_nrtlhv, only: NRTLHV
      class(HV_NRTL), intent(in) :: self
      real(pr), intent(in) :: T, n(:)
      real(pr), intent(in) :: ai(:), daidt(:), daidt2(:)
      real(pr), intent(out) :: D, dDdT, dDdT2, dDi(:), dDidT(:), dDij(:, :)

      real(pr) :: Ge, GeT, GeT2
      real(pr) :: Gen(size(n)), GeTn(size(n)), Gen2(size(n), size(n))

      real(pr) :: B, dBi(size(n)), dBij(size(n), size(n))
      real(pr) :: D1, dD1i(size(n)), dD1ij(size(n), size(n))
      real(pr) :: L

      type(NRTLHV) :: ge_model

      real(pr) :: gii(size(n)), gji(size(n), size(n))
      real(pr) :: bi(size(n))

      integer :: i, j, nc

      ge_model = self%ge

      nc = size(n)
      bi = self%bi
      call self%Bmix(n, bi, B, dBi, dBij)
      call self%D1Mix(n, self%del1, D1, dD1i, dD1ij)
      call lamdba_hv(nc, D1, L=L)

      gii = - ai/bi * L

      do i=1,nc
         do j=1,nc
            if (self%use_kij(i, j)) then
               ge_model%alpha(i, j) = 0
               ge_model%gij(i, j) = -2 * sqrt(bi(i) * bi(j)) / (bi(i) + bi(j)) * &
                  sqrt(gii(i) * gii(j)) * (1 - self%kij(i, j)) - gii(j)
            end if
         end do
      end do

      call ge_model%excess_gibbs( &
         n, T, Ge=Ge, GeT=GeT, GeT2=GeT2, Gen=Gen, GeTn=GeTn, Gen2=Gen2 &
         )

      call DMixHV(n, T,&
         bi, B, dBi, dBij, &
         D1, dD1i, dD1ij, &
         ai, daidt, daidt2, &
         Ge, GeT, GeT2, Gen, GeTn, Gen2,&
         D, dDdT, dDdT2, dDi, dDidT, dDij)
   end subroutine DmixHVNRTL

   ! ===========================================================================
   ! Michelsen's Modified Huron-Vidal 1
   ! ---------------------------------------------------------------------------
   type(MHV) function init_mhv(ge, b, q, lij) result(mixrule)
      class(GeModel), intent(in) :: Ge
      real(pr), intent(in) :: b(:)
      real(pr), intent(in) :: q
      real(pr), optional, intent(in) :: lij(:, :)

      integer :: i, nc

      nc = size(b)

      mixrule%q = q
      mixrule%bi = b
      mixrule%Ge = ge
      if (present(lij)) then
         mixrule%l = lij
      else
         mixrule%l = reshape([(0, i=1, nc**2)], [nc, nc])
      end if
   end function init_mhv

   subroutine BmixMHV(self, n, bi, B, dBi, dBij)
      !! # Repulsive parameter \(B\) mixing rule
      !! Quadratinc mixing rule for the repulsive parameter, using
      !! \( b_{ij} = \frac{b_i + b_j}{2} (1 - l_{ij}) \) as a combining rule.
      !!
      !! # Description
      !! Michelsen's modified Huron-Vidal mixing rule assumes a linear mix of
      !! the repulsive parameter.
      !!
      !! \[B = \sum_i n_i b_i\]
      !!
      !! In this implementation the most known crossed combining rule is used:
      !! \[nB = \sum_i \sum_j \frac{b_i + b_j}{2} (1 - l_{ij})\]
      !! to provide versatility to the used model.
      !!
      !! @warning
      !! This mixing rule is intended to use only with a linear combining
      !! rule, using \(l_{ij}\) could negatively affect the thermodynamic
      !! consistency of the model.
      !! @endwarning
      !!
      !! # Examples
      !!
      !! ```fortran
      !!  A basic code example
      !! ```
      !!
      !! # References
      !!
      use yaeos__models_ar_cubic_mixing_base, only: bmix_linear
      class(MHV), intent(in) :: self
      real(pr), intent(in) :: n(:)
      real(pr), intent(in) :: bi(:)
      real(pr), intent(out) :: B, dBi(:), dBij(:, :)
      call bmix_qmr(n, bi, self%l, b, dbi, dbij)
      ! call bmix_linear(n, bi, b, dbi, dbij)
   end subroutine BmixMHV

   subroutine DmixMHV(self, n, T, &
      ai, daidt, daidt2, &
      D, dDdT, dDdT2, dDi, dDidT, dDij &
      )
      !! # Michelsen Modified Huron-Vidal mixing rule.
      !! Mixing rule at infinite pressure as defined in the book of Michelsen and
      !! Møllerup.
      !!
      !! # Description
      !! At the infinite pressure limit of a cubic equation of state it is possible to
      !! relate teh mixing rule for the attractive term with a excess Gibbs energy
      !! model like NRTL with the expression:
      !!
      !! \[
      !! \frac{D}{RTB}(n, T) = \sum_i n_i \frac{a_i(T)}{b_i} + \frac{1}{q}
      !!  \left(\frac{G^E(n, T)}{RT} + \sum_i n_i \ln \frac{B}{nb_i} \right)
      !! \]
      !!
      !! # Examples
      !!
      !! ```fortran
      !!  type(CubicEoS)
      !! ```
      !!
      !! # References
      !!
      class(MHV), intent(in) :: self
      real(pr), intent(in) :: T, n(:)
      real(pr), intent(in) :: ai(:), daidt(:), daidt2(:)
      real(pr), intent(out) :: D, dDdT, dDdT2, dDi(:), dDidT(:), dDij(:, :)
      real(pr) :: f, fdt, fdt2, fdi(size(n)), fdit(size(n)), fdij(size(n), size(n))

      real(pr) :: b, bi(size(n)), dbi(size(n)), dbij(size(n), size(n))
      real(pr) :: Ge, GeT, GeT2, Gen(size(n)), GeTn(size(n)), Gen2(size(n), size(n))

      real(pr) :: totn !! Total number of moles
      real(pr) :: dot_n_logB_nbi
      real(pr) :: logB_nbi(size(n)) !! \(\ln \frac{B}{n b_i}\)
      real(pr) :: dlogBi_nbi(size(n))
      real(pr) :: d2logBi_nbi(size(n), size(n))

      integer :: i, j, l, nc
      real(pr) :: q

      nc = size(n)
      totn = sum(n)

      q = self%q
      bi = self%bi

      if (.not. solving_volume) then
         call self%ge%excess_gibbs( &
            n, T, Ge=Ge, GeT=GeT, GeT2=GeT2, Gen=Gen, GeTn=GeTn, Gen2=Gen2 &
            )
      else
         call self%ge%excess_gibbs( &
            n, T, Ge=Ge &
            )
      end if
      call self%Bmix(n, bi, B, dBi, dBij)
      logb_nbi = log(B/(totn*bi))
      dot_n_logB_nbi = dot_product(n, logB_nbi)

      do i = 1, nc
         dlogBi_nbi(i) = logB_nbi(i) + sum(n*dBi(i))/B - 1
      end do

      if (.not. solving_volume) then
         do i = 1, nc
            do j = 1, nc
               !TODO: Need to figure out this derivative
               d2logBi_nbi(i, j) = dlogBi_nbi(j) &
                  + (sum(n*dBij(i, j)) + dBi(i))/B &
                  - totn*dBi(i)*dBi(j)/B**2
            end do
         end do

         autodiff: block
            !! Autodiff injection until we can decipher this derivative
            use hyperdual_mod
            type(hyperdual) :: hB
            type(hyperdual) :: hdot_ln_B_nbi
            type(hyperdual) :: hn(nc)

            integer :: ii, jj
            hn = n

            do i = 1, nc
               do j = i, nc
                  hn = n
                  hn(i)%f1 = 1
                  hn(j)%f2 = 1

                  hB = 0._pr
                  do ii=1,nc
                     do jj=1,nc
                        hB = hB &
                           + (hn(ii)*hn(jj)) &
                           * 0.5_pr * (bi(ii) + bi(jj)) * (1._pr - self%l(ii, jj))
                     end do
                  end do
                  hB = hB/sum(hn)

                  hdot_ln_B_nbi = sum(hn*log(hB/(sum(hn)*bi)))

                  d2logBi_nbi(i, j) = hdot_ln_B_nbi%f12
                  d2logBi_nbi(j, i) = hdot_ln_B_nbi%f12
               end do
            end do
         end block autodiff
      end if

      f = sum(n*ai/bi) + (Ge + R*T*dot_n_logB_nbi)/q
      fdt = sum(n*daidt/bi) + (GeT + R*dot_n_logB_nbi)/q
      fdt2 = sum(n*daidt2/bi) + (GeT2)/q

      fdi = ai/bi + (1._pr/q)*(GeN + R*T*(dlogBi_nbi))
      fdit = daidt/bi + (1._pr/q)*(GeTn + R*(dlogBi_nbi))

      do i = 1, nc
         do j = 1, nc
            fdij(i, j) = R*T*(d2logBi_nbi(i, j))
            fdij(i, j) = 1/q*(fdij(i, j) + GeN2(i, j))
            fdij(i, j) = &
               dBi(j)*fdi(i) + B*fdij(i, j) + fdi(j)*dBi(i) + f*dBij(i, j)
         end do
      end do

      dDi = B*fdi + f*dBi
      dDidT = B*fdiT + fdT*dBi

      D = f*B
      dDdT = fdT*B
      dDdT2 = fdT2*B
      dDij = fdij

   end subroutine DmixMHV

   subroutine D1MixMHV(self, n, d1i, D1, dD1i, dD1ij)
      use yaeos__models_ar_cubic_mixing_base, only: d1mix_rkpr
      class(MHV), intent(in) :: self
      real(pr), intent(in) :: n(:)
      real(pr), intent(in) :: d1i(:)
      real(pr), intent(out) :: D1
      real(pr), intent(out) :: dD1i(:)
      real(pr), intent(out) :: dD1ij(:, :)
      call d1mix_rkpr(n, d1i, D1, dD1i, dD1ij)
   end subroutine D1MixMHV
end module yaeos__models_cubic_mixing_rules_huron_vidal