module yaeos__models_ar_gerg2008 use yaeos__constants, only: Ryaeos => R, pr !! Ideal gas constants used on yaeos use yaeos__adiff_hyperdual_ar_api, only: ArModelAdiff use yaeos__models_ar_cubic_implementations, only: SoaveRedlichKwong use yaeos__models_ar_genericcubic, only: CubicEoS use yaeos__models_ar_multifluid_parameters_gerg2008, only: Gerg2008Binary, Gerg2008Pure use hyperdual_mod implicit none type, extends(ArModelAdiff) :: Gerg2008 type(Gerg2008Pure), allocatable :: pures(:) type(Gerg2008Binary), allocatable :: binaries(:, :) type(CubicEoS) :: srk contains procedure :: ar => arfun procedure :: get_v0 => volume_initalizer procedure :: volume => volume end type Gerg2008 type, private :: GERG2008Selector integer :: methane=1 integer :: nitrogen=2 integer :: carbon_dioxide=3 integer :: ethane=4 integer :: propane=5 integer :: nbutane=6 integer :: isobutane=7 integer :: npentane=8 integer :: isopentane=9 integer :: nhexane=10 integer :: nheptane=11 integer :: noctane=12 integer :: nonane=13 integer :: decane=14 integer :: hydrogen=15 integer :: oxygen=16 integer :: carbon_monoxide=17 integer :: water=18 integer :: hydrogen_sulfide=19 integer :: helium=20 integer :: argon=21 end type GERG2008Selector type(GERG2008Selector) :: G2008Components contains type(Gerg2008) function gerg_2008(ids) use yaeos__models_ar_multifluid_parameters_gerg2008, only: get_original_parameters integer, intent(in) :: ids(:) type(Gerg2008Pure) :: pures(size(ids)) type(Gerg2008Binary) :: binaries(size(ids), size(ids)) call get_original_parameters(ids, pures, binaries, gerg_2008%components) gerg_2008%pures = pures gerg_2008%binaries = binaries gerg_2008%srk =SoaveRedlichKwong(gerg_2008%components%Tc, gerg_2008%components%Pc, gerg_2008%components%w) end function gerg_2008 subroutine reducing_functions(self, n, Vr, Tr) class(Gerg2008), intent(in) :: self type(hyperdual), intent(in) :: n(:) type(hyperdual), intent(out) :: Vr type(hyperdual), intent(out) :: Tr type(hyperdual) :: X(size(n)) real(8) :: Vc(size(n)), Tc(size(n)), rho_c(size(n)) real(8) :: Bv(size(n), size(n)), Gv(size(n), size(n)) real(8) :: Bt(size(n), size(n)), Gt(size(n), size(n)) integer :: i, j, nc Vc = self%components%Vc Tc = self%components%Tc Bv = self%binaries%Bv Gv = self%binaries%Gv Bt = self%binaries%Bt Gt = self%binaries%Gt rho_c = 1/Vc X = n / sum(n) nc = size(n) Vr = sum(X ** 2 * Vc) Tr = sum(X ** 2 * Tc) do i=1,nc do j=i+1,nc Vr = Vr + & 2 * X(i) * X(j) * Bv(i, j) * Gv(i, j) & * (X(i) + X(j)) / (Bv(i, j) ** 2 * X(i) + X(j)) & * 1._pr / 8._pr * (rho_c(i) ** (- 1._pr / 3._pr) & + rho_c(j) ** (- 1._pr / 3)) ** 3 Tr = Tr + & 2 * X(i) * X(j) * Bt(i, j) * Gt(i, j) & * (X(i) + X(j)) / (Bt(i, j) ** 2 * X(i) + X(j)) & * sqrt((Tc(i) * Tc(j))) end do end do end subroutine reducing_functions subroutine ar_pure(pure, delta, tau, ar) type(Gerg2008Pure), intent(in) :: pure type(hyperdual), intent(in) :: delta type(hyperdual), intent(in) :: tau type(hyperdual), intent(out) :: ar integer :: i, Kpol, Kexp real(8) :: n_pol(pure%Kpol), d_pol(pure%Kpol), t_pol(pure%Kpol) real(8) :: n_exp(pure%Kexp), d_exp(pure%Kexp), t_exp(pure%Kexp) real(8) :: c_exp(pure%Kexp) Kpol = pure%Kpol Kexp = pure%Kexp n_pol = pure%n(1:Kpol) d_pol = pure%d(1:Kpol) t_pol = pure%t(1:Kpol) n_exp = pure%n(Kpol+1:Kpol+Kexp) d_exp = pure%d(Kpol+1:Kpol+Kexp) t_exp = pure%t(Kpol+1:Kpol+Kexp) c_exp = pure%c ar = sum(n_pol * delta ** d_pol * tau ** t_pol) + & sum(n_exp * delta**d_exp * tau**t_exp * exp(-delta**c_exp)) end subroutine ar_pure subroutine ar_ij(delta, tau, binary, aij) type(hyperdual), intent(in) :: delta type(hyperdual), intent(in) :: tau type(Gerg2008Binary), intent(in) :: binary type(hyperdual), intent(out) :: aij integer :: idx_poly, idx_exp real(pr) :: n_pol(binary%Kpolij), d_pol(binary%Kpolij), t_pol(binary%Kpolij) real(pr) :: n_exp(binary%Kexpij), d_exp(binary%Kexpij), t_exp(binary%Kexpij) real(pr) :: etha(binary%Kexpij), eps(binary%Kexpij), beta(binary%Kexpij), gama(binary%Kexpij) idx_poly = binary%Kpolij idx_exp = binary%Kexpij + idx_poly n_pol = binary%nij(1:idx_poly) d_pol = binary%dij(1:idx_poly) t_pol = binary%tij(1:idx_poly) n_exp = binary%nij(idx_poly+1:idx_exp) d_exp = binary%dij(idx_poly+1:idx_exp) t_exp = binary%tij(idx_poly+1:idx_exp) etha = binary%ethaij(1:binary%Kexpij) eps = binary%epsij(1:binary%Kexpij) beta = binary%betaij(1:binary%Kexpij) gama = binary%gammaij(1:binary%Kexpij) aij = sum(n_pol * delta ** d_pol * tau ** t_pol) + & sum(n_exp * delta**d_exp * tau**t_exp * exp(-etha * (delta - eps) ** 2 - beta * (delta - gama))) end subroutine ar_ij function arfun(self, n, v, t) result(arval) class(Gerg2008) :: self type(hyperdual), intent(in) :: n(:), v, t type(hyperdual) :: arval type(hyperdual) :: Vr, Tr, X(size(n)), rho_r type(hyperdual) :: delta, tau type(hyperdual) :: aij type(hyperdual) :: ar_pures(size(n)) type(Gerg2008Pure) :: pures(size(n)) type(Gerg2008Binary) :: binary real(pr) :: rho_c(size(n)) real(pr) :: Fij(size(n), size(n)) integer :: i, j, nc Fij = self%binaries%Fij pures = self%pures nc = size(n) X = n / sum(n) call reducing_functions(self, n, Vr, Tr) rho_r = 1._pr/Vr delta = (1._pr/(V/sum(n)))/rho_r tau = Tr/T do i=1,nc call ar_pure(pures(i), delta, tau, ar_pures(i)) end do arval = sum(x * ar_pures) do i=1,nc do j=1,nc!i+1,nc if (Fij(i, j) == 0._pr) cycle binary = self%binaries(i, j) call ar_ij(delta, tau, binary, aij) arval = arval + X(i) * X(j) * Fij(i, j) * aij end do end do arval = arval * (sum(n) * ryaeos * t) end function arfun subroutine volume(eos, n, P, T, V, root_type) !! Volume solver routine for the GERG2008. !! !! Solves volume roots using newton method. Given pressure and !! temperature. It will use the SRK equation of state to initialize the !! volume values. !! !! # Description !! This subroutine solves the volume using a newton method. The variable !! `root_type` is used to specify the desired root to solve. The options !! are: `["liquid", "vapor", "stable"]` !! !! # Examples !! !! ```fortran !! eos = PengRobinson76(Tc, Pc, w) !! !! n = [1.0_pr, 1.0_pr] !! T = 300.0_pr !! P = 1.0_pr !! !! call eos%volume(n, P, T, V, root_type="liquid") !! call eos%volume(n, P, T, V, root_type="vapor") !! call eos%volume(n, P, T, V, root_type="stable") !! ``` use yaeos__constants, only: pr, R use yaeos__math, only: newton class(GERG2008), intent(in) :: eos !! Model real(pr), intent(in) :: n(:) !! Moles number vector real(pr), intent(in) :: P !! Pressure [bar] real(pr), intent(in) :: T !! Temperature [K] real(pr), intent(out) :: V !! Volume [L] character(len=*), intent(in) :: root_type !! Desired root-type to solve. Options are: !! `["liquid", "vapor", "stable"]` integer :: max_iters=30 real(pr) :: tol=1e-8 real(pr) :: totnRT, GrL, GrV, Gr real(pr) :: Vliq, Vvap logical :: failed GrL = HUGE(GrL) GrV = HUGE(GrV) totnRT = sum(n) * R * T select case(root_type) case("liquid") Vliq = eos%get_v0(n, P, T) * 1.001_pr call newton(foo, Vliq, tol=tol, max_iters=max_iters, failed=failed) GrL = Gr case("vapor") call eos%srk%volume(n=n, P=P, T=T, V=Vvap, root_type="vapor") call newton(foo, Vvap, tol=tol, max_iters=max_iters, failed=failed) GrV = Gr case("stable") Vliq = eos%get_v0(n, P, T)*1.00001_pr call newton(foo, Vliq, tol=tol, max_iters=max_iters, failed=failed) GrL = Gr Vvap = R * T / P call newton(foo, Vvap, tol=tol, max_iters=max_iters, failed=failed) GrV = Gr end select if (GrL < GrV) then V = Vliq else V = Vvap end if if (failed) V = -1 contains subroutine foo(x, f, df) real(pr), intent(in) :: x real(pr), intent(out) :: f, df real(pr) :: Ar, ArV, ArV2, Pcalc, dPcalcdV, Vin Vin = x call eos%residual_helmholtz(n, Vin, T, Ar=Ar, ArV=ArV, ArV2=ArV2) Pcalc = totnRT / Vin - ArV dPcalcdV = -totnRT / Vin**2 - ArV2 f = Pcalc - p df = dPcalcdV Gr = Ar + P * Vin - totnRT - totnRT * log(P*Vin/(R*T)) end subroutine foo end subroutine volume function volume_initalizer(self, n, p, t) result(v0) class(Gerg2008), intent(in) :: self real(pr), intent(in) :: n(:) real(pr), intent(in) :: p real(pr), intent(in) :: t real(pr) :: v0 v0 = self%srk%get_v0(n, p, t) end function volume_initalizer end module yaeos__models_ar_gerg2008