Skip to content

Commit

Permalink
merged RKPR mix
Browse files Browse the repository at this point in the history
  • Loading branch information
@fedebenelli
fedebenelli committed Dec 27, 2024
1 parent f1cfe8c commit 26dca1a
Show file tree
Hide file tree
Showing 6 changed files with 158 additions and 67 deletions.
14 changes: 8 additions & 6 deletions src/models/residual_helmholtz/cubic/implementations/implementations.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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)
Expand All @@ -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
Expand Down
30 changes: 3 additions & 27 deletions src/models/residual_helmholtz/cubic/mixing_rules/quadratic_mixing.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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(:, :)
Expand Down Expand Up @@ -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.
Expand All @@ -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
Expand Down Expand Up @@ -287,16 +265,14 @@ 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))
aij_hd(j, i) = aij_hd(i, j)
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
Expand Down
23 changes: 12 additions & 11 deletions test/test_cubic_implementations.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down
45 changes: 24 additions & 21 deletions test/test_implementations/ar_models/cubics/test_rkpr.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand All @@ -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 &
Expand All @@ -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)
Expand Down Expand Up @@ -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)
Expand Down
4 changes: 2 additions & 2 deletions test/test_models_ar_cubic_quadrating_mixing.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down
109 changes: 109 additions & 0 deletions test/test_qmrtd.f90
View file
Original file line number Diff line number Diff line change
@@ -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

0 comments on commit 26dca1a

Please sign in to comment.