From 26dca1a9feca8cedbee821e0786036764a6ff481 Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Fri, 27 Dec 2024 18:04:45 -0300 Subject: [PATCH] merged RKPR mix --- .../cubic/implementations/implementations.f90 | 14 ++- .../cubic/mixing_rules/quadratic_mixing.f90 | 30 +---- test/test_cubic_implementations.f90 | 23 ++-- .../ar_models/cubics/test_rkpr.f90 | 45 ++++---- ...test_models_ar_cubic_quadrating_mixing.f90 | 4 +- test/test_qmrtd.f90 | 109 ++++++++++++++++++ 6 files changed, 158 insertions(+), 67 deletions(-) create mode 100644 test/test_qmrtd.f90 diff --git a/src/models/residual_helmholtz/cubic/implementations/implementations.f90 b/src/models/residual_helmholtz/cubic/implementations/implementations.f90 index 494e550c..e4d614e7 100644 --- a/src/models/residual_helmholtz/cubic/implementations/implementations.f90 +++ b/src/models/residual_helmholtz/cubic/implementations/implementations.f90 @@ -295,7 +295,7 @@ type(CubicEoS) function RKPR(tc, pc, w, zc, kij, lij, delta_1, k) result(model) !! !! \[\delta_1 = d_1 + d_2 (d_3 - Z_c)^d_4 + d_5 (d_3 - Z_c) ^ d_6\] !! \[k = (A_1 Z_c + A_0)\omega^2 + (B_1 Z_c + B_0)\omega + (C_1 Z_c + C_0)\] - use yaeos__models_ar_cubic_quadratic_mixing, only: QMR_RKPR + use yaeos__models_ar_cubic_quadratic_mixing, only: QMR use yaeos__models_ar_cubic_alphas, only: AlphaRKPR real(pr), intent(in) :: tc(:) !! Critical Temperature [K] real(pr), intent(in) :: pc(:) !! Critical Pressure [bar] @@ -307,7 +307,7 @@ type(CubicEoS) function RKPR(tc, pc, w, zc, kij, lij, delta_1, k) result(model) real(pr), optional, intent(in) :: k(:) type(AlphaRKPR) :: alpha - type(QMR_RKPR) :: mixrule + type(QMR) :: mixrule type(Substances) :: composition integer :: i, nc @@ -339,9 +339,9 @@ type(CubicEoS) function RKPR(tc, pc, w, zc, kij, lij, delta_1, k) result(model) Zc_eos = 1.168 * Zc if (present(k)) then - alpha%k = k + alpha%k = k else - alpha%k = (A1 * zc + A0)*w**2 + (B1*zc + B0)*w + (C1*Zc + C0) + alpha%k = (A1 * Zc_eos + A0)*w**2 + (B1*zc + B0)*w + (C1*Zc_eos + C0) end if if (present(kij)) then @@ -358,9 +358,9 @@ type(CubicEoS) function RKPR(tc, pc, w, zc, kij, lij, delta_1, k) result(model) model%components = composition if (present(delta_1)) then - model%del1 = delta_1 + model%del1 = delta_1 else - model%del1 = d1 + d2 * (d3 - zc) ** d4 + d5 * (d3 - zc) ** d6 + model%del1 = d1 + d2 * (d3 - zc_eos) ** d4 + d5 * (d3 - zc_eos) ** d6 end if model%del2 = (1._pr - model%del1)/(1._pr + model%del1) @@ -380,6 +380,8 @@ type(CubicEoS) function RKPR(tc, pc, w, zc, kij, lij, delta_1, k) result(model) Psat_i = model%Psat_pure(i, 0.7*Tc(i)) diff = (w(i) - (-1 - log10(Psat_i/Pc(i)))) alpha%k(i) = alpha%k(i) + 0.1*diff + + deallocate(model%alpha) model%alpha = alpha end do end do diff --git a/src/models/residual_helmholtz/cubic/mixing_rules/quadratic_mixing.f90 b/src/models/residual_helmholtz/cubic/mixing_rules/quadratic_mixing.f90 index 1432b957..568ef2d8 100644 --- a/src/models/residual_helmholtz/cubic/mixing_rules/quadratic_mixing.f90 +++ b/src/models/residual_helmholtz/cubic/mixing_rules/quadratic_mixing.f90 @@ -24,12 +24,8 @@ module yaeos__models_ar_cubic_quadratic_mixing procedure :: aij => kij_constant !! Default attractive parameter combining rule procedure :: Dmix !! Attractive parameter mixing rule procedure :: Bmix !! Repulsive parameter mixing rule - procedure :: D1mix => D1mix_constant + procedure :: D1mix => RKPR_D1mix end type QMR - type, extends(QMR) :: QMR_RKPR - contains - procedure :: D1Mix => RKPR_D1mix - end type QMR_RKPR type, extends(QMR) :: QMRTD real(pr), allocatable :: k0(:, :) @@ -147,24 +143,6 @@ subroutine Bmix(self, n, bi, B, dBi, dBij) call bmix_qmr(n, bi, self%l, b, dbi, dbij) end subroutine Bmix - subroutine D1mix_constant(self, n, d1i, D1, dD1i, dD1ij) - !! Constant \(\delta_1\) parameter. - !! - !! Most Cubic EoS keep a constant value for their \(\delta_1\) parameter. - !! This procedure assumes that all the components have the same \(delta_1\) - !! and takes the first value as the one of the mixture. - use yaeos__models_ar_cubic_mixing_base, only: d1mix_rkpr - class(QMR), intent(in) :: self !! Mixing rule - real(pr), intent(in) :: n(:) !! Moles vector - real(pr), intent(in) :: d1i(:) !! \(\delta_1\) parameter - real(pr), intent(out) :: D1 !! Mixture's \(\Delta_1\) - real(pr), intent(out) :: dD1i(:) !! \(\frac{dDelta_1}{dn_i} = 0\) - real(pr), intent(out) :: dD1ij(:, :) !! \(\frac{d^2Delta_1}{dn_{ij}} = 0\) - D1 = d1i(1) - dD1i = 0 - dD1ij = 0 - end subroutine D1mix_constant - 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. @@ -177,7 +155,7 @@ subroutine RKPR_D1mix(self, n, d1i, D1, dD1i, dD1ij) !! \Delta_1 = \sum_i^N n_i \delta_{1i} !! \] !! - class(QMR_RKPR), intent(in) :: self + class(QMR), intent(in) :: self real(pr), intent(in) :: n(:) real(pr), intent(in) :: d1i(:) real(pr), intent(out) :: D1 @@ -287,7 +265,7 @@ subroutine kij_exp_tdep(& nc = size(a) - do i=1,size(a)-1 + do i=1,size(a) aij_hd(i, i) = sqrt(a_hd(i) * a_hd(i)) do j=i+1,size(a) aij_hd(i, j) = sqrt(a_hd(i) * a_hd(j)) * (1._pr - kij_hd(i, j)) @@ -295,8 +273,6 @@ subroutine kij_exp_tdep(& end do end do - aij_hd(size(a), size(a)) = sqrt(a_hd(size(a)) * a_hd(size(a))) - aij = aij_hd%f0 daijdt = aij_hd%f1 daijdt2 = aij_hd%f12 diff --git a/test/test_cubic_implementations.f90 b/test/test_cubic_implementations.f90 index e2f28dcc..4e4036a9 100644 --- a/test/test_cubic_implementations.f90 +++ b/test/test_cubic_implementations.f90 @@ -179,17 +179,18 @@ subroutine test_RKPR(error) real(pr) :: Ar_val, ArV_val, ArV2_val, ArT_val, ArTV_val, ArT2_val real(pr) :: Arn_val(n), ArVn_val(n), ArTn_val(n), Arn2_val(n, n) - Ar_val = -9.6482477518604082 - ArV_val = 8.5009523009964276 - ArT_val = 1.3361597729919862E-002 - ArT2_val = -1.9770475708937217E-005 - ArV2_val = -15.089728149827863 - ArTV_val = -1.2540180323248612E-002 - Arn_val = [-15.232342932903979, -19.399281746656321] - ArVn_val = [12.154119716077535, 16.347846121425317] - ArTn_val = [2.0758309799901429E-002, 2.8106121715415409E-002] - Arn2_val(1, :) = [-11.397459730966577, -12.478402585247725] - Arn2_val(2, :) = [-12.478402585247725, -18.006179159638897] + Ar_val = -9.5995556136857765 + ArV_val = 8.7952184911116689 + ArT_val = 2.1732218215585023E-002 + ArT2_val = -4.7480350864298389E-005 + ArV2_val = -16.171996177263189 + ArTV_val = -2.0750113033467819E-002 + Arn_val = [-15.106022804526534, -19.804239027780842] + ArVn_val = [ 12.785101276300841, 17.623522645787737] + ArTn_val = [ 3.6168578086399804E-002, 4.5188225615222108E-002] + Arn2_val(1, :) = [-11.652771433485331, -13.270385522198149] + Arn2_val(2, :) = [-13.270385522198149, -19.489152947516825] + eos = binary_RKPR() z = [0.3, 0.7] diff --git a/test/test_implementations/ar_models/cubics/test_rkpr.f90 b/test/test_implementations/ar_models/cubics/test_rkpr.f90 index f5256e9f..ae446db2 100644 --- a/test/test_implementations/ar_models/cubics/test_rkpr.f90 +++ b/test/test_implementations/ar_models/cubics/test_rkpr.f90 @@ -34,6 +34,7 @@ subroutine test_rkpr_cons_mixture(error) real(pr) :: eq31, eq33(size(n), size(n)), eq34(size(n)), eq36, eq37 real(pr) :: kij(size(n), size(n)), lij(size(n), size(n)) + real(pr) :: tol n = [1.5, 0.2, 0.7, 2.3] tc = [369.83, 425.12, 507.6, 540.2] @@ -54,7 +55,7 @@ subroutine test_rkpr_cons_mixture(error) ArTV=ArTV, ArV2=ArV2, ArT2=ArT2, ArVn=ArVn, ArTn=ArTn, Arn2=Arn2) call numeric_ar_derivatives(& - model, n, v, t, d_n = 0.0001_pr, d_v = 0.0001_pr, d_t = 0.01_pr, & + model, n, v, t, d_n = 0.0001_pr, d_v = 0.00001_pr, d_t = 0.001_pr, & Ar=Ar_num, ArV=ArV_num, ArT=ArT_num, ArTV=ArTV_num, ArV2=ArV2_num, & ArT2=ArT2_num, Arn=Arn_num, ArVn=ArVn_num, ArTn=ArTn_num, & Arn2=Arn2_num & @@ -64,17 +65,19 @@ subroutine test_rkpr_cons_mixture(error) model, n, v, t, eq31=eq31, eq33=eq33, eq34=eq34, eq36=eq36, eq37=eq37 & ) + tol = 1e-4 + ! Numeric derivatives - call check(error, rel_error(Ar, Ar_num) < 1e-5) - call check(error, rel_error(ArV, ArV_num) < 1e-5) - call check(error, rel_error(ArT, ArT_num) < 1e-5) - call check(error, allclose(Arn, Arn_num, 1e-5_pr)) - call check(error, rel_error(ArV2, ArV2_num) < 1e-5) - call check(error, rel_error(ArT2, ArT2_num) < 1e-5) - call check(error, rel_error(ArTV, ArTV_num) < 1e-5) - call check(error, allclose(ArVn, ArVn_num, 1e-5_pr)) - call check(error, allclose(ArTn, ArTn_num, 1e-5_pr)) - call check(error, maxval(rel_error(Arn2, Arn2_num)) < 1e-5) + call check(error, rel_error(Ar, Ar_num) < tol) + call check(error, rel_error(ArV, ArV_num) < tol) + call check(error, rel_error(ArT, ArT_num) < tol) + call check(error, allclose(Arn, Arn_num, tol)) + call check(error, rel_error(ArV2, ArV2_num) < tol) + call check(error, rel_error(ArT2, ArT2_num) < tol) + call check(error, rel_error(ArTV, ArTV_num) < tol) + call check(error, allclose(ArVn, ArVn_num, tol)) + call check(error, allclose(ArTn, ArTn_num, tol)) + call check(error, maxval(rel_error(Arn2, Arn2_num)) < tol) ! Consistency tests call check(error, abs(eq31) <= 1e-13) @@ -118,16 +121,16 @@ subroutine test_rkpr_cons_mixture(error) ) ! Numeric derivatives - call check(error, rel_error(Ar, Ar_num) < 1e-5) - call check(error, rel_error(ArV, ArV_num) < 1e-5) - call check(error, rel_error(ArT, ArT_num) < 1e-5) - call check(error, allclose(Arn, Arn_num, 1e-5_pr)) - call check(error, rel_error(ArV2, ArV2_num) < 1e-5) - call check(error, rel_error(ArT2, ArT2_num) < 1e-5) - call check(error, rel_error(ArTV, ArTV_num) < 1e-5) - call check(error, allclose(ArVn, ArVn_num, 1e-5_pr)) - call check(error, allclose(ArTn, ArTn_num, 1e-5_pr)) - call check(error, maxval(rel_error(Arn2, Arn2_num)) < 1e-5) + call check(error, rel_error(Ar, Ar_num) < tol) + call check(error, rel_error(ArV, ArV_num) < tol) + call check(error, rel_error(ArT, ArT_num) < tol) + call check(error, allclose(Arn, Arn_num, tol)) + call check(error, rel_error(ArV2, ArV2_num) < tol) + call check(error, rel_error(ArT2, ArT2_num) < tol) + call check(error, rel_error(ArTV, ArTV_num) < tol) + call check(error, allclose(ArVn, ArVn_num, tol)) + call check(error, allclose(ArTn, ArTn_num, tol)) + call check(error, maxval(rel_error(Arn2, Arn2_num)) < tol) ! Consistency tests call check(error, abs(eq31) <= 1e-13) diff --git a/test/test_models_ar_cubic_quadrating_mixing.f90 b/test/test_models_ar_cubic_quadrating_mixing.f90 index c376ccf4..4dd8d4ea 100644 --- a/test/test_models_ar_cubic_quadrating_mixing.f90 +++ b/test/test_models_ar_cubic_quadrating_mixing.f90 @@ -20,10 +20,10 @@ end subroutine collect_suite subroutine test_QMR_RKPR(error) use yaeos__constants, only: pr - use yaeos__models_ar_cubic_quadratic_mixing, only: QMR_RKPR + use yaeos__models_ar_cubic_quadratic_mixing, only: QMR type(error_type), allocatable, intent(out) :: error - type(QMR_RKPR) :: mixrule + type(QMR) :: mixrule integer, parameter :: n = 3 real(pr) :: del1(3) = [0.2, 0.5, 0.6] diff --git a/test/test_qmrtd.f90 b/test/test_qmrtd.f90 new file mode 100644 index 00000000..0d2e6b92 --- /dev/null +++ b/test/test_qmrtd.f90 @@ -0,0 +1,109 @@ +program main + use yaeos + implicit none + type(CubicEoS) :: model + type(QMRTD) :: mixrule + + integer, parameter :: nc=2 + real(pr) :: k0(nc, nc)=0, kinf(nc, nc)=0, Tref(nc, nc)=0 + + real(pr) :: tc(nc), pc(nc), w(nc) + + real(pr) :: D, dDT, dDT2, dDi(nc), dDij(nc, nc), dDidT(nc) + real(pr) :: ai(nc), daidt(nc), daidt2(nc) + + real(pr) :: n(nc), T, V, P, f1, f2, f3, f4, dx=1E-3 + integer :: i, j + + tc = [126.2, 568.7] + pc = [33.98, 24.9] + w = [9.01E-002, 0.492] + + k0(1, 2) = 0!-0.429315E+00 + k0(2, 1) = k0(1, 2) + kinf(1, 2) = 0.485065E+01 + kinf(2, 1) = kinf(1, 2) + + Tref(1, 2) = Tc(1) + Tref(2, 1) = Tref(1, 2) + + mixrule = QMRTD(k=kinf, k0=k0, Tref=Tref, l=0*k0) + + model = SoaveRedlichKwong(Tc, Pc, w, kij=kinf) + n = [0.2, 0.8] + T = 150 + + + print *, "===================================================================" + print *, "NORMAL MR" + call model%alpha%alpha(T/model%components%Tc, ai, daidt, daidt2) + ai = ai*model%ac + daidt = daidt*model%ac/model%components%Tc + daidt2 = daidt2*model%ac/model%components%Tc**2 + call model%mixrule%Dmix(n, T, ai=ai, daidt=daidt, daidt2=daidt2, D=D, dDdT=dDT, dDdT2=dDT2, dDi=dDi, dDidT=dDidT, dDij=dDij) + print *, D, dDT, dDT2 + print *, dDi + print *, dDidT + print *, dDij + print *, "===================================================================" + + + call model%set_mixrule(mixrule) + + call model%alpha%alpha(T/model%components%Tc, ai, daidt, daidt2) + ai = ai*model%ac + daidt = daidt*model%ac/model%components%Tc + daidt2 = daidt2*model%ac/model%components%Tc**2 + + + call model%mixrule%Dmix(n, T, ai=ai, daidt=daidt, daidt2=daidt2, D=D, dDdT=dDT, dDdT2=dDT2, dDi=dDi, dDidT=dDidT, dDij=dDij) + print *, D + + f2 = get_D(n, T+dx) + f1 = get_D(n, T-dx) + + print *, (f2-f1)/(2*dx), dDT + print *, (f2 - 2*D + f1)/dx**2, dDT2 + + do i = 1, nc + n = [0.2, 0.8] + n(i) = n(i) + dx + f2 = get_D(n, T) + n(i) = n(i) - 2*dx + f1 = get_D(n, T) + n(i) = n(i) + dx + print *, (f2-f1)/(2*dx), dDi(i) + end do + + dx = 0.01 + + f1 = get_D(n + [dx, dx], T) + f2 = get_D(n + [dx, -dx], T) + f3 = get_D(n + [-dx, dx], T) + f4 = get_D(n + [-dx, -dx], T) + print *, (f1 - f2 - f3 + f4)/(4*dx**2), dDij(1, 2), dDij(2, 1) + + dx = 0.01 + f1 = get_D(n + [dx, 0._pr], T) + f2 = get_D(n - [dx, 0._pr], T) + print *, (f1 - 2*D + f2)/(dx**2), dDij(1, 1) + + f1 = get_D(n + [0._pr, dx], T) + f2 = get_D(n - [0._pr, dx], T) + print *, (f1 - 2*D + f2)/(dx**2), dDij(2, 2) + + +contains + real(pr) function get_D(n, T) result(D) + real(pr) :: n(nc), T + + real(pr) :: ai(nc), daidt(nc), daidt2(nc) + real(pr) :: dDT, dDT2, dDi(nc), dDij(nc, nc), dDidT(nc) + + call model%alpha%alpha(T/model%components%Tc, ai, daidt, daidt2) + ai = ai*model%ac + daidt = daidt*model%ac/model%components%Tc + daidt2 = daidt2*model%ac/model%components%Tc**2 + call model%mixrule%Dmix(n, T, ai=ai, daidt=daidt, daidt2=daidt2, D=D, dDdT=dDT, dDdT2=dDT2, dDi=dDi, dDidT=dDidT, dDij=dDij) + end function get_D +end program main