module yaeos__models_ar_cubic_cubic_mixing !! Cubic Mixing Rules for Cubic EoS. use yaeos__constants, only: pr, solving_volume use yaeos__substance, only: substances use yaeos__models_ar_genericcubic, only: CubicMixRule use yaeos__models_ar_cubic_mixing_base, only: bmix_qmr implicit none private public :: CMR public :: CMRTD public :: kijk_exp_tdep type, extends(CubicMixRule) :: CMR !! Cubic Mixing Rule (CMR) derived type. Classic Van der Waals mixing !! rules. !! !! QMR depends on binary interaction parameters, on a Cubic EoS !! the mixture is obtained by the combination of an attractive and !! repulsive parameter matrices. !! !! By default the attractive parameter matrix is calculated with: !! \[a_{ijk} = \sqrt[3]{a_i a_j a_k}(1 - k_{ijk})\] !! generating the \(a_{ijk}\) matrix, but this procedure can be overriden !! replacing the `aijk` pointer procedure. real(pr), allocatable :: k(:, :, :) !! Attractive Binary Interatction parameter matrix real(pr), allocatable :: l(:, :, :) !! Repulsive Binary Interatction parameter matrix contains procedure :: aijk => kijk_constant !! Default attractive parameter combining rule procedure :: Dmix !! Attractive parameter mixing rule procedure :: Bmix !! Repulsive parameter mixing rule procedure :: D1mix => RKPR_D1mix end type CMR type, extends(CMR) :: CMRTD real(pr), allocatable :: k0(:, :, :) real(pr), allocatable :: Tref(:, :, :) contains procedure :: aijk => kijk_exp_tdep end type CMRTD abstract interface subroutine get_aijk(& self, T, & ai, daidt, daidt2, & a, dadt, dadt2 & ) !! Combining rule for the attractive parameter. !! !! From previously calculated attractive parameters calculate the !! \(a_{ijk}\) matrix and it's corresponding derivatives. import pr, CMR class(CMR), intent(in) :: self real(pr), intent(in) :: T real(pr), intent(in) :: ai(:), daidt(:), daidt2(:) real(pr), intent(out):: a(:, :), dadt(:, :), dadt2(:, :) end subroutine get_aijk end interface contains subroutine Dmix(self, n, V, T, & ai, daidt, daidt2, & D, & dDdV, dDdT, dDdV2, dDdT2, dDi, dDdTV, dDidV, dDidT, dDij & ) !! Attractive parameter mixing rule with quadratic mix. !! !! Takes the all the pure components attractive parameters and their !! derivatives with respect to temperature and mix them with the !! Van der Waals quadratic mixing rule: !! !! \[ !! D = \sum_i \sum_j \sum_k n_i n_j n_k a_{ijk} = n^3 a_{mix} / n !! \] !! !! Inside the routine the \(a_{ijk}\) matrix is calculated using the !! procedure contained in the `CMR` object, this procedures defaults !! to the common combining rule: !! \(a_{ijk} = \sqrt[3]{a_i a_j a_k} (1 - k_{ijk}) \) !! !! The procedure can be overloaded by a common one that respects the !! interface [[get_aijk(interface)]] !! !! ```fortran !! type(CMR) :: my_mixing_rule !! my_mixing_rule%aij => new_aij_procedure !! ``` use yaeos__models_ar_cubic_mixing_base, only: CMR_Dmix class(CMR), intent(in) :: self !! Mixing rule object. 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) :: ai(:) !! Pure components attractive parameters \(a_i\) real(pr), intent(in) :: daidt(:) !! \(\frac{da_i}{dT}\) real(pr), intent(in) :: daidt2(:) !! \(\frac{d^2a_i}{dT^2}\) 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}}\) integer :: i, j, nc real(pr) :: a(size(ai), size(ai), size(ai)) real(pr) :: dadt(size(ai), size(ai), size(ai)) real(pr) :: dadt2(size(ai), size(ai), size(ai)) nc = size(ai) call self%aijk(T, ai, daidt, daidt2, a, dadt, dadt2) call CMR_Dmix(& n, V, T, & a, dadt, dadt2, & D=D, & dDdV=dDdV, & dDdT=dDdT, & dDdT2=dDdT2, & dDdV2=dDdV2, & dDdTV=dDdTV, & dDi=dDi, & dDidV=dDidV, & dDidT=dDidT, & dDij=dDij & ) end subroutine Dmix subroutine Bmix(self, n, bi, B, dBi, dBij) !! Mixture repulsive parameter. !! !! Calculate the mixture's repulsive parameter and it's derivatives !! with respect to composition: !! !! \[ !! n^2B = \sum_i \sum_j \sum_k n_i n_j n_k !! \frac{b_i + b_j + b_k}{3} (1 - l_{ijk}) !! \] !! use yaeos__models_ar_cubic_mixing_base, only: CMR_Bmix class(CMR), intent(in) :: self !! Mixing rule object. real(pr), intent(in) :: n(:) !! Moles vector. real(pr), intent(in) :: bi(:) !! Pure components repulsive parameters. real(pr), intent(out) :: B !! Mixture repulsive parameter. real(pr), intent(out) :: dBi(:) !! \(\frac{dB}{dn_i}\) real(pr), intent(out) :: dBij(:, :) !!\(\frac{d^2B}{dn_{ij}}\) real(pr) :: bijk(size(n), size(n), size(n)) integer :: i, j, l, nc nc = size(n) do i=1,nc do j=i,nc do l=i,nc bijk(i, j, l) = (bi(i) + bi(j) + bi(l))/3._pr * (1 - self%l(i, j, l)) bijk(i, l, j) = bijk(i, j, l) bijk(j, i, l) = bijk(i, j, l) bijk(j, l, i) = bijk(i, j, l) bijk(l, i, j) = bijk(i, j, l) bijk(l, j, i) = bijk(i, j, l) end do end do end do call CMR_Bmix(n, bijk, b, dbi, dbij) end subroutine Bmix subroutine RKPR_D1mix(self, n, d1i, D1, dD1i, dD1ij) use yaeos__models_ar_cubic_mixing_base, only: d1mix_rkpr !! RKPR \(\delta_1\) parameter mixing rule. !! !! The RKPR EoS doesn't have a constant \(\delta_1\) value for each !! component, so a proper mixing rule should be provided. A linear !! combination is used. !! !! \[ !! \Delta_1 = \sum_i^N n_i \delta_{1i} !! \] !! class(CMR), 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 RKPR_D1mix subroutine kijk_constant(& self, T, ai, daidt, daidt2, & a, dadt, dadt2 & ) !! Combining rule that uses constant \(k_{ijk}\) values. !! !! \[ !! a_{ijk} = \sqrt[3]{a_i a_j a_k} (1 - k_{ijk}) !! ] use yaeos__models_ar_cubic_mixing_base, only: CMR_aijk class(CMR), intent(in) :: self real(pr), intent(in) :: T !! Temperature [K] real(pr), intent(in) :: ai(:) !! Pure components attractive parameters (\a_i\) real(pr), intent(in) :: daidt(:) !! \(\frac{da_i}{dT}\) real(pr), intent(in) :: daidt2(:) !! \(\frac{d^2a_i}{dT^2}\) real(pr), intent(out) :: a(:, :, :) !! \(a_{ijk}\) Matrix real(pr), intent(out) :: dadt(:, :, :) !! \(\frac{da_{ijk}{dT}\) real(pr), intent(out) :: dadt2(:, :, :)!! \(\frac{d^2a_{ijk}{dT^2}\) integer :: i, j, l real(pr) :: k(size(ai), size(ai), size(ai)) real(pr) :: zeros(size(ai), size(ai), size(ai)) k = self%k zeros = 0 call CMR_aijk(ai, daidt, daidt2, k, zeros, zeros, a, dadt, dadt2) end subroutine kijk_constant subroutine kijk_exp_tdep(& self, T, ai, daidt, daidt2, & a, dadt, dadt2 & ) !! Combining rule that uses \(k_{ijk}\) as a function of temperature. !! !! \[ !! k_{ijk} = k_{ijk} + k_{ijk}^0 \exp \left(-T/T^{*}\) !! \] !! !! \[ !! a_{ijk} = \sqrt[3]{a_i a_j a_k} (1 - k_{ijk}) !! \] use yaeos__models_ar_cubic_mixing_base, only: CMR_aijk class(CMRTD), intent(in) :: self real(pr), intent(in) :: T !! Temperature [K] real(pr), intent(in) :: ai(:) !! Pure components attractive parameters (\a_i\) real(pr), intent(in) :: daidt(:) !! \(\frac{da_i}{dT}\) real(pr), intent(in) :: daidt2(:) !! \(\frac{d^2a_i}{dT^2}\) real(pr), intent(out) :: a(:, :, :) !! \(a_{ijk}\) Matrix real(pr), intent(out) :: dadt(:, :, :) !! \(\frac{da_{ijk}{dT}\) real(pr), intent(out) :: dadt2(:, :, :)!! \(\frac{d^2a_{ijk}{dT^2}\) integer :: i, j, l, nc real(pr) :: c(size(ai), size(ai), size(ai)) real(pr) :: k(size(ai), size(ai), size(ai)) real(pr) :: dkdt(size(ai), size(ai), size(ai)) real(pr) :: dkdt2(size(ai), size(ai), size(ai)) nc = size(ai) k = 0 dkdt = 0 dkdt2 = 0 where(self%k0 /= 0 .or. self%k /= 0) c = self%k0 * exp(-T/self%Tref) k = self%k + c dkdt = -c / self%Tref dkdt2 = c / self%Tref**2 end where call CMR_aijk(ai, daidt, daidt2, k, dkdt, dkdt2, a, dadt, dadt2) end subroutine kijk_exp_tdep end module yaeos__models_ar_cubic_cubic_mixing