相平衡大作业(7)

附：原程序（MATLAB6.5语言程序）

%使用 Soase (RKS) 状态方程计算泡点温度及组成

%1：物性数据Propene 丙烯 Tc1=365k,Pc1=4.62Mpa,w1=0.148

%2：物性数据isobutane 异丁烷 Tc2=408.1k,Pc2=3.648Mpa,w2=0.176 tic %开始计时

Tc1=365;Pc1=4.620;w1=0.148; Tc2=408.1;Pc2=3.648;w2=0.176; %总压P=2.0Mpa P=2.000;

m1=0.48+1.574*w1-0.176*w1^2; m2=0.48+1.574*w2-0.176*w2^2; aTc1=0.42748*(8.314*Tc1)^2/Pc1; aTc2=0.42748*(8.314*Tc2)^2/Pc2; b1=0.0866*8.314*Tc1/Pc1; b2=0.0866*8.314*Tc2/Pc2; for x1=0:0.05:1 x2=1-x1;

for T=310:0.01:378.10 Tr1=T/Tc1; Tr2=T/Tc2;

%童景山书 aTr=[1+m(1-Tr1^0.5)]^2; aTr1=(1+m1*(1-Tr1^0.5))^2; aTr2=(1+m2*(1-Tr2^0.5))^2; a1=aTc1*aTr1; a2=aTc2*aTr2;

%烃烃二元体系混合规则,二元交互作用参数kij=0 a12=(a1*a2)^0.5; b12=(b1+b2)/2;

a=a1*x1.^2+2*a12*x1*x2+a2*x2.^2; b=b1*x1.^2+2*b12*x1*x2+b2*x2.^2;

A=P*a/(8.314*T)^2; B=P*b/T/8.314;

%ai的平均值用avi表示 av1=2*(x1*a1+x2*a12); av2=2*(x1*a12+x2*a2); %bi的平均值用bvi表示 bv1=2*(x1*b1+x2*b12)-b; bv2=2*(x1*b12+x2*b2)-b; H1=A-B-B^2; H2=A*B;

%求z^3-z^2+z*H1-H2=0的根 eq=[1,-1,H1,-H2]; z=roots(eq);

zv=max(z); zL=min(z);

%求逸度系数L,V,液相为L，汽相为V

L1=exp(bv1*(zL-1)/b-log(zL-B)-A/B*(av1/a-bv1/b)*log(1+B/zL)); L2=exp(bv2*(zL-1)/b-log(zL-B)-A/B*(av2/a-bv2/b)*log(1+B/zL)); V1=exp(bv1*(zv-1)/b-log(zv-B)-A/B*(av1/a-bv1/b)*log(1+B/zv)); V2=exp(bv2*(zv-1)/b-log(zv-B)-A/B*(av2/a-bv2/b)*log(1+B/zv)); %求汽相组成y1,y2 k1=L1/V1; k2=L2/V2; y1=k1*x1; y2=k2*x2; y=y1+y2; res=y-1;

if abs(res)<0.0001 x1

T y1 y2 L1 L2 V1 V2 zv zL

break; end end end

toc %终止计时

t=toc %计算共用时间（seconds）