RKS方程计算压缩因子

%RKS方程计算压缩因子子程序

function Z=RKSZ(x,Tc,Pc,w,T,P,c,qy)

% c=8;
R=8.314;
PI=3.1415926;

for i=1:c %//RKS方程参数
m=0.480+1.574*w(i)-0.176*w(i)*w(i);
W=1+m*(1-sqrt(T/Tc(i)));%---sqrt(ai(T))见郭天民书251页
a(i)=0.42748*W*W*R*R*Tc(i)*Tc(i)/Pc(i);%---αi(T)*ai(Tci)---αi(T)-- α-alpha, a--ai
b(i)=0.08664*R*Tc(i)/Pc(i);
end

am=0;bm=0;A=0;B=0;m=zeros(c,1);p=0;q=0;r=0;ang=0;zi=zeros(3,1);zj=0;
for i=1:c %混合参数
am=am+x(i)*sqrt(a(i));
bm=bm+x(i)*b(i);%--B
end
am=am^2;


A=am*P/(R*R*T*T);%--A*
B=bm*P/(R*T);%--B*
p=A-B-B*B;%一次项系数
q=-A*B;%常数项


xishu=[1 -1 p q];
zi=roots(xishu); %直用matlab函数求解最方便
if qy==1

if zi(1)>zi(2)
Z=zi(2);
else
Z=zi(1);
end

if Z>zi(3)
Z=zi(3);
end %求出三个根中最小值－液相压缩因子

else
if zi(1)>zi(2)
Z=zi(1);
else
Z=zi(2);
end

if ZZ=zi(3);
end %求出三个根中最大值－气相压缩因子
end