In [1]:
import numpy as np
RGAS=0.08314472
In [3]:
RGAS=0.08314472
# Critical constants must be given in K and bar
# b will be in L/mol and ac in bar*(L/mol)**2
# PARAMETER (A0=0.0017,B0=1.9681,C0=-2.7238)
A0 = 0.0017
B0 = 1.9681
C0 =-2.7238
# PARAMETER (A1=-2.4407,B1=7.4513,C1=12.504)
A1 = -2.4407
B1 = 7.4513
C1 = 12.504
# dimension D(6)
# CHARACTER*16 comp
# COMMON /Tcdc/ Tc,dc
# COMMON /ABd1/ a,b,del1
D = np.array([0.428363,18.496215,0.338426,0.660,789.723105,2.512392])
NC=2
NIN=1
nout=9
"""
acá encontré esto que lo tenía aparte, como un ejemplo de ese tipo de corrida con un RHOLSat
3 3 0 ICALC,NMODEL
304.21 73.83 0.2236 Tc, Pc, omega
270.0 21.4626 T(K), RhoLsat (L/mol)
"""
Out[3]:
In [13]:
T = 270.0 #148.0 # 110.0
P = 10.0 # 1.00
# METHANE (1)
# 190.564 45.99 0.01155 0.115165 1.00000173664 tc, pc, ohm, vc, zrat
# 2.3277 0.029962 0.932475 1.49541 ac, b, delta1, k
Tc, Pc, ohm, vc, zrat = 190.564, 45.99, 0.01155, 0.115165, 1.00000173664
ac, b, d, rk = 2.3277, 0.029962, 0.932475, 1.49541
Tr = T / Tc
print("Tr_cal = {0}".format(Tr))
#d = delta1
#dc = 1/vc
OM = ohm
Zc = zrat
T = 270.0 #148.0 # 110.0
P = 10.0 # 1.00
# acá encontré esto que lo tenía aparte, como un ejemplo de ese tipo de corrida con un RHOLSat
# 3 3 0 ICALC,NMODEL
# 304.21 73.83 0.2236 Tc, Pc, omega
# 270.0 21.4626 T(K), RhoLsat (L/mol)
Tc, Pc, omega = 304.21, 73.83, 0.2236
OM = omega
# 270.0 21.4626 T(K), RhoLsat (L/mol)
T = 270.0
RhoLsat = 21.4626
del1 = 2.0
d1 = (1 + del1 ** 2) / (1 + del1)
y = 1 + (2 * (1 + del1)) ** (1.0 / 3) + (4 / (1 + del1)) ** (1.0 / 3)
OMa = (3 * y * y + 3 * y * d1 + d1 ** 2 + d1 - 1.0) / (3 * y + d1 - 1.0) ** 2
OMb = 1 / (3 * y + d1 - 1.0)
Zc = y / (3 * y + d1 - 1.0)
dc = Pc / Zc /(RGAS*Tc)
Vceos = 1.0/dc
rk = (A1 * Zc + A0) * OM**2 + (B1 * Zc + B0) * OM + (C1 * Zc + C0) #! initial guess for k parameter
Tr = 0.7 #D0 ! Change here to use another Pv than the one at Tr 0.7
Pvdat = Pc * 10** -(1.0 + OM)
a = ac * (3/(2 + Tr))**rk
P = Pvdat*1.5
print(ac, rk, a, Pvdat, dc, Vceos)
In [8]:
# SUBROUTINE VaporPressure(Tr,PVini,Pv,RHOL,RHOV,phiL)
def vapor_pressure_cal(Tr,PVini,Pvdat):
"""
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
PARAMETER (ERRMAX=1.D-8)
COMMON /Tcdc/ Tc,dc
"""
dphi = 0.01
Pold = PVini
P_new = Pvdat
n=1
T=Tr*Tc
ITER_p = 0
while 1:
# 30 call VCALC(1,T,P,V)
V_liquido = VCALC(1,T,Pold,4.0)
RHOL = 1/V_liquido
# call VCALC(-1,T,P,V) #! SOLVE for vapor density
V_vapor = VCALC(-1,T,Pold,4.0) #! SOLVE for vapor density
RHOV = 1/V_vapor
print("RHOL = {0}".format(RHOL))
print("RHOV = {0}".format(RHOV))
print("ITER_p = {0}".format(ITER_p))
"""
#if(RHOL.LT.0.9*dc) then
if RHOL < 0.9*dc:
P = 1.01 * P
print("P = 1.01 * P")
# go to 30
# else if(RHOV.GT.dc) then
elif RHOV > dc:
P = 0.99 * P
print("P = 0.99 * P")
# go to 30
# end if
"""
ITER_v = 0
while RHOL < 0.9*dc or RHOV > dc:
V_liquido = VCALC(1,T,Pold,4.0)
RHOL = 1/V_liquido
V_vapor = VCALC(-1,T,Pold,4.0)
RHOV = 1/V_vapor
if RHOL < 0.9*dc:
Pold = 1.01 * Pold
print("P = 1.01 * P")
# go to 30
# else if(RHOV.GT.dc) then
elif RHOV > dc:
Pold = 0.99 * Pold
print("P = 0.99 * P")
print("RHOL = {0} y RHOV = {1}".format(RHOL,RHOV))
print("otra vez P")
if ITER_v > 500:
break
ITER_v += 1
# call FUG_CALC(T,P,1/RHOL,phi)
phi = FUG_CALC(T,Pold,1/RHOL)
phiL = phi#[0]
# call FUG_CALC(T,P,V,phi)
phi = FUG_CALC(T,Pold,1/RHOV)
phiV = phi#[0]
dphiold = dphi
dphi = phiV - phiL
print("phiL = {0}".format(phiL))
print("phiV = {0}".format(phiV))
print("dphi = {0}".format(dphi))
# if(dphiold.eq.0.0D0.or.Tr.gt.0.975) then
if abs(dphiold) == 0.0 or Tr > 0.975:
P_k_new = Pold * (phiL/phiV)
print("P = P * (phiL/phiV)")
print("(phiL/phiV) = {0}".format((phiL/phiV)))
print("Tr = {0}".format(Tr))
# else
else:
#P = Plast - dphi*(Plast - Pold)/(dphi - dphiold)
#P = P - 0.5*dphi*(P - Pold)/(dphi - dphiold)
P_k_new = P_new - dphi*(P_new - Pold)/(dphi - dphiold)
print("P = P - dphi*(P - Pold)/(dphi - dphiold)")
print("P_new - Pold = {0}".format(P_new - Pold))
print("(dphi - dphiold)= {0}".format(dphi - dphiold))
print("dphi*(P_new - Pold)/(dphi - dphiold)= {0}".format(dphi*(P_new - Pold)/(dphi - dphiold)))
# end if
# c n=n+1
# GO TO 30
# END IF
# P_k_new = P_new - dphi*(P_new - Pold)/(dphi - dphiold)
print("P = {0}".format(P_k_new))
# IF (ABS(dphi).gt.ERRMAX) THEN
ERRMAX = 1e-4
if abs(dphi) > ERRMAX:
# Pold = Plast
# Plast = P
Pold = P_new
P_new = P_k_new
print("Pold = {0}".format(Pold))
print("P_new = {0}".format(P_new))
else:
break
if ITER_p >= 150:
break
ITER_p += 1
PV = P_k_new
print("Presión saturación = {0} {1}".format(PV, "Bar"))
# c WRITE (31,*) ' n=',n
return
# END
# Tr = 0.7
PVini = 4.4
Pvdat = Pc * 10 **-(1.0 + OM)
Trho = Tr
Pvdat = Pc * 10 **-((1./Trho-1.0)*7*(1.0+OM)/3)
a=ac*(3/(2+Trho))**rk
#rk = (A1*Zc+A0)*OM**2+(B1*Zc+B0)*OM+(C1*Zc+C0) # ! initial guess for k parameter
print("kr = {0}".format(rk))
#Tr = 0.7D0 ! Change here to use another Pv than the one at Tr 0.7
#Pvdat=Pc*10**-(1.0D0+OM)
a = ac * (3/(2+Tr))**rk
print("PVini = {0}".format(PVini))
print("Pvdat = {0}".format(Pvdat))
print("dc = {0}".format(dc))
#PVini = Pvdat
vapor_pressure_cal(Tr,PVini,Pvdat)
In [10]:
# subroutine vdWg_Derivs(NDER,T,V,F,F_V,F_2V,F_N)
def vdWg_Derivs_cal(NDER,T,V,d):
"""THE SUBROUTINE CALCULATES THE CONTRIBUTION TO THE RESIDUAL,
REDUCED HELMHOLZ ENERGY (F) AND
ITS FIRST AND SECOND DERIVATIVE WRT V
INPUT:
NDER: indicates which derivatives are required.
1 is used for density calculation and 2 for fugacity
NDER = 1: CALCULATES F, F_V AND F_2V
NDER = 2: CALCULATES F AND F_N
T: TEMPERATURE
V: VOLUME (ML/MOL) or (ML) for checking n-derivatives
OUTPUT: NDER
F: 5 A^RES/RT CONTRIBUTION (DIMENSIONLESS) or (MOLES)
F_V: 5 1ST V-DERIVATIVE OF F
F_2V: 1ST V-DERIVATIVE OF F_V (*V**2)
F_N: 1ST N-DERIVATIVE OF F
"""
# IMPLICIT DOUBLE PRECISION (A-H,O-Z)
# PARAMETER (RGAS=0.08314472d0)
# COMMON /ABd1/ a,b,d
C = (1 - d)/(1 + d)
aRT = a / (RGAS*T)
ETA = 0.25 * b / V
# SUMC = c * b + V
SUMC = C * b + V
SUMD = d * b + V
REP = -np.log(1 - 4 * ETA)
# print(1 - 4 * ETA)
#asc_1 = np.log(SUMD / SUMC)
#ATT = aRT * np.log(SUMD / SUMC) / (b * (C - D))
ATT = aRT
ATTV = aRT / SUMC / SUMD
REPV = 1/(1 - 4 * ETA) - 1
REP2V = 1 / (1 - 4 * ETA)**2 - 1
# ATT2V = aRT * V**2 * (1 / SUMD **2 - 1 / SUMC**2) / (b * (C - D))
ATT2V = aRT * V**2 * (1 / SUMD **2 - 1 / SUMC**2) / (b * (C - d))
F = REP + ATT
F_V = (-REPV /V + ATTV)
# IF (NDER.EQ.1) THEN
if NDER == 1:
F_2V = REP2V - ATT2V
calculo_1 = "F_V"
calculo_2 = "F_2V"
#print(F_2V, calculo_2)
return F, F_V, F_2V, calculo_1, calculo_2
# ELSE
else:
F_N = REP + ATT - V*F_V
calculo = "F_N"
#print(F_N, calculo)
return F_N, calculo
# END IF
In [4]:
#ac = 2422.44
#b = 3.90091
#Tr = 0.7
#NDER, T, V, d = 1, 340.15, 4.003345, 0.24
#ac,b,d,rk = 2422.44, 3.90091, 1.80323, 8.60102
#a = ac * (3 / (2 + Tr))**rk
# METHANE (1)
# 190.564 45.99 0.01155 0.115165 1.00000173664 tc, pc, ohm, vc, zrat
# 2.3277 0.029962 0.932475 1.49541 ac, b, delta1, k
Tc, pc, ohm, dc, zrat = 190.564, 45.99, 0.01155, 0.115165, 1.00000173664
ac, b, delta1, k = 2.3277, 0.029962, 0.932475, 1.49541
Tr = T / Tc
d = delta1
a = ac * (3 / (2 + Tr))**rk
V = 4.0
vdWg = vdWg_Derivs_cal(1,T,V,d)
F_cal = vdWg[0]
F_V_cal = vdWg[1]
F_V2_cal = vdWg[2]
print(F_cal, F_V_cal, F_V2_cal)
In [14]:
# SUBROUTINE VCALC(ITYP,T,P,V)
def VCALC(ITYP,T,P,V):
"""
C
C ROUTINE FOR CALCULATION OF VOLUME, GIVEN PRESSURE
C
C INPUT:
C
C ITYP: TYPE OF ROOT DESIRED
C T: TEMPERATURE
C P: PRESSURE
C
C OUTPUT:
C
C V: VOLUME
C
"""
# IMPLICIT DOUBLE PRECISION (A-H,O-Z)
# PARAMETER (RGAS=0.08314472d0)
# LOGICAL FIRST_RUN
# COMMON /ABd1/ a,b,d1
FIRST_RUN = True
VCP = b * 1.0
S3R = 1.0/VCP
ITER = 0
ZETMIN = 0.00
ZETMAX = 0.99
# IF (ITYP .GT. 0) THEN
if ITYP > 0.0:
ZETA = 0.5
# print("ZETA_init_ITYP_1 = {0}".format(ZETA))
# ELSE
else:
# C..............IDEAL GAS ESTIMATE
ZETA = min(0.5,(VCP * P)/(RGAS * T))
#ZETA = 0.93 * VCP * P/(RGAS * T) # No funciona
#ZETA = 0.5
#print("valores_init_ = {0}".format(VCP * P/(RGAS * T)))
#print("ZETA_init_ITYP_-1 = {0}".format(ZETA))
# END IF
# 100 CONTINUE
while 1:
# C WRITE(*,*)'ZETA',ZETA
V = VCP/ZETA
#print("ITER = {0}".format(ITER))
#print("ZETA = {0}".format(ZETA))
#print("V = {0}".format(V))
# CALL vdWg_Derivs(1,T,V,F,F_V,F_2V,F_N)
vdWg = vdWg_Derivs_cal(1,T,V,d)
F_cal = vdWg[0]
F_V_cal = vdWg[1]
F_V2_cal = vdWg[2]
F = F_cal
F_V = F_V_cal
F_2V = F_V2_cal
PCALC = RGAS * T * (1/V - F_V)
#print("PCALC = {0}".format(PCALC))
# IF (PCALC .GT. P) THEN
if PCALC > P:
ZETMAX = ZETA
# ELSE
else:
ZETMIN = ZETA
# END IF
# c write(*,*)'VCALC V=',V
AT = F - np.log(V) + V*P/(T*RGAS)
DER = RGAS*T*(F_2V+1.0)*S3R
DEL = -(PCALC-P)/DER
#print("DEL = {0}".format(DEL))
#print("V = {0}".format(V))
# ZETA = ZETA + max(min(DEL,0.1),-1.0)
ZETA = ZETA + 1.0 * max(min(DEL,0.1),-0.1)
# IF (ZETA .GT. ZETMAX .OR. ZETA .LT. ZETMIN)
if ZETA > ZETMAX or ZETA < ZETMIN:
#& ZETA = .5D0*(ZETMAX+ZETMIN)
ZETA = 0.5*(ZETMAX + ZETMIN)
# IF (ABS(DEL) .GT. 1D-10) GOTO 100
if abs(DEL) < 1e-10: #1e-10:
break
if ITER >= 150:
break
ITER += 1
"""
# IF (ITYP .EQ. 0 ) THEN
if ITYP == 0:
# C FIRST RUN WAS VAPOUR; RERUN FOR LIQUID
# IF (FIRST_RUN) THEN
if FIRST_RUN:
VVAP = V
AVAP = AT
FIRST_RUN = False
ZETA = 0.5
# GOTO 100
# ELSE
else:
# IF (AT .GT. AVAP) V = VVAP
if AT > AVAP:
V = VVAP
# ENDIF
# ENDIF
"""
return V
# END
"""
# ITYP = 1
T = 148.0 # 110.0
P = 10.0 # 1.00
# METHANE (1)
# 190.564 45.99 0.01155 0.115165 1.00000173664 tc, pc, ohm, vc, zrat
# 2.3277 0.029962 0.932475 1.49541 ac, b, delta1, k
Tc, Pc, ohm, vc, zrat = 190.564, 45.99, 0.01155, 0.115165, 1.00000173664
ac, b, d, rk = 2.3277, 0.029962, 0.932475, 1.49541
Tr = T / Tc
#d = delta1
a = ac * (3 / (2 + Tr))**rk
dc = 1/vc
OM = ohm
"""
Vo = 1.0 * b
#print(Vo)
V_cal = VCALC(1,T,P,Vo)
print("V_cal = {0}".format(V_cal))
In [15]:
V_cal_L = VCALC(1,T,P,Vo)
print("V_cal = {0}".format(V_cal_L))
V_cal_V = VCALC(-1,T,P,Vo)
print("V_cal = {0}".format(V_cal_V))
In [7]:
print("V_cal_L = {0}".format(V_cal_L))
print("V_cal_V = {0}".format(V_cal_V))
#V_cal_L = 0.04985655425118943
#V_cal_V = 3.463117452110996
In [ ]:
# METHANE (1)
# 190.564 45.99 0.01155 0.115165 1.00000173664 tc, pc, ohm, vc, zrat
# 2.3277 0.029962 0.932475 1.49541 ac, b, delta1, k
Tc, pc, ohm, dc, zrat = 190.564, 45.99, 0.01155, 0.115165, 1.00000173664
ac, b, delta1, k = 2.3277, 0.029962, 0.932475, 1.49541
d = delta1
In [ ]:
# V = VCP/ZETA
ZATA_p = b/0.04985655425118943
ZATA_p
In [8]:
# SUBROUTINE FUG_CALC(T,P,V,phi)
def FUG_CALC(T,P,V):
# IMPLICIT DOUBLE PRECISION (A-H,O-Z)
# PARAMETER (RGAS=0.08314472d0)
RT = RGAS*T
Z = P*V/RT
# CALL vdWg_Derivs(2,T,V,F,F_V,F_2V,F_N)
vdWg = vdWg_Derivs_cal(2,T,V,d)
F_N = vdWg[0]
PHI = np.exp(F_N)/Z
# print("PHI = {0}".format(PHI))
return PHI
# END
#FUG_CALC(T,P,V)
#V_cal_L = 0.04985655425118943
#V_cal_V = 3.463117452110996
FUG_CALC(T,P,V_cal_L)
FUG_CALC(T,P,V_cal_V)
Out[8]:
In [8]:
Tr
Out[8]: