我使用Python对使用VdW混合规则下的Peng-Robinson状态方程进行了复现，在计算逸度系数时发现了报错问题，之后我寻找报错原因时，发现问题出现在压缩因子的计算上，文献中对压缩因子的定义是Z=PV/RT( 原文献：M.-T. Lee, S.-T. Lin / Fluid Phase Equilibria 254 (2007) 28–34)，由于是混合物，因此要用到混合物的第二维里系数来计算压缩因子。
混合物的第二维里系数计算公式为Bm=y1^2B11+2y1y2B12+y2^2B22，B11，B22分别是纯物质1和纯物质2的第二维里系数；B12代表混合物的性质，用以下经验式来计算：Bij=RTcij/Pcij*B(0)+wij*B(1)。
计算各临界参数的混合规则：Tcij=(1-kij)(TciTcj)^0.5，Vcij=((Vci^(1/3)+Vcj^(1/3))/2)^3，Zcij=(Zci+Zcj)/2，Pcij=(ZcijRTcij)/Vcij，
wij=(wi+wj)/2。(详见马沛生，李永红主编《化工热力学》(通用型) 第二版，第2章，26-27页)
计算压缩因子：Z=1+(BmP/RT)(问题出现在这里，我将混合物的各临界参数计算代入后，临界参数计算是没有错误的，单位也换算成统一的单位了，但是压缩因子的结果出乎意料，竟然达到了负几万的数值，难怪无法计算出逸度系数！)
我目前只会Python的基础知识，代码写的很长，友友们可以帮忙看一下代码的问题出现在哪里吗？谢谢友友们的解答！
要复现的文献的公式
复现文献公式的Python代码如下
#对M.-T.Lee,S.-T.Lin/FluidPhaseEquilibria254(2007)28–34文献中的PR-vdW状态方程复现
#n-Pentane(正戊烷)为组分1，1-Butanol(1-丁醇)为组分2
import numpy as np
import math
#J. Chem. Eng. Data 2005,50,1520-1524
#组分1和组分2的基本物性参数
R=8.314
Tc1=469.7
Tc2=563.3
Pc1=3.37*pow(10,6)
Pc2=4.414*pow(10,6)
Zc1=0.268
Zc2=0.258
w1=0.249
w2=0.589
Vc1=311*pow(10,-6)
Vc2=274*pow(10,-6)
#使用PR-vdW状态方程对J. Chem. Eng. Data 2005,50,1520-1524的数据的复现
#文献数据
T=468.15
P=np.array([8.33,12.02,15.23,18.52,21.88,25.32,28.68,30.91,32.06,32.57,32.86]) #压力单位是bar
x1=np.array([0,0.052,0.14,0.258,0.402,0.535,0.721,0.861,0.938,0.979,1])
y1=np.array([0,0.149,0.35,0.478,0.608,0.711,0.809,0.885,0.943,0.979,1])
x2=1-x1
y2=1-y1
P=P*pow(10,5)
#ki的计算
def ki(wi):
ki=0.37464+1.54226*wi-0.26992*pow(wi,2)
return ki
#αi(T)的计算
def αi_T(ki,T,Tci):
αi_T=pow(1+ki*(1-math.sqrt(T/Tci)),2)
return αi_T
#ai(T)的计算
def ai_T(R,Tci,Pci,αi_T):
ai_T=0.45724*pow(R,2)*pow(Tci,2)*αi_T/Pci
return ai_T
#bi(T)的计算
def bi(R,Tci,Pci):
bi=0.07780*R*Tci/Pci
return bi
#aij的计算
def aij(kij,ai,aj):
aij=math.sqrt(ai*aj)*(1-kij)
return aij
#混合物的临界温度
def Tcij(kij,Tci,Tcj):
Tcij=(1-kij)*math.sqrt(Tci*Tcj)
return Tcij
#对比温度
def Trij(T,Tcij):
Trij=T/Tcij
return Trij
#混合物的偏心因子
def wij(wi,wj):
wij=(w1+w2)/2
return wij
#混合物的临界压缩因子
def Zcij(Zci,Zcj):
Zcij=(Zci+Zcj)/2
return Zcij
#混合物的临界体积
def Vcij(Vci,Vcj):
Vcij=pow((pow(Vci,1/3)+pow(Vcj,1/3))/2,3)
return Vcij
#混合物的临界压力
def Pcij(Zcij,R,Tcij,Vcij):
Pcij=Zcij*R*Tcij/Vcij
return Pcij
#第二维里系数中B0的计算
def B0(Trij):
B0=0.083-0.422/pow(Trij,1.6)
return B0
#第二维里系数中B1的计算
def B1(Trij):
B1=0.139-0.172/pow(Trij,4.2)
return B1
#第二混合维里系数的计算
def Bij(R,Tcij,Pcij,wij,B0,B1):
Bij=R*Tcij*(B0+w1*B1)/Pcij
return Bij
#混合物总的第二维里系数的计算
def Bm(xi,xj,Bii,Bij,Bjj):
Bm=pow(xi,2)*Bii+2*xi*xj*Bij+pow(xj,2)*Bjj
return Bm
#压缩因子的计算
def Z(Bm,P,R,T):
Z=1+(Bm*P/R*T)
return Z
#混合物中常数项am的计算
def am(xi,xj,aii,aij,ajj):
am=pow(xi,2)*aii+2*xi*xj*aij+pow(xj,2)*ajj
return am
#混合物中常数项bm的计算
def bm(xi,xj,bi,bj):
bm=xi*bi+xj*bj
return bm
#逸度系数的计算
def phi(xi,xj,P,T,R,aii,aij,bi,Z,am,bm):
Φi=np.exp(bi/bm*(Z-1)-np.log(Z-(bm*P/R*T))-(am/(2*math.sqrt(2)*bm*R*T))*(2*(xi*aii+xj*aij)/am-bi/bm)*np.log(Z*R*T+(1+math.sqrt(2)*bm*P)/(Z*R*T+(1-math.sqrt(2))*bm*P)))
return Φi
#k1和k2的计算
k1=ki(w1)
k2=ki(w2)
#α1(T)和α2(T)的计算
α1_T=αi_T(k1,T,Tc1)
α2_T=αi_T(k2,T,Tc2)
#a1,a2和a12的计算
a1=ai_T(R,Tc1,Pc1,α1_T)
a2=ai_T(R,Tc2,Pc2,α2_T)
k12=0
a12=aij(k12,a1,a2)
#b1和b2的计算
b1=bi(R,Tc1,Pc1)
b2=bi(R,Tc2,Pc2)
#对比温度Tr1,Tr2和Tr12的计算
Tr1=Trij(T,Tc1)
Tr2=Trij(T,Tc2)
k12=0
Tc12=Tcij(k12,Tc1,Tc2)
Tr12=Trij(T,Tc12)
#Pc12的计算
Zc12=Zcij(Zc1,Zc2)
Vc12=Vcij(Vc1,Vc2)
Pc12=Pcij(Zc12,R,Tc12,Vc12)
#第二维里系数纯组分1，纯组分2以及组分1，2的计算
#B11的计算
B0_1=B0(Tr1)
B1_1=B1(Tr1)
B11=Bij(R,Tc1,Pc1,w1,B0_1,B1_1)
#B22的计算
B0_2=B0(Tr2)
B1_2=B1(Tr2)
B22=Bij(R,Tc2,Pc2,w2,B0_2,B1_2)
#B12的计算
B0_12=B0(Tr12)
B1_12=B1(Tr12)
B12=Bij(R,Tc12,Pc12,w2,B0_12,B1_12)
#总的第二维里系数的计算
Bm=Bm(x1,x2,B11,B12,B22)
#am的计算
a11=a1
a12=a12
a22=a2
am=am(x1,x2,a11,a12,a22)
#bm的计算
bm=bm(x1,x2,b1,b2)
#压缩因子的计算
Z=Z(Bm,P,R,T)
#组分1的液相逸度系数的计算
Φ1=phi(x1,x2,P,T,R,a11,a12,b1,Z,am,bm)
