혼합물의 기체-액체 상평형 (vapor-liquid phase equilibrium)¶
The notebook executes VLE_mixture.py that uses DATA.txt as input¶
In [1]:
%ls DATA.txt *.py
C 드라이브의 볼륨: Windows10
볼륨 일련 번호: 70CE-288B
C:\DATA\Ch\Chang\Theory\SAFT PC-SAFT\Mixture\2020 PC-SAFT 디렉터리
C:\DATA\Ch\Chang\Theory\SAFT PC-SAFT\Mixture\2020 PC-SAFT 디렉터리
2021-02-15 오후 01:47 1,068 DATA.txt
2018-03-07 오전 04:43 18,456 pcsaft_mix.py
2018-04-19 오후 04:31 42,322 VLE_mixture.py
3개 파일 61,846 바이트
0개 디렉터리 80,717,217,792 바이트 남음
In [2]:
!type DATA.txt
[ ethanol and water binary mixture ] # EOS : 0 for PR 1 for PC-SAFT 1 # ncomp 2 ----------------------------------------------- ethanol (polar) # Mw (g/mol) Tc (K) Pc (bar) omega Tnb (K) 46.069 513.9 61.4 0.644 351.4 # m sig (A) epsk (K) kAB eAB (K) dipole(D) xp 2.2049 3.2774 187.24 0.03363 2652.7 1.7 0.29466 ----------------------------------------------- water (polar) # Mw (g/mol) Tc (K) Pc (bar) omega Tnb (K) 18.015 647.3 221.2 0.344 373.15 # m sig (A) epsk (K) kAB eAB (K) dipole(D) xp 1.0405 2.9657 175.15 0.08924 2706.7 1.85 0.66245 ----------------------------------------------- # k12 k13 .. k1N, k23 ... 0. # i j kAiBj eAiBj (K) -- cross association override ------------------------------------ end of cross association # type of calculation: 0(Z), 1(phi), 2(BublP) 3(DewP) 4(BublT) 5(DewT) 6(Flash) 4 # T (K) P (bar) z1 z2 ... zN 351.4 1.013 -0.5 0.5 # scaling fraction for initial liquid volume 0.99 --EOF------------------------------------------
In [3]:
%matplotlib inline
In [4]:
%run VLE_mixture.py
=========================================
Vapor-Liquid Equilibrium of Mixture
=========================================
[ ethanol and water binary mixture ]
using PC-SAFT eos ...
Component 1 = ethanol (polar)
Mw (g/mol) Tc (K) Pc (bar) omega Tb (K)
46.069 513.90 61.400 0.644 351.40
m sig (A) epsk (K) kAB eAB (K) dipole (D) xp
2.2049 3.2774 187.24 0.033630 2652.70 1.700 0.29466
Component 2 = water (polar)
Mw (g/mol) Tc (K) Pc (bar) omega Tb (K)
18.015 647.30 221.200 0.344 373.15
m sig (A) epsk (K) kAB eAB (K) dipole (D) xp
1.0405 2.9657 175.15 0.089240 2706.70 1.850 0.66245
Kij's are
0.0000 0.0000
0.0000 0.0000
cross association override in PC-SAFT
i j kAiBj eAiBj (K)
The number of overrides is = 0
checking associating components in PC-SAFT ...
component 1 is associating
component 2 is associating
The number of associating components = 2
Bubble T calculation at P = 1.013 (bar)
scaling fraction for initial liquid volume = 0.99
-- Experimental data --
T (K) P (bar) x1 y1
press Enter to proceed
<< Txy - Bubble point temperature calculation for given P, x >>
P = 1.013 (bar)
T (K) x1 y1 iteration
-----------------------------------------
374.55 0.0000 0.0000 4
372.88 0.0050 0.0586 4
371.37 0.0100 0.1089 4
370.05 0.0150 0.1520 4
368.84 0.0200 0.1896 5
367.76 0.0250 0.2224 5
366.80 0.0300 0.2511 4
365.91 0.0350 0.2766 4
365.11 0.0400 0.2992 4
364.38 0.0450 0.3195 4
363.71 0.0500 0.3377 4
363.10 0.0550 0.3540 4
362.55 0.0600 0.3688 4
362.02 0.0650 0.3824 4
361.56 0.0700 0.3945 4
361.11 0.0750 0.4058 4
360.71 0.0800 0.4161 4
360.36 0.0850 0.4253 4
360.00 0.0900 0.4341 4
359.70 0.0950 0.4420 4
359.40 0.1000 0.4494 4
359.11 0.1050 0.4564 4
358.87 0.1100 0.4627 4
358.62 0.1150 0.4688 4
358.39 0.1200 0.4744 4
358.20 0.1250 0.4795 4
357.99 0.1300 0.4845 4
357.81 0.1350 0.4891 4
357.65 0.1400 0.4933 4
357.48 0.1450 0.4974 4
357.32 0.1500 0.5014 4
357.17 0.1550 0.5051 4
357.05 0.1600 0.5085 4
356.92 0.1650 0.5118 4
356.78 0.1700 0.5151 4
356.67 0.1750 0.5181 4
356.56 0.1800 0.5209 4
356.46 0.1850 0.5236 3
356.36 0.1900 0.5263 4
356.25 0.1950 0.5290 4
356.16 0.2000 0.5314 3
356.08 0.2050 0.5337 3
355.99 0.2100 0.5361 3
355.91 0.2150 0.5384 3
355.83 0.2200 0.5405 3
355.76 0.2250 0.5426 3
355.68 0.2300 0.5447 3
355.61 0.2350 0.5467 3
355.55 0.2400 0.5486 3
355.48 0.2450 0.5506 3
355.41 0.2500 0.5525 3
355.35 0.2550 0.5543 3
355.29 0.2600 0.5561 3
355.23 0.2650 0.5579 3
355.17 0.2700 0.5597 3
355.11 0.2750 0.5615 3
355.05 0.2800 0.5632 3
355.00 0.2850 0.5649 3
354.94 0.2900 0.5667 3
354.89 0.2950 0.5684 3
354.84 0.3000 0.5700 3
354.78 0.3050 0.5717 3
354.73 0.3100 0.5734 3
354.68 0.3150 0.5751 3
354.63 0.3200 0.5767 3
354.58 0.3250 0.5784 3
354.53 0.3300 0.5801 3
354.49 0.3350 0.5818 3
354.44 0.3400 0.5834 3
354.39 0.3450 0.5851 3
354.35 0.3500 0.5868 3
354.30 0.3550 0.5885 3
354.25 0.3600 0.5901 3
354.21 0.3650 0.5918 3
354.16 0.3700 0.5936 3
354.12 0.3750 0.5953 3
354.08 0.3800 0.5970 3
354.04 0.3850 0.5986 3
354.00 0.3900 0.6004 3
353.96 0.3950 0.6021 3
353.91 0.4000 0.6039 3
353.87 0.4050 0.6057 3
353.83 0.4100 0.6075 3
353.79 0.4150 0.6093 3
353.74 0.4200 0.6111 3
353.70 0.4250 0.6130 3
353.66 0.4300 0.6148 3
353.62 0.4350 0.6167 3
353.58 0.4400 0.6186 3
353.54 0.4450 0.6205 3
353.50 0.4500 0.6224 3
353.47 0.4550 0.6243 3
353.43 0.4600 0.6263 3
353.39 0.4650 0.6282 3
353.35 0.4700 0.6302 3
353.31 0.4750 0.6322 3
353.27 0.4800 0.6343 3
353.24 0.4850 0.6363 3
353.20 0.4900 0.6384 3
353.16 0.4950 0.6405 3
353.12 0.5000 0.6426 3
353.09 0.5050 0.6447 3
353.05 0.5100 0.6469 3
353.02 0.5150 0.6490 3
352.98 0.5200 0.6512 3
352.95 0.5250 0.6534 3
352.91 0.5300 0.6557 3
352.88 0.5350 0.6579 3
352.84 0.5400 0.6602 3
352.81 0.5450 0.6625 3
352.77 0.5500 0.6649 3
352.74 0.5550 0.6672 3
352.71 0.5600 0.6696 3
352.68 0.5650 0.6720 3
352.64 0.5700 0.6744 3
352.61 0.5750 0.6769 3
352.58 0.5800 0.6794 3
352.55 0.5850 0.6819 3
352.52 0.5900 0.6844 3
352.49 0.5950 0.6869 3
352.46 0.6000 0.6895 3
352.43 0.6050 0.6921 3
352.40 0.6100 0.6947 3
352.37 0.6150 0.6974 3
352.34 0.6200 0.7001 3
352.32 0.6250 0.7027 3
352.29 0.6300 0.7055 3
352.26 0.6350 0.7082 3
352.23 0.6400 0.7110 3
352.21 0.6450 0.7138 3
352.18 0.6500 0.7167 3
352.16 0.6550 0.7196 3
352.13 0.6600 0.7225 3
352.11 0.6650 0.7254 3
352.08 0.6700 0.7284 3
352.05 0.6750 0.7314 3
352.03 0.6800 0.7344 3
352.01 0.6850 0.7374 3
351.99 0.6900 0.7405 3
351.96 0.6950 0.7436 3
351.94 0.7000 0.7467 2
351.92 0.7050 0.7499 2
351.90 0.7100 0.7531 2
351.88 0.7150 0.7563 2
351.86 0.7200 0.7596 2
351.84 0.7250 0.7628 2
351.82 0.7300 0.7662 2
351.81 0.7350 0.7695 2
351.79 0.7400 0.7729 2
351.77 0.7450 0.7763 2
351.76 0.7500 0.7797 2
351.74 0.7550 0.7832 2
351.72 0.7600 0.7868 2
351.71 0.7650 0.7903 2
351.69 0.7700 0.7939 2
351.68 0.7750 0.7975 2
351.66 0.7800 0.8011 2
351.65 0.7850 0.8048 2
351.64 0.7900 0.8085 2
351.63 0.7950 0.8123 2
351.61 0.8000 0.8161 2
351.61 0.8050 0.8199 2
351.60 0.8100 0.8237 2
351.58 0.8150 0.8277 2
351.57 0.8200 0.8316 2
351.57 0.8250 0.8355 2
351.56 0.8300 0.8395 2
351.55 0.8350 0.8436 2
351.54 0.8400 0.8477 2
351.54 0.8450 0.8518 2
351.53 0.8500 0.8559 2
351.53 0.8550 0.8601 2
351.52 0.8600 0.8644 2
351.52 0.8650 0.8686 2
351.52 0.8700 0.8729 2
351.51 0.8750 0.8773 2
351.51 0.8800 0.8817 2
351.51 0.8850 0.8861 2
351.51 0.8900 0.8906 2
351.50 0.8950 0.8951 2
351.51 0.9000 0.8996 2
351.51 0.9050 0.9042 2
351.51 0.9100 0.9089 2
351.51 0.9150 0.9136 2
351.52 0.9200 0.9183 2
351.52 0.9250 0.9230 2
351.52 0.9300 0.9278 2
351.52 0.9350 0.9327 2
351.53 0.9400 0.9376 2
351.53 0.9450 0.9425 1
351.54 0.9500 0.9475 2
351.55 0.9550 0.9526 2
351.55 0.9600 0.9577 2
351.56 0.9650 0.9628 2
351.57 0.9700 0.9680 2
351.58 0.9750 0.9732 2
351.60 0.9800 0.9784 2
351.61 0.9850 0.9838 2
351.61 0.9900 0.9891 2
351.63 0.9950 0.9945 2
351.64 1.0000 1.0000 2
Press Enter to quit
VLE_mixture.py¶
In [5]:
!type VLE_mixture.py
# VLE_mixture.py | 2016-9 ~ 2018-3 J. Chang
# Vapor-Liquid Equilibrium of Mixture
from math import *
import matplotlib.pyplot as plt
from pcsaft_mix import *
def Zpr( V, T, z ) : # returns compressibility factor of Peng-Robinson equation of state
global R, Tc, Pc, omega, ncP, Keps
# V molar volume (m3/mol)
# unit of pressure (bar)
a = [0]*ncP
b = [0]*ncP
for i in range(1,ncP) :
kapa = 0.37464+1.54226*omega[i]-0.26992*omega[i]**2
alpha = (1+kapa*( 1- sqrt(T/Tc[i]) ))**2
a[i] = 0.45724 * R**2 * Tc[i]**2 / Pc[i] * alpha
b[i] = 0.07780 * R * Tc[i] / Pc[i]
amix = 0
bmix = 0
for i in range(1,ncP) :
bmix += z[i] * b[i]
for j in range(1,ncP) :
if i == j : aij = a[i]
else : aij = sqrt( a[i]*a[j]) * (1 - Keps[i][j])
amix += z[i]*z[j]*aij
Z = V/( V-bmix ) - amix/(R*T) * V / ( V*(V+bmix) + bmix*(V-bmix) )
return Z
def phi_pr( V, T, z ) : # returns fugacity coefficients
global R, Tc, Pc, omega, ncP, Keps
a = [0]*ncP
b = [0]*ncP
for i in range(1,ncP) :
kapa = 0.37464+1.54226*omega[i]-0.26992*omega[i]**2
alpha = (1+kapa*( 1- sqrt(T/Tc[i]) ))**2
a[i] = 0.45724 * R**2 * Tc[i]**2 / Pc[i] * alpha
b[i] = 0.07780 * R * Tc[i] / Pc[i]
amix = 0
bmix = 0
za = [0]*ncP
for i in range(1,ncP) :
bmix += z[i] * b[i]
for j in range(1,ncP) :
if i == j : aij = a[i]
else : aij = sqrt( a[i]*a[j]) * (1 - Keps[i][j])
amix += z[i]*z[j]*aij
za[i] += z[j]*aij
Z = V/( V-bmix ) - amix/(R*T) * V / ( V*(V+bmix) + bmix*(V-bmix) )
root2 = sqrt(2)
phi = [0]*ncP
for i in range(1,ncP) :
lnphi = b[i]/bmix *(Z - 1) - log( Z - bmix*Z/V ) - amix / (2*root2*bmix*R*T) \
* ( 2*za[i]/amix - b[i]/bmix ) * log( (1+(1+root2)*bmix/V) / (1+(1-root2)*bmix/V) )
phi[i] = exp( lnphi )
return phi
def PR_ab( T, z ) : # returns coefficients a and b of PR equation of state
# used to find analytic solution of PR eos
global R, Tc, Pc, omega, ncP, Keps
a = [0]*ncP
b = [0]*ncP
for i in range(1,ncP) :
kapa = 0.37464+1.54226*omega[i]-0.26992*omega[i]**2
alpha = (1+kapa*( 1- sqrt(T/Tc[i]) ))**2
a[i] = 0.45724 * R**2 * Tc[i]**2 / Pc[i] * alpha
b[i] = 0.07780 * R * Tc[i] / Pc[i]
amix = 0
bmix = 0
for i in range(1,ncP) :
bmix += z[i] * b[i]
for j in range(1,ncP) :
if i == j : aij = a[i]
else : aij = sqrt( a[i]*a[j]) * (1 - Keps[i][j])
amix += z[i]*z[j]*aij
return amix, bmix
def Fugacity( V, T, z, P ) : # calls fugacity routine of equation of state, and return fugacity coefficients
global eos, ncP
global m, sig, epsk, kAB, eAB, nc, Keps, ncas, dipol, xp
# fugacity coefficient
if eos == 'PR' :
phi = phi_pr( V, T, z )
elif eos == 'PC-SAFT' :
phi = pcsaft( 'Phi', V, T, z, m, sig, epsk, kAB, eAB, nc, Keps, ncas, dipol, xp )
# fugacity
fug = [0]*ncP
for i in range(1,ncP) : fug[i] = phi[i] * z[i] * P
return phi, fug
def Peos( V, T, z ) : # calls equation of state, and returns Peos(V)
global eos, R
global m, sig, epsk, kAB, eAB, nc, Keps, ncas, dipol, xp
if eos == 'PR' : Z = Zpr( V, T, z )
elif eos == 'PC-SAFT' : Z = pcsaft( 'Z', V, T, z, m, sig, epsk, kAB, eAB, nc, Keps, ncas, dipol, xp )
Peos = (R*T/V) * Z
return Peos, Z
def Peos_P( V, P, T, z ) : # calls equation of state, and returns Peos(V) - P
global eos, R
global m, sig, epsk, kAB, eAB, nc, Keps, ncas, dipol, xp
if eos == 'PR' : Z = Zpr( V, T, z )
elif eos == 'PC-SAFT' : Z = pcsaft( 'Z', V, T, z, m, sig, epsk, kAB, eAB, nc, Keps, ncas, dipol, xp )
Peos = (R*T/V) * Z
return Peos - P
def newton__Peos_P( V0, P, T, z ) : # returns molar volume that satisfies Peos(V) - P = 0 by Newton method
# V0 is initial guess of molar volume
# returns None when convergence error occurs
MAXIT = 100 # maximum 100 iterations
relax = 1. # 0.3
#VV = [ V0 ] # Keep the record
V = V0
dV = V * 1e-5
for i in range(MAXIT) :
fx = Peos_P( V, P, T, z )
if abs( fx/P ) < 1e-5 : break # converged
dfdV = ( Peos_P( V+dV, P, T, z ) - fx ) / dV
if abs( dfdV * V/P ) < 1e-5 : # convergence error
#print( VV ) # check how bad convergence records are
return None
delV = - fx / dfdV
#if abs(delV/V) > 0.1 : relax *= 0.5 # when large change occurs (>10%), reduce relax
V += delV * relax
#VV.append( V )
return V
def Find_vol_PR( P, T, z, phase ) : # returns molar volumes of PR equation of state using analytical formula
global R
pi = 3.1415926536
a, b = PR_ab( T, z )
# using math formula for cubic equation x^3 + a2 * x^2 + a1 * x + a0 = 0
a2 = b-R*T/P
a1 = -3*b*b + (a-2*b*R*T)/P
a0 = b*b*b + (b*b*R*T-a*b)/P
q = ( 3*a1 - a2*a2 ) / 9
r = ( 9*a1*a2 - 27*a0 - 2*a2*a2*a2 ) / 54
D = q*q*q + r*r
if D < 0 : # three real roots, q < 0
# The following formula are for D < 0, there exist three real roots
th = acos( r / sqrt(-q*q*q) )
V1 = 2*sqrt(-q) * cos( th/3 ) - a2 / 3
V2 = 2*sqrt(-q) * cos( th/3 + 2/3*pi ) - a2 / 3
V3 = 2*sqrt(-q) * cos( th/3 + 4/3*pi ) - a2 / 3
# vapor volume is the largest, and liquid volume is the smallest one
Vvap = V1
if V2 > Vvap : Vvap = V2
if V3 > Vvap : Vvap = V3
Vliq = V1
if V2 < Vliq : Vliq = V2
if V3 < Vliq : Vliq = V3
else : # D > 0, single real root
s = r + sqrt(D)
if s > 0 : s = s**(1/3) # s arg > 0
else : s = -(-s)**(1/3) # s arg < 0
t = r - sqrt(D)
if t > 0 : t = t**(1/3) # t arg > 0
else : t = -(-t)**(1/3) # t arg < 0
Vvap = -a2 / 3 + s + t # one real root, two complex roots not physical
Vliq = Vvap
if phase == 'vapor' :
Vol = Vvap
else :
Vol = Vliq
return Vol
def Find_vol( P, T, z, phase ) : # returns molar volumes of general equations of state
global eos, R, Tc, Pc, ncP, scalVl0, SavedResult
if phase == 'vapor' and SavedResult[ 'existVvap' ] == False : # for a single calculation or the first calculation in scanning analysis
# analytic solution of PR eos is used as initial guess
Vol = Find_vol_PR( P, T, z, phase )
if eos == 'PR' : return Vol # solution already at hand # return and exit
# find numerical solution for general equation of state
Vol = newton__Peos_P( Vol, P, T, z )
# SavedResult[ 'existVvap' ] = True
# SavedResult[ 'Vvap' ] = Vol
elif phase == 'vapor' and SavedResult[ 'existVvap' ] == True :
# when previously saved results are available, use it as initial guess
Vol = SavedResult[ 'Vvap' ]
Vol = newton__Peos_P( Vol, P, T, z )
# SavedResult[ 'Vvap' ] = Vol
elif phase == 'liquid' and SavedResult[ 'existVliq' ] == False : # for a single calculation or the first calculation in scanning analysis
# analytic solution of PR eos is used as initial guess
Vol = Find_vol_PR( P, T, z, phase )
if eos == 'PR' : return Vol # solution already at hand # return and exit
# find numerical solution for general equation of state
Vol = newton__Peos_P( Vol*scalVl0, P, T, z ) # use different guess value of liquid volume than PR eos using scaling fraction
# SavedResult[ 'existVliq' ] = True
# SavedResult[ 'Vliq' ] = Vol
elif phase == 'liquid' and SavedResult[ 'existVliq' ] == True :
# when previously saved results are available, use it as initial guess
Vol = SavedResult[ 'Vliq' ]
Vol = newton__Peos_P( Vol, P, T, z )
# SavedResult[ 'Vliq' ] = Vol
return Vol
def BubbleP_ini( T, x ) : # returns initial guesses of P and y
global Tc, Pc, Tb, ncP
Ps = [0]*ncP
# Guess bubble pressure from vapor pressure and x's
# Estimated vapor pressure from Tc, Pc and Tb by using ln(Ps) = A - B / T
# print()
# print("Estimated (hypothetical) vapor pressures from Tc, Pc, Tb " )
P = 0.
for i in range(1,ncP) :
B = ( log(Pc[i]) - log(1.013) ) / ( 1/Tb[i] - 1/Tc[i] )
A = log(1.013) + B / Tb[i]
Ps[i] = exp( A - B / T )
# print(" Component %d, Psat = %9.3f (bar) " %( i, Ps[i] ) )
P += x[i] * Ps[i] # assuming ideal solution
y = [0]*ncP
for i in range(1,ncP) : y[i] = x[i] * Ps[i] / P # assuming ideal solution
# initial guesses
# print(" P = %9.3f (bar) <-- guess value from estimated vapor pressures" % P )
# print(" y = ", end=" " )
# for i in range(1,ncP) :
# print("%8.5f" % y[i], end=" ")
# print(" <-- initial guess \n")
return P, y
def BubbleP( T, x, P, y ) : # P, y are initial guesses
global ncP
Maxit = 3000
for it in range(Maxit) :
# Get fugacity of liquid and volume
Vliq = Find_vol( P, T, x, 'liquid' )
phiL, fL = Fugacity( Vliq, T, x, P )
# Get fugacity of vapor and volume
Vvap = Find_vol( P, T, y, 'vapor' )
phiV, fV = Fugacity( Vvap, T, y, P )
# Adjust y composition
ynew = [0]*ncP
sumy = 0
for i in range(1,ncP) :
ynew[i] = x[i] * ( phiL[i] / phiV[i] )
# ynew[i] = y[i] * ( fL[i] / fV[i] ) # Division by zero for pure substance
sumy += ynew[i]
# Adjust P
Pnew = P * sumy
# Convergence check for fugacity equality and constancy of pressure
Converged = True
for i in range(1,ncP) :
if abs( fL[i] - fV[i] ) > 1e-5 : Converged = False
if abs( Pnew - P ) < 1e-5 and Converged : break
# Update y and P
P = Pnew
y = ynew
# Error check for single phase
if abs(Vvap-Vliq)/Vvap < 1e-5 :
it = None
return P, y, Vliq, Vvap, it
# Print running results
# print(" iter = %2d P = %8.3f " %( it, P ), end=" " )
# print(" y = ", end=" " )
# for i in range(1,ncP) : print("%8.5f" %y[i], end=" ")
# print()
return P, y, Vliq, Vvap, it
def DewP_ini( T, y ) : # returns initial guesses of P and x
global Tc, Pc, Tb, ncP
# Guess bubble pressure from vapor pressure and x's
# Estimated vapor pressure from Tc, Pc and Tb by using ln(Ps) = A - B / T
# print("Estimated (hypothetical) vapor pressures from Tc, Pc, Tb " )
Ps = [0]*ncP
sum = 0.
for i in range(1,ncP) :
B = ( log(Pc[i]) - log(1.013) ) / ( 1/Tb[i] - 1/Tc[i] )
A = log(1.013) + B / Tb[i]
Ps[i] = exp( A - B / T )
# print(" Component %d, Psat = %9.3f (bar) " %( i, Ps[i] ) )
sum += y[i] / Ps[i] # assuming ideal solution
P = 1 / sum
x = [0]*ncP
for i in range(1,ncP) : x[i] = y[i] * P / Ps[i] # assuming ideal solution
return P, x
def DewP( T, y, P, x ) : # P, x are initial guesses
global ncP
Maxit = 3000
for it in range(Maxit) :
# Get fugacity of vapor
Vvap = Find_vol( P, T, y, 'vapor' )
phiV, fV = Fugacity( Vvap, T, y, P )
# Get fugacity of liquid
Vliq = Find_vol( P, T, x, 'liquid' )
phiL, fL = Fugacity( Vliq, T, x, P )
# Adjust x composition
xnew = [0]*ncP
sumx = 0
for i in range(1,ncP) :
xnew[i] = y[i] * ( phiV[i] / phiL[i] )
sumx += xnew[i]
# Adjust P
Pnew = P / sumx
# Convergence check for fugacity equality and constancy of pressure
Converged = True
for i in range(1,ncP) :
if abs( fL[i] - fV[i] ) > 1e-5 : Converged = False
if abs( Pnew - P ) < 1e-5 and Converged : break
# Update x and P
P = Pnew
x = xnew
# Error check for single phase
if abs(Vvap-Vliq)/Vvap < 1e-5 :
it = None
return P, x, Vliq, Vvap, it
return P, x, Vliq, Vvap, it
def BubbleT_ini( P, x ) : # returns initial guesses of T and y
global Tc, Pc, Tb, ncP, Bp
Ap = [0]*ncP
Bp = [0]*ncP
Ps = [0]*ncP
# Guess bubble temperature from vapor pressure and x's
# Estimated boiling temperature from Tc, Pc and Tb by using ln(Ps) = A - B / T
# print("\n")
# print("Estimated boiling temperatures from Tc, Pc, Tb " )
T = 0.
for i in range(1,ncP) :
B = ( log(Pc[i]) - log(1.013) ) / ( 1/Tb[i] - 1/Tc[i] )
A = log(1.013) + B / Tb[i]
Tbp = B / ( A - log(P) )
Ap[i] = A
Bp[i] = B
# print(" Component %d, Tb = %9.2f (K) " %( i, Tbp ))
T += x[i] * Tbp # mole fraction average
# initial guesses
# print(" T = %9.2f (K) <-- guess value as mole fraction average" % T )
# Estimated vapor pressure from Tc, Pc and Tb by using ln(Ps) = A - B / T
# print("Estimated (hypothetical) vapor pressures at this T " )
y = [0]*ncP
for i in range(1,ncP) :
Ps[i] = exp( Ap[i] - Bp[i] / T )
# print(" Component %d, Psat = %9.3f (bar) " %( i, Ps[i] ))
y[i] = x[i] * Ps[i] / P # assuming ideal solution
# print(" y = ", end=" " )
# for i in range(1,ncP) : print("%8.5f" % y[i], end=" ")
# print(" <-- initial guess \n")
return T, y
def BubbleT( P, x, T, y ) : # T, y are initial guesses
global ncP, Bp
B = 0 # parameter to adjust temperature
for i in range(1,ncP) : B += x[i]*Bp[i] # mole fraction average
Maxit = 3000
for it in range(Maxit) :
# Get fugacity of liquid
Vliq = Find_vol( P, T, x, 'liquid' )
phiL, fL = Fugacity( Vliq, T, x, P )
# Get fugacity of vapor
Vvap = Find_vol( P, T, y, 'vapor' )
phiV, fV = Fugacity( Vvap, T, y, P )
# Adjust y composition
ynew = [0]*ncP
sumy = 0
for i in range(1,ncP) :
ynew[i] = x[i] * ( phiL[i] / phiV[i] )
sumy += ynew[i]
# Adjust T
# ln(Ps) = A - B/T -> d ln(Ps) = B dT/T^2
# Pnew = P * sumy -> - ln sumy = B/T^2 dT
delT = -log(sumy) * T**2 / B
# Convergence check for fugacity equality and constancy of pressure
Converged = True
for i in range(1,ncP) :
if abs( fL[i] - fV[i] ) > 1e-5 : Converged = False
if abs( delT ) < 1e-3 and Converged : break
# Update y and T
T = T + delT
y = ynew
# Error check for single phase
if abs(Vvap-Vliq)/Vvap < 1e-5 :
it = None
return T, y, Vliq, Vvap, it
# Print running results
# print(" iter = %2d T = %8.3f " %( it, T ), end=" " )
# print(" y = ", end=" " )
# for i in range(1,ncP) :
# print("%8.5f" % y[i], end=" ")
# print()
return T, y, Vliq, Vvap, it
def DewT_ini( P, y ) : # returns initial guesses of T and x
global Tc, Pc, Tb, ncP, Bp
Ap = [0]*ncP
Bp = [0]*ncP
Ps = [0]*ncP
# Guess bubble temperature from vapor pressure and x's
# Estimated boiling temperature from Tc, Pc and Tb by using ln(Ps) = A - B / T
# print("\n")
# print("Estimated boiling temperatures from Tc, Pc, Tb " )
T = 0.
for i in range(1,ncP) :
B = ( log(Pc[i]) - log(1.013) ) / ( 1/Tb[i] - 1/Tc[i] )
A = log(1.013) + B / Tb[i]
Tbp = B / ( A - log(P) )
Ap[i] = A
Bp[i] = B
# print(" Component %d, Tb = %9.2f (K) " %( i, Tbp ))
T += y[i] * Tbp # mole fraction average
# initial guesses
# print(" T = %9.2f (K) <-- guess value as mole fraction average" % T )
# Estimated vapor pressure from Tc, Pc and Tb by using ln(Ps) = A - B / T
# print("Estimated (hypothetical) vapor pressures at this T " )
x = [0]*ncP
for i in range(1,ncP) :
Ps[i] = exp( Ap[i] - Bp[i] / T )
# print(" Component %d, Psat = %9.3f (bar) " %( i, Ps[i] ))
x[i] = y[i] * P / Ps[i] # assuming ideal solution
# print(" x = ", end=" " )
# for i in range(1,ncP) : print("%8.5f" % x[i], end=" ")
# print(" <-- initial guess \n")
return T, x
def DewT( P, y, T, x ) : # T, x are initial guesses
global ncP, Bp
B = 0 # parameter to adjust temperature
for i in range(1,ncP) : B += y[i]*Bp[i] # mole fraction average
Maxit = 3000
for it in range(Maxit) :
# Get fugacity of vapor
Vvap = Find_vol( P, T, y, 'vapor' )
phiV, fV = Fugacity( Vvap, T, y, P )
# Get fugacity of liquid
Vliq = Find_vol( P, T, x, 'liquid' )
phiL, fL = Fugacity( Vliq, T, x, P )
# Adjust x composition
xnew = [0]*ncP
sumx = 0
for i in range(1,ncP) :
xnew[i] = y[i] * ( phiV[i] / phiL[i] )
sumx += xnew[i]
# Adjust T
# ln(Ps) = A - B/T -> d ln(Ps) = B dT/T^2
# Pnew = P / sumx -> ln sumx = B/T^2 dT
delT = log(sumx) * T**2 / B
# Convergence check for fugacity equality and constancy of pressure
Converged = True
for i in range(1,ncP) :
if abs( fL[i] - fV[i] ) > 1e-5 : Converged = False
if abs( delT ) < 1e-3 and Converged : break
# Update x and T
T = T + delT
x = xnew
# Error check for single phase
if abs(Vvap-Vliq)/Vvap < 1e-5 :
it = None
return T, x, Vliq, Vvap, it
# Print running results
# print(" iter = %2d T = %8.3f " %( it, T ), end=" " )
# print(" x = ", end=" " )
# for i in range(1,ncP) :
# print("%8.5f" % x[i], end=" ")
# print()
return T, x, Vliq, Vvap, it
def Flash( P, T, z, x, y, v ) : # x, y, v are initial guesses
global ncP
Maxit = 3000
K = [0]*ncP
for it in range(Maxit) :
# Get fugacity of liquid
Vliq = Find_vol( P, T, x, 'liquid' )
phiL, fL = Fugacity( Vliq, T, x, P )
# Get fugacity of vapor
Vvap = Find_vol( P, T, y, 'vapor' )
phiV, fV = Fugacity( Vvap, T, y, P )
# K value
for i in range(1,ncP) : K[i] = phiL[i] / phiV[i]
# Find vapor fraction v
vold = v
for iv in range(100) :
F = 0
dFdv = 0
for i in range(1,ncP) :
F += z[i] * (K[i]-1) / ( 1 + v * (K[i]-1) )
dFdv += - z[i] * (K[i]-1)**2 / ( 1 + v * (K[i]-1) )**2
delv = - F / dFdv
if abs(delv) < 1e-5 : break
v += delv
# Convergence check for fugacity equality and constancy of v
Converged = True
for i in range(1,ncP) :
if abs( fL[i] - fV[i] ) > 1e-5 : Converged = False
if abs( v - vold ) < 1e-5 and Converged : break
# update x, y
for i in range(1,ncP) :
x[i] = z[i] / ( 1 + v * (K[i]-1) )
y[i] = K[i] * x[i]
# Print running results
# print(" iter = %2d v = %8.5f " %( it, v ), end=" " )
# print(" x, y = ", end=" " )
# for i in range(1,ncP) : print("%8.4f %8.4f" %(x[i],y[i]), end=" ")
# print()
return x, y, K, Vliq, Vvap, v, it
def Read_data__initialize() :
global eos, R, component, Mw, Tc, Pc, omega, Tb, nc, ncP, Keps, scalVl0, plot_title, SavedResult
global Texp, Pexp, Xexp, Yexp, isTPxySet
global m, sig, epsk, kAB, eAB, ncas, dipol, xp
R = 8.314e-5 # (m3.bar/mol/K)
# Read data
f = open("Data.txt",'r')
mixture = f.readline()
f.readline() # skip comment line
eos = int( f.readline() )
f.readline()
nc = int( f.readline() ) # number of components
f.readline() # --------------------
print( mixture )
if eos == 0 : eos = 'PR'
elif eos == 1 : eos = 'PC-SAFT'
print("using %s eos ... \n" % eos )
ncP = nc+1 # used for set the size of list
component = [0]*ncP
Mw = [0]*ncP
Tc = [0]*ncP
Pc = [0]*ncP
omega = [0]*ncP
Tb = [0]*ncP
m = [0]*ncP
sig = [0]*ncP
epsk = [0]*ncP
kAB = [ [0]*ncP for i in range(ncP) ]
eAB = [ [0]*ncP for i in range(ncP) ]
dipol= [0]*ncP
xp = [0]*ncP
for i in range(1,ncP) :
component[i] = f.readline().replace("\n","")
f.readline()
Mw[i], Tc[i], Pc[i], omega[i], Tb[i] = map( float, f.readline().split() )
f.readline()
pcsaftDATA = f.readline().split()
if len(pcsaftDATA) == 5 :
m[i], sig[i], epsk[i], kAB[i][i], eAB[i][i] = map( float, pcsaftDATA )
elif len(pcsaftDATA) == 7 : # polar component
m[i], sig[i], epsk[i], kAB[i][i], eAB[i][i], dipol[i], xp[i] = map( float, pcsaftDATA )
f.readline() # -----------------------------------------------------
# Echo Property
print("Component %d = %s " %( i, component[i] ) )
print(" Mw (g/mol) Tc (K) Pc (bar) omega Tb (K) ")
print(" %8.3f %8.2f %8.3f %7.3f %8.2f" %( Mw[i], Tc[i], Pc[i], omega[i], Tb[i] ) )
print(" m sig (A) epsk (K) kAB eAB (K) dipole (D) xp ")
print(" %8.4f %8.4f %8.2f %11.6f %8.2f %8.3f %9.5f \n"
%( m[i], sig[i], epsk[i], kAB[i][i], eAB[i][i], dipol[i], xp[i] ) )
# check diameter in (m)
if sig[i] > 1e-5 : sig[i] *= 1e-10 # convert (A) to (m)
# Binary parameters in the order of k12, k13, .. k1N, k23, k24, ..
f.readline()
Kij = list( map( float, f.readline().split() ) )
# Setup kij
# make 2d Zero array for Keps[i][j]
Keps = [ [0]*ncP for i in range(ncP) ]
for i in range(1,ncP) :
for j in range(i+1,ncP) :
Keps[i][j] = Kij.pop(0) # pop from the first
Keps[j][i] = Keps[i][j]
print("Kij's are ")
for i in range(1,ncP) :
for j in range(1,ncP) :
print("%9.4f" % Keps[i][j], end=" ")
print('')
print("")
# fill cross association by Walbach-Sandler combining rules
for i in range(1,ncP) :
for j in range(i+1,ncP) :
eAB[i][j] = (eAB[i][i]+eAB[j][j])/2
eAB[j][i] = eAB[i][j]
kAB[i][j] = sqrt( kAB[i][i]*kAB[j][j] ) * ( sqrt(sig[i]*sig[j]) * 2/(sig[i]+sig[j]) )**3
kAB[j][i] = kAB[i][j]
# Read cross association override Walbach-Sandler combining rules (if provided)
print( "cross association override in PC-SAFT \n i j kAiBj eAiBj (K)" )
f.readline()
novride = 0
while True :
try :
args = f.readline().split() # i j kAiBj eAiBj
i = int( args[0] )
j = int( args[1] )
kAB[i][j] = float( args[2] )
eAB[i][j] = float( args[3] )
kAB[j][i] = kAB[i][j]
eAB[j][i] = eAB[i][j]
print("%5d %5d %11.6f %8.2f " %( i, j, kAB[i][j], eAB[i][j] ) )
novride += 1
except :
break
print("The number of overrides is = %d \n" % novride )
# check the number of associating components including self- and induced- associations
if eos == 'PC-SAFT' :
print("checking associating components in PC-SAFT ..." )
ncas = 0 # number of associating components
for i in range(1,ncP) :
for j in range(1,ncP) :
if kAB[i][j] > 1e-10 :
print(" component %d is associating " %i )
ncas += 1
break # stop checking and examine next component
print("The number of associating components = %d \n" % ncas )
# read T, P, compositions
z = [0]*ncP
f.readline()
TypCal = int( f.readline() ) # 0 (Z), 1 (phi), 2 (BubbleP), 3 (DewP),4 (BubbleT), 5 (DewT), 6 (Flash)
#
f.readline()
args = list( map( float, f.readline().split() ))
T = args.pop(0)
P = args.pop(0)
for i in range(1,ncP) :
z[i] = args.pop(0)
# Echo type of calculation
if TypCal == 2 : # 2 (BubbleP)
print("Bubble P calculation at T = %7.2f (K) %7.2f (C) \n" % ( T, T-273.15 ) )
plot_title = "Bubble P calculation at T = %7.2f (K) \n" % T \
+ component[1] + " and " + component[2] + " \n" # for binary mixture
elif TypCal == 3 : # 3 (DewP)
print(" Dew P calculation at T = %7.2f (K) %7.2f (C) \n" % ( T, T-273.15 ) )
plot_title = " Dew P calculation at T = %7.2f (K) \n" % T \
+ component[1] + " and " + component[2] + " \n" # for binary mixture
elif TypCal == 4 : # 4 (BubbleT)
print("Bubble T calculation at P = %8.3f (bar) \n" % P )
plot_title = "Bubble T calculation at P = %8.3f (bar) \n" % P \
+ component[1] + " and " + component[2] + " \n" # for binary mixture
elif TypCal == 5 : # 5 (DewT)
print(" Dew T calculation at P = %8.3f (bar) \n" % P )
plot_title = " Dew T calculation at P = %8.3f (bar) \n" % P \
+ component[1] + " and " + component[2] + " \n" # for binary mixture
elif TypCal == 6 : # 6 (Flash)
print(" Flash calculation at T = %7.2f (K) and P = %8.3f (bar) \n " %( T, P ) )
# scaling fraction for initial liquid volume taken from PR eos
f.readline()
scalVl0 = float( f.readline() ) # 1.0 # 0.95 # 0.9
# scalVl0 = float( input("enter scaling fraction for initial liquid volume (e.g. 0.95) = ") )
print("scaling fraction for initial liquid volume = %7.2f \n" % scalVl0 )
# read experimental data for binary mixture if provided
Texp = []
Pexp = []
Xexp = []
Yexp = []
expdataheader = f.readline()
if expdataheader.find("y") > -1 : # if header contains "y" for gas composition
isTPxySet = True
isTPx_Set = False
nexp = 4
else :
isTPxySet = False
isTPx_Set = True
nexp = 3
print("-- Experimental data -- \n T (K) P (bar) x1 y1 ")
while True :
try :
line = f.readline().replace("\t"," ") # convert tab into space
#print(line)
args = line.split() # T, P, x, y
for i in range(nexp) :
args[i] = float( args[i] ) # T, P, x, y into numeric data
if expdataheader.find("(C)") > -1 : args[0] += 273.15 # C -> K
if expdataheader.find("(F)") > -1 : args[0] = (args[0]-32)/1.8 + 273.15 # F -> K
if expdataheader.find("(mmHg)")> -1 : args[1] *= 1.013 / 760 # mmHg -> bar
if expdataheader.find("(psi)") > -1 : args[1] *= 1.013 / 14.7 # psi -> bar
if expdataheader.find("(MPa)") > -1 : args[1] *= 10.0 # MPa -> bar
if expdataheader.find("(kPa)") > -1 : args[1] /= 100.0 # kPa -> bar
if isTPxySet : print(" %7.2f %8.3f %8.4f %8.4f " %( args[0], args[1], args[2], args[3] ) )
if isTPx_Set : print(" %7.2f %8.3f %8.4f " %( args[0], args[1], args[2] ) )
Texp.append( args[0] )
Pexp.append( args[1] )
Xexp.append( args[2] )
if isTPxySet : Yexp.append( args[3] )
#input("Press Enter to do .. ")
except :
break
# saving result for successive calculations
SavedResult = { }
SavedResult[ 'existVvap' ] = False
SavedResult[ 'existVliq' ] = False
return TypCal, T, P, z
#--[ Main program ]-----------------------------------------------------
global eos, R, component, Mw, Tc, Pc, omega, Tb, nc, ncP, Keps, plot_title, SavedResult
global Texp, Pexp, Xexp, Yexp, isTPxySet
global m, sig, epsk, kAB, eAB, ncas, dipol, xp
# unit conversion
m3TOcm3 = 1e6
print("=========================================")
print(" Vapor-Liquid Equilibrium of Mixture ")
print("=========================================")
TypCal, T, P, z = Read_data__initialize()
input("press Enter to proceed ")
print()
if TypCal == -1 : # test calculation for pressure
V = float( input("enter V (cm3/mol)? ") )
print("<< Calculation for given T, V, z's >> \n " )
print(" T = %10.3f (K) " % T )
print(" V = %10.4f (cm3/mol) " % V )
print(" z's = ", end=" " )
for i in range(1,ncP) :
print("%8.5f" %z[i], end=" ")
print("\n")
V = V * 1e-6
P, Z = Peos( V, T, z )
print(" P = %10.4f (bar) " % P )
print(" Z = %10.4f " % Z )
elif TypCal == 0 : # 0 (Z)
print("<< Calculation for given T, P, z's >> \n " )
# T, P, z's
print(" T = %10.3f (K) " % T )
print(" P = %10.4f (bar) " % P )
print(" z's = ", end=" " )
for i in range(1,ncP) :
print("%8.5f" %z[i], end=" ")
print("\n")
Vvap = Find_vol( P, T, z, 'vapor' )
Vliq = Find_vol( P, T, z, 'liquid' )
Zvap = P * Vvap / ( R * T )
Zliq = P * Vliq / ( R * T )
print(" compressibility factor molar volume (cm^3/mol) ")
print(" vapor liquid vapor liquid ")
print(" %8.4f %8.4f %10.2f %10.2f " % (Zvap, Zliq, Vvap*m3TOcm3, Vliq *m3TOcm3 ) )
elif TypCal == 1 : # 1 (phi)
print("<< Calculation for given T, P, z's >> \n " )
# T, P, z's
print(" T = %10.2f (K) " % T )
print(" P = %10.4f (bar) " % P )
print(" z's = ", end=" " )
for i in range(1,ncP) :
print("%8.5f" %z[i], end=" ")
print("\n")
Vvap = Find_vol( P, T, z, 'vapor' )
Vliq = Find_vol( P, T, z, 'liquid' )
phiV, fV = Fugacity( Vvap, T, z, P )
phiL, fL = Fugacity( Vliq, T, z, P )
print(" Vapor phase Liquid phase ")
print(" Comp fugacity (bar) phi fugacity (bar) phi " )
for i in range(1,ncP) :
print(" %d %10.4f %10.4f %10.4f %10.4f" \
%( i, fV[i], phiV[i], fL[i], phiL[i] ) )
elif TypCal == 2 : # 2 (BubbleP)
print("<< Pxy - Bubble point pressure calculation for given T, x >> \n " )
print(" T = %9.2f (K) " % T )
x = z[:]
WholeRange = False # single point calculation
if nc == 2 and z[1]+z[2] < 1e-10 : WholeRange = True # for example, zero compositions are given
if not WholeRange :
print(" x = ", end=" " )
for i in range(1,ncP) :
print("%8.4f" %x[i], end=" ")
print()
P, y = BubbleP_ini( T, x )
P, y, Vliq, Vvap, it = BubbleP( T, x, P, y )
print("\n Bubble point pressure = %10.5f (bar) \n" % P )
print(" comp liquid vapor K " )
for i in range(1,ncP) :
print(" %d %8.4f %8.4f %8.4f " %( i, x[i], y[i], y[i]/x[i] ) )
print()
Zvap = P * Vvap / ( R * T )
Zliq = P * Vliq / ( R * T )
print(" Z %8.4f %8.4f \n" %( Zliq, Zvap ) )
print(" Number of iterations = %d " % it )
else : # WholeRange from x = 0 to 1 only for binary mixture
# list for plot
PP = []
XX = []
YY = []
nstep = 200
X = [ i/nstep for i in range(nstep+1) ]
# initial guess at x1 = 0
x[1] = X[0]
x[2] = 1 - x[1]
P, y = BubbleP_ini( T, x )
print()
print(" P (bar) x1 y1 iteration")
print("-----------------------------------------")
for x[1] in X :
x[2] = 1 - x[1]
# P, y, Vliq, Vvap, it = BubbleP( T, x, P, y )
try : P, y, Vliq, Vvap, it = BubbleP( T, x, P, y )
except : break
# Check error for single phase close to critical point
if it == None : break
if it > 999 : break # Convergence error occurs
print("%10.5f %8.4f %8.4f %d " %( P, x[1], y[1], it ) )
SavedResult[ 'existVvap' ] = True
SavedResult[ 'existVliq' ] = True
SavedResult[ 'Vvap' ] = Vvap
SavedResult[ 'Vliq' ] = Vliq
PP.append( P )
XX.append( x[1] )
YY.append( y[1] )
# plot Pxy diagram
plt.plot( XX, PP, linestyle='-' )
plt.plot( YY, PP, linestyle='--')
plt.plot( Xexp, Pexp, 'bo')
if isTPxySet : plt.scatter( Yexp, Pexp, s=40, facecolors='none', edgecolors ='r' )
plt.xlabel(r'$x_1, y_1$')
plt.xlim(0,1)
plt.ylabel('P (bar)')
plt.title( plot_title, fontsize=11 )
plt.show()
elif TypCal == 3 : # 3 (DewP)
print("<< Pxy - Dew point pressure calculation for given T, y >> \n " )
print(" T = %9.2f (K) " % T )
y = z[:]
WholeRange = False # single point calculation
if nc == 2 and z[1]+z[2] < 1e-10 : WholeRange = True
if not WholeRange :
print(" y = ", end=" " )
for i in range(1,ncP) :
print("%8.4f" % y[i], end=" ")
print("\n")
P, x = DewP_ini( T, y )
P, x, Vliq, Vvap, it = DewP( T, y, P, x )
print("\n Dew point pressure = %10.5f (bar) \n" % P )
print(" comp vapor liquid K " )
for i in range(1,ncP) :
print(" %d %8.4f %8.4f %8.4f " %( i, y[i], x[i], y[i]/x[i] ) )
print()
Zvap = P * Vvap / ( R * T )
Zliq = P * Vliq / ( R * T )
print(" Z %8.4f %8.4f \n" %( Zvap, Zliq ) )
print(" Number of iterations = %d " % it )
else : # WholeRange from x = 0 to 1 only for binary mixture
# list for plot
PP = []
XX = []
YY = []
nstep = 100
Y = [ i/nstep for i in range(nstep+1) ]
# initial guess at y1 = 0
y[1] = Y[0]
y[2] = 1 - y[1]
P, x = DewP_ini( T, y )
print()
print(" P (bar) x1 y1 iteration")
print("-----------------------------------------")
for y[1] in Y :
y[2] = 1 - y[1]
P, x, Vliq, Vvap, it = DewP( T, y, P, x )
# Check error for single phase close to critical point
if it == None : break
if it > 999 : break # Convergence error occurs
print("%10.5f %8.4f %8.4f %d " %( P, x[1], y[1], it ) )
SavedResult[ 'existVvap' ] = True
SavedResult[ 'existVliq' ] = True
SavedResult[ 'Vvap' ] = Vvap
SavedResult[ 'Vliq' ] = Vliq
PP.append( P )
XX.append( x[1] )
YY.append( y[1] )
# plot Pxy diagram
plt.plot( XX, PP, linestyle='-' )
plt.plot( YY, PP, linestyle='--')
plt.plot( Xexp, Pexp, 'bo')
if isTPxySet : plt.scatter( Yexp, Pexp, s=40, facecolors='none', edgecolors ='r' )
plt.xlabel(r'$x_1, y_1$')
plt.xlim(0,1)
plt.ylabel('P (bar)')
plt.title( plot_title, fontsize=11 )
plt.show()
elif TypCal == 4 : # 4 (BubbleT)
print("<< Txy - Bubble point temperature calculation for given P, x >> \n " )
print(" P = %9.3f (bar) " % P )
x = z[:]
WholeRange = False # single point calculation
if nc == 2 and z[1]+z[2] < 1e-10 : WholeRange = True
if not WholeRange :
print(" x = ", end=" " )
for i in range(1,ncP) :
print("%8.4f" % x[i], end=" ")
print()
T, y = BubbleT_ini( P, x )
T, y, Vliq, Vvap, it = BubbleT( P, x, T, y )
print("\n Bubble point temperature = %8.2f (K) \n" % T )
print(" comp liquid vapor K " )
for i in range(1,ncP) :
print(" %d %8.4f %8.4f %8.4f " %( i, x[i], y[i], y[i]/x[i] ) )
print()
Zvap = P * Vvap / ( R * T )
Zliq = P * Vliq / ( R * T )
print(" Z %8.4f %8.4f \n" %( Zliq, Zvap ) )
print(" Number of iterations = %d " % it )
else : # WholeRange from x = 0 to 1 only for binary mixture
# list for plot
TT = []
XX = []
YY = []
nstep = 200
X = [ i/nstep for i in range(nstep+1) ]
# initial guess at x1 = 0
x[1] = X[0]
x[2] = 1 - x[1]
T, y = BubbleT_ini( P, x )
print()
print(" T (K) x1 y1 iteration")
print("-----------------------------------------")
for x[1] in X :
x[2] = 1 - x[1]
T, y, Vliq, Vvap, it = BubbleT( P, x, T, y )
# Check error for single phase close to critical point
if it == None : break
if it > 999 : break # Convergence error occurs
print(" %8.2f %8.4f %8.4f %d " %( T, x[1], y[1], it ) )
SavedResult[ 'existVvap' ] = True
SavedResult[ 'existVliq' ] = True
SavedResult[ 'Vvap' ] = Vvap
SavedResult[ 'Vliq' ] = Vliq
TT.append( T )
XX.append( x[1] )
YY.append( y[1] )
# plot Pxy diagram
plt.plot( XX, TT, linestyle='-' )
plt.plot( YY, TT, linestyle='--')
plt.plot( Xexp, Texp, 'bo')
if isTPxySet : plt.scatter( Yexp, Texp, s=40, facecolors='none', edgecolors ='r' )
plt.xlabel(r'$x_1, y_1$')
plt.xlim(0,1)
plt.ylabel('T (K)')
plt.title( plot_title, fontsize=11 )
plt.show()
elif TypCal == 5 : # 5 (DewT)
print("<< Txy - Dew point temperature calculation for given P, y >> \n " )
print(" P = %9.3f (bar) " % P )
y = z[:]
WholeRange = False # single point calculation
if nc == 2 and z[1]+z[2] < 1e-10 : WholeRange = True
if not WholeRange :
print(" y = ", end=" " )
for i in range(1,ncP) :
print("%8.4f" % y[i], end=" ")
print()
T, x = DewT_ini( P, y )
T, x, Vliq, Vvap, it = DewT( P, y, T, x )
print("\n Dew point temperature = %8.2f (K) \n" % T )
print(" comp vapor liquid K " )
for i in range(1,ncP) :
print(" %d %8.4f %8.4f %8.4f " %( i, y[i], x[i], y[i]/x[i] ) )
print()
Zvap = P * Vvap / ( R * T )
Zliq = P * Vliq / ( R * T )
print(" Z %8.4f %8.4f \n" %( Zvap, Zliq ) )
print(" Number of iterations = %d " % it )
else : # WholeRange from x = 0 to 1 only for binary mixture
# list for plot
TT = []
XX = []
YY = []
nstep = 100
Y = [ i/nstep for i in range(nstep+1) ]
# initial guess at y1 = 0
y[1] = Y[0]
y[2] = 1 - y[1]
T, x = DewT_ini( P, y )
print()
print(" T (K) x1 y1 iteration")
print("-----------------------------------------")
for y[1] in Y :
y[2] = 1 - y[1]
T, x, Vliq, Vvap, it = BubbleT( P, y, T, x )
# Check error for single phase close to critical point
if it == None : break
if it > 999 : break # Convergence error occurs
print(" %8.2f %8.4f %8.4f %d " %( T, x[1], y[1], it ) )
SavedResult[ 'existVvap' ] = True
SavedResult[ 'existVliq' ] = True
SavedResult[ 'Vvap' ] = Vvap
SavedResult[ 'Vliq' ] = Vliq
TT.append( T )
XX.append( x[1] )
YY.append( y[1] )
# plot Pxy diagram
plt.plot( XX, TT, linestyle='-' )
plt.plot( YY, TT, linestyle='--')
plt.plot( Xexp, Texp, 'bo')
if isTPxySet : plt.scatter( Yexp, Texp, s=40, facecolors='none', edgecolors ='r' )
plt.xlabel(r'$x_1, y_1$')
plt.xlim(0,1)
plt.ylabel('T (K)')
plt.title( plot_title, fontsize=11 )
plt.show()
elif TypCal == 6 : # 6 (Flash)
print("<< Flash calculation for given T, P, z >> \n " )
print(" T = %9.2f (K) " % T )
print(" P = %9.3f (bar) " % P )
print(" z = ", end=" " )
for i in range(1,ncP) : print("%8.4f" % z[i], end=" ")
print("\n")
# Get bubble pressure and dew pressure, and check whether Pdew < P < Pbubble
Pb, yb = BubbleP_ini( T, z )
Pb, yb, Vliq, Vvap, it = BubbleP( T, z, Pb, yb )
print(" bubble P = %9.3f (bar) " % Pb )
Pd, xd = DewP_ini( T, z )
Pd, xd, Vliq, Vvap, it = DewP( T, z, Pd, xd )
print(" dew P = %9.3f (bar) \n" % Pd )
if P < Pd or P > Pb :
print(" Pressure is not in valid range! \n Taking mean value ... \n" )
P = ( Pd + Pb ) / 2
print(" P = %9.3f (bar) " % P )
else :
print(" Pressure is valid: P(dew) < P < P(bubble)")
# initial guess
v = ( Pb - P ) / ( Pb - Pd )
x = [0]*ncP
y = [0]*ncP
for i in range(1,ncP) : x[i] = (1-v) * z[i] + v *xd[i]
for i in range(1,ncP) : y[i] = (1-v) *yb[i] + v * z[i]
x, y, K, Vliq, Vvap, v, it = Flash( P, T, z, x, y, v )
print()
print(" comp feed liquid vapor K " )
for i in range(1,ncP) :
print(" %d %8.4f %8.4f %8.4f %8.4f " %( i, z[i], x[i], y[i], K[i] ) )
print()
Zvap = P * Vvap / ( R * T )
Zliq = P * Vliq / ( R * T )
print(" fraction 1.0000 %8.4f %8.4f \n" %( 1-v, v ) )
print(" Z %8.4f %8.4f \n" %( Zliq, Zvap ) )
print(" Number of iterations = %d \n" % it )
input("Press Enter to quit ")
pcsaft_mix.py¶
In [6]:
!type pcsaft_mix.py
# PC-SAFT mixture
from math import *
import numpy as np
from scipy import linalg
#-----------------------------------------------------------------------
def pcsaft( ZorPhi, V, T, x, m, sig, epsk, kAB, eAB, nc, Keps, ncas, dipol, xp ) :
ncP = nc + 1
Pi = 3.1415926536
Nav = 6.02214179e23
rhom = 1 / V # molar density (mol/m3)
rho = rhom * Nav # number density (molecules/m3)
M = 0 # mean chain length
d = [0]*ncP # effective segmemt diameter (m)
t = [0]*4
g = [0]*ncP
roDg = [0]*ncP
for i in range(1,ncP) :
M += x[i] * m[i]
d[i] = sig[i] * ( 1 - 0.12*exp( -3*epsk[i]/T ) )
t[0] += Pi/6*rho* x[i] * m[i]
t[1] += Pi/6*rho* x[i] * m[i] * d[i]
t[2] += Pi/6*rho* x[i] * m[i] * d[i]**2
t[3] += Pi/6*rho* x[i] * m[i] * d[i]**3
# Hard-sphere term m * (Zhs-1) ...
Zhc = M * ( t[3]/(1-t[3]) + 3*t[1]*t[2]/t[0]/(1-t[3])**2 \
+ t[2]**3 *( 3 - t[3] ) /t[0]/(1-t[3])**3 )
# Hard chain term
for i in range(1,ncP) :
g[i] = 1/(1-t[3]) + d[i]/2 * 3*t[2]/(1-t[3])**2 \
+ (d[i]/2)**2 * 2*t[2]**2/(1-t[3])**3
roDg[i] = t[3]/(1-t[3])**2 \
+ d[i]/2 * ( 3*t[2]/(1-t[3])**2 + 6*t[2]*t[3]/(1-t[3])**3 ) \
+ (d[i]/2)**2 * ( 4*t[2]**2/(1-t[3])**3 + 6*t[2]**2 *t[3]/(1-t[3])**4 )
Zhc += - x[i] * (m[i]-1) * roDg[i] / g[i]
# Dispersion term
y = t[3] # packing fraction
a0 = [ 0.9105631445, 0.6361281449, 2.6861347891, -26.547362491, 97.759208784, -159.59154087, 91.297774084]
a1 = [ -0.3084016918, 0.1860531159, -2.5030047259, 21.419793629, -65.255885330, 83.318680481, -33.746922930]
a2 = [ -0.0906148351, 0.4527842806, 0.5962700728, -1.7241829131, -4.1302112531, 13.776631870, -8.6728470368]
b0 = [ 0.7240946941, 2.2382791861, -4.0025849485, -21.003576815, 26.855641363, 206.55133841, -355.60235612]
b1 = [ -0.5755498075, 0.6995095521, 3.8925673390, -17.215471648, 192.67226447, -161.82646165, -165.20769346]
b2 = [ 0.0976883116, -0.2557574982, -9.1558561530, 20.642075974, -38.804430052, 93.626774077, -29.666905585]
a = [ ]
b = [ ]
I1 = 0
I2 = 0
DyI1 = 0
DyI2 = 0
for i in range(7): # [0,1,2 ... 6]
a.append( a0[i] + (M-1)/M*a1[i] + (M-1)*(M-2)/M**2 *a2[i] )
b.append( b0[i] + (M-1)/M*b1[i] + (M-1)*(M-2)/M**2 *b2[i] )
I1 += a[i] * y**i # used for Adisp
I2 += b[i] * y**i
DyI1 += (i+1)* a[i] * y**i
DyI2 += (i+1)* b[i] * y**i
C1 = 1/( 1 + M*(8*y-2*y**2)/(1-y)**4 - (M-1)*(20*y-27*y**2+12*y**3-2*y**4)/(1-y)**2/(2-y)**2 )
# C2 = DC1 = (d/dy) C1
C2 = -C1*C1 * ( M*(-4*y**2+20*y+8)/(1-y)**5 - (M-1)*(2*y**3+12*y**2-48*y+40)/(1-y)**3/(2-y)**3 )
m2ekTsig3 = 0
m2ekT2sig3 = 0
for i in range(1,ncP) :
for j in range(1,ncP) :
sig3 = ( ( sig[i] + sig[j] ) / 2 )**3
ekT = sqrt( epsk[i] * epsk[j] ) * (1-Keps[i][j]) / T
m2ekTsig3 += x[i]*x[j]*m[i]*m[j]* ekT * sig3
m2ekT2sig3 += x[i]*x[j]*m[i]*m[j]* ekT**2 * sig3
Zdisp = - 2*Pi*rho *DyI1 *m2ekTsig3 - Pi*rho *M *( C1*DyI2 + C2*y*I2 ) * m2ekT2sig3
# Association term
Zassoc = 0
# ncas = 2 # for test
# single associaing components
if ncas == 1 :
# find the only associating component I
for i in range(1,ncP) :
if eAB[i][i] > 1e-10 : I = i
F = exp( eAB[I][I] / T ) - 1
K = sig[I]**3 * kAB[I][I]
Del = g[I] * F * K
if x[I] < 1e-10 : # to avoid math error for pure substance
X = 1
Zassoc = 0
else :
X = ( -1 + sqrt( 1 + 4*rho*x[I]*Del) ) / ( 2*rho*x[I]*Del )
roDDel = roDg[I] * F * K
DX = -x[I] * X**3 / ( 1 + rho*x[I]*X**2 *Del ) * ( Del + roDDel )
Zassoc = x[I] * 2* ( 1/X - 1/2 ) * rho * DX
#print(X,rho*DX) # for check
# multiple associaing components
elif ncas >= 2 :
F = [ [0]*ncP for i in range(ncP) ]
K = [ [0]*ncP for i in range(ncP) ]
DEL = [ [0]*ncP for i in range(ncP) ]
roDDEL = [ [0]*ncP for i in range(ncP) ]
for i in range(1,ncP) :
for j in range(1,ncP) :
F[i][j] = exp( eAB[i][j] / T ) - 1
K[i][j] = ( ( sig[i]+sig[j] )/2 )**3 * kAB[i][j]
Gij = 1/(1-t[3]) + d[i]*d[j]/(d[i]+d[j]) * 3*t[2]/(1-t[3])**2 \
+ ( d[i]*d[j]/(d[i]+d[j]) )**2 * 2*t[2]**2/(1-t[3])**3
DEL[i][j] = Gij * F[i][j] * K[i][j]
roDGij = t[3]/(1-t[3])**2 \
+ d[i]*d[j]/(d[i]+d[j]) * ( 3*t[2]/(1-t[3])**2 + 6*t[2]*t[3]/(1-t[3])**3 ) \
+ ( d[i]*d[j]/(d[i]+d[j]) )**2 * ( 4*t[2]**2/(1-t[3])**3 + 6*t[2]**2 *t[3]/(1-t[3])**4 )
roDDEL[i][j] = roDGij * F[i][j] * K[i][j]
# Find X's by iteration (successive substitution)
X = [1]*ncP # defaul value X = 1 for non-associating components
# guess of X for associating component
for i in range(1,ncP) :
if eAB[i][i] > 1e-10 and x[i] > 1e-10 :
X[i] = ( -1 + sqrt( 1 + 4*rho*x[i]*DEL[i][i] ) ) / ( 2*rho*x[i]*DEL[i][i] )
method = 0 # successive substitution
#method = 1 # Newton method
if method == 0 : # successive substitution
while 1 :
#print(X[1],X[2]) # to check
Xold = X[:]
err = 0
for i in range(1,ncP) :
sum = 1
for j in range(1,ncP) :
sum += rho * x[j] * X[j] * DEL[i][j]
X[i] = 1 / sum
err += abs( X[i]-Xold[i] )
if err/ncP < 1e-5 : break
elif method == 1 : # Newton method
dX = [0]*ncP
while 1 :
#print(X) # to check
A = [ [0]*nc for i in range(nc) ] # from 0, 1, 2 ...
B = [0]*nc
for i in range(1,ncP) :
i_ = i - 1
A[i_][i_] = 1
B[i_] = 1 - X[i]
for j in range(1,ncP) :
A[i_][i_] += rho * x[j] * X[j] * DEL[i][j]
B[i_] += - rho * x[j] * X[i] * X[j] * DEL[i][j]
j_ = j - 1
A[i_][j_] += rho * x[j] * X[i] * DEL[i][j]
# Solve matrix equation
u = linalg.solve(A, B)
dX[1:] = u[:]
# update X and check error
err = 0
for i in range(1,ncP) :
X[i] += dX[i]
err += abs( dX[i] )
if err/ncP < 1e-5 : break
# end if
# Find DX # = dX / d rho
DX = [0]*ncP
# method = 0 # by iteration
method = 1 # by matrix
if method == 0 : # method = 0 by iteration (successive substitution)
while 1 :
DXold = DX[:]
err = 0
for i in range(1,ncP) :
sum = 0
for j in range(1,ncP) :
sum += x[j] * X[j] * ( DEL[i][j] + roDDEL[i][j] )
if j == i : continue
sum += rho * x[j] * DX[j] * DEL[i][j]
DX[i] = - X[i]**2 / ( 1 + rho * x[i] * X[i]**2 *DEL[i][i] ) * sum
err += abs( DX[i]-DXold[i] )
if err/ncP < 1e-6 : break
elif method == 1 : # method == 1 by solving matrix
A = [ [0]*nc for i in range(nc) ] # from 0, 1, 2 ...
B = [0]*nc
for i in range(1,ncP) :
i_ = i - 1
A[i_][i_] = 1
B[i_] = 0
for j in range(1,ncP) :
j_ = j - 1
A[i_][j_] += X[i]**2 * rho * x[j] * DEL[i][j]
B[i_] += - X[i]**2 * x[j]*X[j]*( DEL[i][j] + roDDEL[i][j] )
# Solve matrix equation
u = linalg.solve(A, B)
DX[1:] = u[:]
# end if method
Zassoc = 0
for i in range(1,ncP) :
Zassoc += x[i]* 2* ( 1/X[i] - 1/2 ) * rho * DX[i]
#print(X[i],rho*DX[i]) # for check
#-- end of association --
# polar term - Jog and Chapman
Zpolar = 0
nc_pol = 0 # number of polar components
mu2k = [0]*ncP
for i in range(1,ncP) :
if dipol[i] > 1e-10 :
nc_pol += 1
mu2k[i] = dipol[i]**2 / 1.380572e26 # dipol^2 (D^2) / k (m3 K)
if nc_pol > 0 :
dd = [ [0]*ncP for i in range(ncP) ]
for i in range(1,ncP) :
for j in range(1,ncP) :
dd[i][j] = ( d[i] + d[j] ) / 2
rosta = 6/Pi * t[3]
I2p = ( 1 -0.3618*rosta -0.3205*rosta**2 +0.1078*rosta**3 ) / ( 1 -0.5236*rosta )**2
I3p = ( 1 +0.62378*rosta -0.11658*rosta**2 ) / ( 1 -0.59056*rosta +0.20059*rosta**2 )
DI2p = ( 0.6854 -0.83043848*rosta +0.3234*rosta**2 -0.05644408*rosta**3 ) / ( 1 -0.5236*rosta )**3
DI3p = ( 1.21434 -0.63434*rosta -0.0562765454*rosta**2 ) / ( 1 -0.59056*rosta +0.20059*rosta**2 )**2
sumij = 0
for i in range(1,ncP) :
for j in range(1,ncP) :
sumij += x[i]*x[j] *m[i]*m[j] *xp[i]*xp[j] *mu2k[i]*mu2k[j] / dd[i][j]**3
A2p = - 2 * Pi / 9 * rho / T**2 * sumij * I2p
roDA2p = A2p - 2 * Pi / 9 * rho / T**2 * sumij * rosta * DI2p
sumijl = 0
for i in range(1,ncP) :
for j in range(1,ncP) :
for l in range(1,ncP) :
sumijl += x[i]*x[j]*x[l] *m[i]*m[j]*m[l] *xp[i]*xp[j]*xp[l] \
*mu2k[i]*mu2k[j]*mu2k[l] /( dd[i][j]*dd[j][l]*dd[i][l] )
A3p = 5 * Pi**2 / 162 * rho**2 / T**3 * sumijl * I3p
roDA3p = 2*A3p + 5 * Pi**2 / 162 * rho**2 / T**3 * sumijl * rosta * DI3p
if abs( A2p ) < 1e-10 : # nonpolar pure component
Zpolar = 0
else :
Zpolar = 2/( 1 - A3p/A2p ) * roDA2p - 1/( 1 - A3p/A2p )**2 * ( roDA2p - roDA3p )
#-- end of polar term
Z = 1 + Zhc + Zdisp + Zassoc + Zpolar
if ZorPhi == 'Z' : return Z # done with Z
# Fugacity coefficient
Ahs = 1/t[0] * ( 3*t[1]*t[2]/(1-t[3]) + t[2]**3 /t[3]/(1-t[3])**2 \
+ ( t[2]**3/t[3]**2 - t[0] )* log(1-t[3]) )
Ahc = M * Ahs
for i in range(1,ncP) : Ahc += -x[i] * (m[i]-1) * log( g[i] )
tx = [ [0]*ncP for i in range(4) ]
for k in range(1,ncP) :
tx[0][k] = Pi/6 *rho *m[k]
tx[1][k] = Pi/6 *rho *m[k] * d[k]
tx[2][k] = Pi/6 *rho *m[k] * d[k]**2
tx[3][k] = Pi/6 *rho *m[k] * d[k]**3
Ahsx = [0]*ncP
for k in range(1,ncP) :
Ahsx[k] = -tx[0][k]/t[0]*Ahs + 1/t[0] * ( 3*( tx[1][k]*t[2] + t[1]*tx[2][k] )/(1-t[3]) \
+ 3*t[1]*t[2]*tx[3][k]/(1-t[3])**2 + 3*t[2]**2 *tx[2][k]/t[3]/(1-t[3])**2 \
+ t[2]**3 *tx[3][k] *(3*t[3]-1) /t[3]**2 /(1-t[3])**3 \
+ ( ( 3*t[2]**2 *tx[2][k]*t[3] - 2*t[2]**3 *tx[3][k] )/t[3]**3 - tx[0][k] ) \
* log(1-t[3]) + ( t[0] - t[2]**3 / t[3]**2 ) * tx[3][k]/(1-t[3]) )
gx = [ [0]*ncP for i in range(ncP) ]
for i in range(1,ncP) :
for k in range(1,ncP) :
gx[i][k] = tx[3][k]/(1-t[3])**2 \
+ d[i]/2 *( 3*tx[2][k]/(1-t[3])**2 + 6*t[2]*tx[3][k]/(1-t[3])**3 ) \
+ (d[i]/2)**2 * ( 4*t[2]*tx[2][k]/(1-t[3])**3 + 6*t[2]**2 *tx[3][k]/(1-t[3])**4 )
Ahcx = [0]*ncP
for k in range(1,ncP) :
Ahcx[k] = m[k] * Ahs + M * Ahsx[k] - (m[k]-1) * log( g[k] ) # with correction
for i in range(1,ncP) : Ahcx[k] += -x[i] * (m[i]-1) / g[i] * gx[i][k]
# add to fugacity coefficient
lnphi = [0]*ncP
for k in range(1,ncP) :
lnphi[k] = Z-1 - log(Z) + Ahc + Ahcx[k]
for j in range(1,ncP) : lnphi[k] += -x[j] * Ahcx[j]
# Dispersion term
Adisp = -2*Pi *rho *I1 * m2ekTsig3 -Pi *rho *M *C1 *I2 *m2ekT2sig3
m2ekTsig3x = [0]*ncP
m2ekT2sig3x = [0]*ncP
for k in range(1,ncP) :
for j in range(1,ncP) :
sig3 = ( ( sig[k] + sig[j] ) / 2 )**3
ekT = sqrt( epsk[k] * epsk[j] ) * (1-Keps[k][j]) / T
m2ekTsig3x[k] += 2*m[k]*x[j]*m[j]* ekT * sig3
m2ekT2sig3x[k] += 2*m[k]*x[j]*m[j]* ekT**2 * sig3
C1x = [0]*ncP
for k in range(1,ncP) :
C1x[k] = C2 * tx[3][k] - C1**2 *( m[k]*(8*y -2*y**2)/(1-y)**4 \
- m[k]*(20*y -27*y**2 +12*y**3 -2*y**4)/(1-y)**2/(2-y)**2 )
ax = [ [0]*ncP for i in range(7) ]
bx = [ [0]*ncP for i in range(7) ]
for i in range(7) :
for k in range(1,ncP) :
ax[i][k] = m[k]/M**2 *a1[i] + m[k]/M**2 *(3-4/M) *a2[i]
bx[i][k] = m[k]/M**2 *b1[i] + m[k]/M**2 *(3-4/M) *b2[i]
I1x = [0]*ncP
I2x = [0]*ncP
for k in range(1,ncP) :
for i in range(7) :
I1x[k] += a[i] * i * tx[3][k] * y**(i-1) + ax[i][k] * y**i
I2x[k] += b[i] * i * tx[3][k] * y**(i-1) + bx[i][k] * y**i
Adispx = [0]*ncP
for k in range(1,ncP) :
Adispx[k] = - 2*Pi*rho * ( I1x[k] * m2ekTsig3 + I1 * m2ekTsig3x[k] ) \
- Pi*rho *( ( m[k]*C1*I2 + M*C1x[k]*I2 + M*C1*I2x[k] ) *m2ekT2sig3 \
+ M *C1 *I2 *m2ekT2sig3x[k] )
# add to fugacity coefficient
for k in range(1,ncP) :
lnphi[k] += Adisp + Adispx[k]
for j in range(1,ncP) : lnphi[k] += -x[j] * Adispx[j]
# Association term
# single associaing components of I
if ncas == 1 :
Aassoc = x[I] * 2* ( log(X) - X/2 + 1/2 )
Delx = [0]*ncP
Xx = [0]*ncP
for k in range(1,ncP) :
Delx[k] = gx[I][k] * F * K
Xx[k] = -X**2/( 1 + rho*x[I]*X**2*Del)*( rho*x[I]*X*Delx[k] )
Xx[I] += -X**2/( 1 + rho*x[I]*X**2*Del)*( rho*X*Del )
Aassocx = [0]*ncP
for k in range(1,ncP) :
Aassocx[k] = x[I] * 2* ( 1/X - 1/2 ) * Xx[k]
Aassocx[I] += 2 * ( log(X) - X/2 + 1/2 )
# add to fugacity coefficient
for k in range(1,ncP) :
lnphi[k] += Aassoc + Aassocx[k]
for j in range(1,ncP) : lnphi[k] += -x[j] * Aassocx[j]
# multiple associaing components
elif ncas >= 2 :
DELx = [ [ [0]*ncP for j in range(ncP) ] for i in range(ncP) ]
for i in range(1,ncP) :
for j in range(1,ncP) :
for k in range(1,ncP) :
Gijxk = tx[3][k]/(1-t[3])**2 \
+ d[i]*d[j]/(d[i]+d[j]) *( 3*tx[2][k]/(1-t[3])**2 + 6*t[2]*tx[3][k]/(1-t[3])**3 ) \
+ ( d[i]*d[j]/(d[i]+d[j]) )**2 * ( 4*t[2]*tx[2][k]/(1-t[3])**3 + 6*t[2]**2 *tx[3][k]/(1-t[3])**4 )
DELx[i][j][k] = Gijxk * F[i][j] * K[i][j]
# Find (X)x
Xx = [ [0]*ncP for i in range(ncP) ]
# method = 0 # by iteration
method = 1 # by matrix
if method == 0 : # method = 0 by iteration (successive substitution)
while 1 :
Xxold = Xx[:]
err = 0
for i in range(1,ncP) :
for k in range(1,ncP) :
sum = rho * X[k] * DEL[i][k]
for j in range(1,ncP) :
sum += rho * x[j] * X[j] * DELx[i][j][k]
if j == i : continue
sum += rho * x[j] * Xx[j][k] * DEL[i][j]
Xx[i][k] = - X[i]**2 / ( 1 + rho *x[i] *X[i]**2 *DEL[i][i] ) * sum
err += abs( Xx[i][k] - Xxold[i][k] )
if err/ncP/ncP < 1e-6 : break
elif method == 1 : # by solving matrix
A = [ [0]*nc**2 for i in range(nc**2) ] # from 0, 1, 2 ...
B = [0]*nc**2
for i in range(1,ncP) :
for k in range(1,ncP) :
ik_ = (i-1)*nc + k - 1
A[ik_][ik_] = 1
B[ik_] = - X[i]**2 * rho * X[k]* DEL[i][k]
for j in range(1,ncP) :
jk_ = (j-1)*nc + k - 1
A[ik_][jk_] += X[i]**2 * rho * x[j] * DEL[i][j]
B[ik_] += - X[i]**2 * rho * x[j] * X[j] * DELx[i][j][k]
# Solve matrix equation
u = linalg.solve(A, B)
# patch solution
for i in range(1,ncP) :
for k in range(1,ncP) :
ik_ = (i-1)*nc + k - 1
Xx[i][k] = u[ik_]
# end if method
Aassoc = 0
for i in range(1,ncP) :
Aassoc += x[i] * 2* ( log(X[i]) - X[i]/2 + 1/2 )
Aassocx = [0]*ncP
for k in range(1,ncP) :
for i in range(1,ncP) :
Aassocx[k] += x[i] * 2* ( 1/X[i] - 1/2 ) * Xx[i][k]
Aassocx[k] += 2* ( log(X[k]) - X[k]/2 + 1/2 )
# add to fugacity coefficient
for k in range(1,ncP) :
lnphi[k] += Aassoc + Aassocx[k]
for j in range(1,ncP) : lnphi[k] += -x[j] * Aassocx[j]
#-- end of association --
# polar term - Jog and Chapman
if nc_pol > 0 and abs( A2p ) > 1e-10 :
Apolar = A2p / ( 1 - A3p / A2p )
dI2px = [0]*ncP
dI3px = [0]*ncP
A2px = [0]*ncP
A3px = [0]*ncP
Apolx = [0]*ncP
for k in range(1,ncP) :
dI2px[k] = DI2p * rho * m[k] * d[k]**3
dI3px[k] = DI3p * rho * m[k] * d[k]**3
for k in range(1,ncP) :
sum = 0
for i in range(1,ncP) :
sum += x[i] *m[i]*m[k] *xp[i]*xp[k] *mu2k[i]*mu2k[k] / dd[i][k]**3
A2px[k] = - 2*Pi/9 * rho / T**2 * ( 2 * sum * I2p + sumij * dI2px[k] )
for k in range(1,ncP) :
sum = 0
for i in range(1,ncP) :
for j in range(1,ncP) :
sum += x[i]*x[j] *m[i]*m[j]*m[k] *xp[i]*xp[j]*xp[k] *mu2k[i]*mu2k[j]*mu2k[k] \
/( dd[i][j]*dd[j][k]*dd[i][k] )
A3px[k] = 5*Pi**2 /162 *rho**2 /T**3 * ( 3 *sum *I3p + sumijl * dI3px[k] )
for k in range(1,ncP) :
Apolx[k] = 2/(1 - A3p/A2p ) *A2px[k] - 1/(1 - A3p/A2p )**2 *( A2px[k] - A3px[k] )
# add to fugacity coefficient
for k in range(1,ncP) :
lnphi[k] += Apolar + Apolx[k]
for j in range(1,ncP) : lnphi[k] += -x[j] * Apolx[j]
#-- end of polar term
phi = [0]*ncP
for k in range(1,ncP) : phi[k] = exp( lnphi[k] )
return phi
