diff --git a/src/extra/fluids.f90 b/src/extra/fluids.f90 index ebdf430f2..971d36af2 100644 --- a/src/extra/fluids.f90 +++ b/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,7 +69,6 @@ type(FixtureFluid) function oil_gao() mixrule = QMR(k=kij, l=lij) - call ar_model%set_mixrule(mixrule) oil_gao%nc = nc @@ -77,4 +76,23 @@ type(FixtureFluid) function oil_gao() 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 diff --git a/src/models/models.f90 b/src/models/models.f90 index 72d7fcd13..9bec1ef28 100644 --- a/src/models/models.f90 +++ b/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 diff --git a/src/models/residual_helmholtz/cubic/generic_cubic.f90 b/src/models/residual_helmholtz/cubic/generic_cubic.f90 index 75f183916..c400dc7d4 100644 --- a/src/models/residual_helmholtz/cubic/generic_cubic.f90 +++ b/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 diff --git a/src/models/residual_helmholtz/cubic/mixing_rules/base.f90 b/src/models/residual_helmholtz/cubic/mixing_rules/base.f90 index 574dd7c6c..49723405c 100644 --- a/src/models/residual_helmholtz/cubic/mixing_rules/base.f90 +++ b/src/models/residual_helmholtz/cubic/mixing_rules/base.f90 @@ -172,7 +172,7 @@ subroutine DmixHV(n, T, & ) !! # `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 @@ subroutine DmixHV(n, T, & !! 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 @@ subroutine DmixHV(n, T, & 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 diff --git a/src/models/residual_helmholtz/cubic/mixing_rules/cubic_mixing.f90 b/src/models/residual_helmholtz/cubic/mixing_rules/cubic_mixing.f90 new file mode 100644 index 000000000..616a05e8d --- /dev/null +++ b/src/models/residual_helmholtz/cubic/mixing_rules/cubic_mixing.f90 @@ -0,0 +1,267 @@ +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)) + 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 diff --git a/test/test_cmr.f90 b/test/test_cmr.f90 new file mode 100644 index 000000000..4e5f0d3a2 --- /dev/null +++ b/test/test_cmr.f90 @@ -0,0 +1,443 @@ +module tests_cmr + use yaeos + use testing_aux, only: assert + integer, parameter :: nc = 3 +contains + subroutine test_aijk_numdiff + use yaeos__models_ar_cubic_mixing_base, only: CMR_aijk + use yaeos__extra_fluids, only: co2_h2o_isop, FixtureFluid + real(pr) :: k(nc, nc, nc), dkdt(nc, nc, nc), dkdt2(nc, nc, nc) + type(FixtureFluid) :: fluid + + integer :: i, j, l + real(pr) :: T, dt=1e-3 + real(pr) :: ai(nc), daidt(nc), daidt2(nc) + real(pr) :: a(nc, nc, nc), dadt(nc, nc, nc), dadt2(nc, nc, nc) + real(pr) :: df, df2 + + fluid = co2_h2o_isop() + + k(1, 1, 2) = 0.1 + k(1, 2, 1) = 0.1 + k(2, 1, 1) = 0.1 + + k(1, 2, 3) = 0.3 + k(2, 1, 3) = 0.3 + k(3, 2, 1) = 0.3 + + k(2, 2, 1) = 0.2 + k(2, 1, 2) = 0.2 + k(1, 2, 2) = 0.2 + dkdt = 0 + dkdt2 = 0 + + T = 200 + + associate(model => fluid%ar_model) + select type(model) + type is (CubicEoS) + call model%calculate_attractive_parameters(T, ai, daidt, daidt2) + call CMR_aijk(& + ai, daidt, daidt2, k, dkdt, dkdt2, & + a, dadt, dadt2 & + ) + end select + end associate + + do i=1,nc + do j=1,nc + do l=1,nc + df = (f(T+dT, i, j, l) - f(T-dT, i, j, l))/(2*dt) + + df2 = (f(T+dT, i, j, l) - 2*f(T, i, j, l) + f(T-dT, i, j, l))/(dT**2) + call assert(abs((df - dadt(i, j, l))/df) < 1e-5, "First order analytic and numerical derivatives should be similar") + call assert(abs((df2 - dadt2(i, j, l))/df) < 1e-5, "Second order analytic and numerical derivatives should be similar") + end do + end do + end do + + contains + real(pr) function f(T, i, j, l) + real(pr), intent(in) :: T + integer, intent(in) :: i, j, l + real(pr) :: ai(nc), daidt(nc), daidt2(nc) + real(pr) :: a(nc, nc, nc), dadt(nc, nc, nc), dadt2(nc, nc, nc) + real(pr) :: eps=1e-5_pr + + associate(model => fluid%ar_model) + select type(model) + type is (CubicEoS) + call model%calculate_attractive_parameters(T, ai, daidt, daidt2) + call CMR_aijk(& + ai, daidt, daidt2, k, dkdt, dkdt2, & + a, dadt, dadt2 & + ) + end select + end associate + f = a(i, j, l) + end function f + end subroutine test_aijk_numdiff + + subroutine test_aijk_kexpt_numdiff + use yaeos__models, only: CMRTD + use yaeos__models_ar_cubic_cubic_mixing, only: kijk_exp_tdep + use yaeos__extra_fluids, only: co2_h2o_isop, FixtureFluid + real(pr) :: k(nc, nc, nc) + real(pr) :: k0(nc, nc, nc) + real(pr) :: tref(nc, nc, nc) + type(FixtureFluid) :: fluid + type(CMRTD) :: mixrule + + integer :: i, j, l + real(pr) :: n(nc), T, dt=1e-3 + real(pr) :: ai(nc), daidt(nc), daidt2(nc) + real(pr) :: a(nc, nc, nc), dadt(nc, nc, nc), dadt2(nc, nc, nc) + real(pr) :: df, df2 + + fluid = co2_h2o_isop() + + k = 0 + k0 = 0 + tref = 0 + + k(1, 1, 2) = 0.1 + k(1, 2, 1) = 0.1 + k(2, 1, 1) = 0.1 + + k(1, 2, 3) = 0.3 + k(2, 1, 3) = 0.3 + k(3, 2, 1) = 0.3 + + k(2, 2, 1) = 0.2 + k(2, 1, 2) = 0.2 + k(1, 2, 2) = 0.2 + + k0(1, 1, 2) = 0.01 + k0(1, 2, 1) = 0.01 + k0(2, 1, 1) = 0.01 + + k0(1, 2, 3) = 0.03 + k0(2, 1, 3) = 0.03 + k0(3, 2, 1) = 0.03 + + k0(2, 2, 1) = 0.02 + k0(2, 1, 2) = 0.02 + k0(1, 2, 2) = 0.02 + + tref(1, 1, 2) = 100 + tref(1, 2, 1) = 100 + tref(2, 1, 1) = 100 + + tref(1, 2, 3) = 300 + tref(2, 1, 3) = 300 + tref(3, 2, 1) = 300 + + tref(2, 2, 1) = 200 + tref(2, 1, 2) = 200 + tref(1, 2, 2) = 200 + + mixrule = CMRTD(k=k, k0=k0, tref=tref, l=0*k0) + + T = 200 + + associate(model => fluid%ar_model) + select type(model) + type is (CubicEoS) + call model%calculate_attractive_parameters(T, ai, daidt, daidt2) + call kijk_exp_tdep(mixrule, & + T, ai, daidt, daidt2, & + a, dadt, dadt2 & + ) + end select + end associate + + do i=1,nc + do j=1,nc + do l=1,nc + df = (f(T+dT, i, j, l) - f(T-dT, i, j, l))/(2*dt) + df2 = (f(T+dT, i, j, l) - 2*f(T, i, j, l) + f(T-dT, i, j, l))/(dT**2) + call assert(abs((df - dadt(i, j, l))/df) < 1e-5, "First order analytic and numerical derivatives should be similar") + call assert(abs((df2 - dadt2(i, j, l))/df) < 1e-5, "Second order analytic and numerical derivatives should be similar") + end do + end do + end do + + contains + real(pr) function f(T, i, j, l) + real(pr), intent(in) :: T + integer, intent(in) :: i, j, l + real(pr) :: eps=1e-5_pr + real(pr) :: ai(nc), daidt(nc), daidt2(nc) + real(pr) :: a(nc, nc, nc), dadt(nc, nc, nc), dadt2(nc, nc, nc) + + associate(model => fluid%ar_model) + select type(model) + type is (CubicEoS) + call model%calculate_attractive_parameters(T, ai, daidt, daidt2) + call kijk_exp_tdep(mixrule, & + T, ai, daidt, daidt2, & + a, dadt, dadt2 & + ) + end select + end associate + + f = a(i, j, l) + end function f + end subroutine test_aijk_kexpt_numdiff + + subroutine test_D_kexpt_numdiff + use yaeos__models, only: CMRTD, CMR + use yaeos__models_ar_cubic_cubic_mixing, only: kijk_exp_tdep + use yaeos__extra_fluids, only: co2_h2o_isop, FixtureFluid + real(pr) :: k(nc, nc, nc) + real(pr) :: k0(nc, nc, nc) + real(pr) :: tref(nc, nc, nc) + type(FixtureFluid) :: fluid + type(CMR) :: mixrule + + integer :: i, j, l + real(pr) :: n(nc), V=0, T, dt=1e-3, dni=1e-9, dn(nc) + real(pr) :: ai(nc), daidt(nc), daidt2(nc) + real(pr) :: a(nc, nc, nc), dadt(nc, nc, nc), dadt2(nc, nc, nc) + real(pr) :: df, df2 + real(pr) :: f1, f2, f3, f4 + + real(pr) :: D + real(pr) :: dDdV + real(pr) :: dDdT + real(pr) :: dDdT2 + real(pr) :: dDdV2 + real(pr) :: dDdTV + real(pr) :: dDi(nc), dDi_num(nc) + real(pr) :: dDidV(nc) + real(pr) :: dDidT(nc) + real(pr) :: dDij(nc, nc), dDij_num(nc, nc) + + fluid = co2_h2o_isop() + + k = 0 + k0 = 0 + tref = 0 + + k(1, 1, 2) = 0.1 + k(1, 2, 1) = 0.1 + k(2, 1, 1) = 0.1 + + k(1, 2, 3) = 0.3 + k(2, 1, 3) = 0.3 + k(3, 2, 1) = 0.3 + + k(2, 2, 1) = 0.2 + k(2, 1, 2) = 0.2 + k(1, 2, 2) = 0.2 + + k0(1, 1, 2) = 0.01 + k0(1, 2, 1) = 0.01 + k0(2, 1, 1) = 0.01 + + k0(1, 2, 3) = 0.03 + k0(2, 1, 3) = 0.03 + k0(3, 2, 1) = 0.03 + + k0(2, 2, 1) = 0.02 + k0(2, 1, 2) = 0.02 + k0(1, 2, 2) = 0.02 + + tref(1, 1, 2) = 100 + tref(1, 2, 1) = 100 + tref(2, 1, 1) = 100 + + tref(1, 2, 3) = 300 + tref(2, 1, 3) = 300 + tref(3, 2, 1) = 300 + + tref(2, 2, 1) = 200 + tref(2, 1, 2) = 200 + tref(1, 2, 2) = 200 + + n = fluid%z0 + + ! mixrule = CMRTD(k=k, k0=k0, tref=tref, l=0*k0) + mixrule = CMR(k=0*k, l=0*k0) + + T = 200 + + associate(model => fluid%ar_model) + select type(model) + type is (CubicEoS) + call model%set_mixrule(mixrule) + call model%calculate_attractive_parameters(T, ai, daidt, daidt2) + call model%mixrule%Dmix(& + n, V, T, & + ai, daidt, daidt2, & + D=D, & + dDdV=dDdV, & + dDdT=dDdT, & + dDdT2=dDdT2, & + dDdV2=dDdV2, & + dDdTV=dDdTV, & + dDi=dDi, & + dDidV=dDidV, & + dDidT=dDidT, & + dDij=dDij & + ) + end select + end associate + + df = (f(n, T+dT) - f(n, T-dT))/(2*dt) + df2 = (f(n, T+dT) - 2*f(n, T) + f(n, T-dT))/(dT**2) + + do i=1,nc + dn = 0 + dn(i) = dni + dDi_num(i) = (f(n + dn, T) - f(n - dn, T))/(2*dni) + do j=1,nc + dn = 0 + dn(i) = dni + dn(j) = dni + f1 = f(n + dn, T) + + dn = 0 + dn(i) = dni + dn(j) = -dni + f2 = f(n + dn, T) + + dn = 0 + dn(i) = -dni + dn(j) = dni + f3 = f(n + dn, T) + + dn = 0 + dn(i) = dni + dn(j) = dni + f4 = f(n + dn, T) + + dDij_num(i, j) = (f1 - f2 - f3 + f4) / (2*dni**2) + + end do + end do + + print *, "" + + print *, dDdT, df + print *, dDdT2, df2 + print *, dDi + print *, dDi_num + + print *, dDij + print *, dDij_num + + contains + real(pr) function f(n, T) + real(pr), intent(in) :: n(nc) + real(pr), intent(in) :: T + real(pr) :: ai(nc), daidt(nc), daidt2(nc) + real(pr) :: a(nc, nc, nc), dadt(nc, nc, nc), dadt2(nc, nc, nc) + + real(pr) :: D !! Mixture attractive parameter \(n^2a_{mix}\) + real(pr) :: dDdV !! \(\frac{dD}{dT}\) + real(pr) :: dDdT !! \(\frac{dD}{dV}\) + real(pr) :: dDdT2 !! \(\frac{d^2D}{dT^2}\) + real(pr) :: dDdV2 !! \(\frac{d^2D}{dV^2}\) + real(pr) :: dDdTV !! \(\frac{d^2D}{dTV\) + real(pr) :: dDi(nc) !! \(\frac{dD}{dn_i}\) + real(pr) :: dDidV(nc) !! \(\frac{d^2D}{dVn_i}\) + real(pr) :: dDidT(nc) !! \(\frac{d^2D}{dTn_i}\) + real(pr) :: dDij(nc, nc)!! \(\frac{d^2D}{dn_{ij}}\) + + associate(model => fluid%ar_model) + select type(model) + type is (CubicEoS) + call model%calculate_attractive_parameters(T, ai, daidt, daidt2) + call model%mixrule%Dmix(& + n, V, T, & + ai, daidt, daidt2, & + D=f, & + dDdV=dDdV, & + dDdT=dDdT, & + dDdT2=dDdT2, & + dDdV2=dDdV2, & + dDdTV=dDdTV, & + dDi=dDi, & + dDidV=dDidV, & + dDidT=dDidT, & + dDij=dDij & + ) + end select + end associate + end function f + end subroutine test_D_kexpt_numdiff + + subroutine test_compare_with_qmr + use yaeos, only: QMR, CMR + use yaeos__extra_fluids, only: co2_h2o_isop, FixtureFluid + type(FixtureFluid) :: fluid + + real(pr) :: k(nc, nc, nc), l(nc, nc, nc) + real(pr) :: n(3), V, T + real(pr) :: D_qmr + real(pr) :: dDdV_qmr, dDdT_qmr, dDdV2_qmr, dDdT2_qmr, dDi_qmr, dDdTV_qmr + real(pr) :: dDidV_qmr(3), dDidT_qmr(3), dDij_qmr(3,3) + + real(pr) :: dDdV_cmr, dDdT_cmr, dDdV2_cmr, dDdT2_cmr, dDi_cmr, dDdTV_cmr + real(pr) :: dDidV_cmr(3), dDidT_cmr(3), dDij_cmr(3,3) + + fluid = co2_h2o_isop() + + T = 200 + V = 2 + n = fluid%z0 + + k = 0 + + associate(model => fluid%ar_model) + select type(model) + type is (CubicEoS) + associate (mr => model%mixrule) + select type (mr) + type is (QMR) + + k(1, 1, 2) = mr%k(1, 2) + k(1, 2, 1) = k(1, 1, 2) + k(1, 2, 1) = k(1, 1, 2) + + k(1, 1, 3) = mr%k(1, 3) + k(1, 3, 1) = k(1, 1, 3) + k(1, 3, 1) = k(1, 1, 3) + + k(2, 2, 3) = mr%k(2, 3) + k(2, 3, 2) = k(2, 2, 3) + k(2, 3, 2) = k(2, 2, 3) + + l(1, 1, 2) = mr%l(1, 2) + l(1, 2, 1) = l(1, 1, 2) + l(1, 2, 1) = l(1, 1, 2) + + l(1, 1, 3) = mr%l(1, 3) + l(1, 3, 1) = l(1, 1, 3) + l(1, 3, 1) = l(1, 1, 3) + + l(2, 2, 3) = mr%l(2, 3) + l(2, 3, 2) = l(2, 2, 3) + l(2, 3, 2) = l(2, 2, 3) + + end select + end associate + end select + end associate + + end subroutine test_compare_with_qmr +end module tests_cmr + +program test_cmr + use tests_cmr, only: test_aijk_numdiff + use tests_cmr, only: test_aijk_kexpt_numdiff + use tests_cmr, only: test_D_kexpt_numdiff + + use testing_aux, only: test_title + + print *, test_title("Cubic Mixing Rules") + call test_aijk_numdiff + call test_aijk_kexpt_numdiff + call test_D_kexpt_numdiff +end program test_cmr