第三题
%第三题 求气相摩尔分量
function [V]=GPMF(T,P)
A=[ 1 1.297 33.2 -0.22 65 0.637503;
2 4.600 190.6 0.008 99 0.093293;
3 4.884 305.4 0.098 148 0.017343;
4 22.05 647.3 0.344 56 0.168547;
5 11.28 405.6 0.250 72.5 0.002714;
6 8.937 373.2 0.100 98.5 0.012105;
7 8.309 324.6 0.12 81 0.000000283;
8 3.384 460.4 0.227 306 0.068494;];
Bi=A(:,2);%Pci
Ci=A(:,3);%Tci
D=A(:,4);%Wi
Z=A(:,6);%Zi
B=P./Bi;%公式2
C=T./Ci;%公式3
%公式1
K=exp(5.37*(1+D).*(1-1./C))./B;%Ki
%V的初值
V0=0.9;
while(1)
%公式4
F=sum((Z.*(K-1))./(1+((K-1)*V0)));
%公式5
F1=sum(-(Z.*(K-1).*(K-1))./((1+(K-1)*V0).*(1+(K-1)*V0)));
%公式6
V=V0-F/F1;
%公式7
dif=abs((V-V0)/V);
if (dif<0.001)
break
else
V0=V;
end
end
end
第四题
%第四题 求气相组成yV、液相组成xL（摩尔分数）
function [YV,XL]=MF(T,P,V)
A=[ 1 1.297 33.2 -0.22 65 0.637503;
2 4.600 190.6 0.008 99 0.093293;
3 4.884 305.4 0.098 148 0.017343;
4 22.05 647.3 0.344 56 0.168547;
5 11.28 405.6 0.250 72.5 0.002714;
6 8.937 373.2 0.100 98.5 0.012105;
7 8.309 324.6 0.12 81 0.000000283;
8 3.384 460.4 0.227 306 0.068494;];
Bi=A(:,2);%Pci
Ci=A(:,3);%Tci
D=A(:,4);%Wi
Z=A(:,6);%Zi
B=P./Bi;%公式2
C=T./Ci;%公式3
%公式1
K=exp(5.37*(1+D).*(1-1./C))./B;%Ki
%公式8
Yi=Z.*K./(1+(K-1).*V);
%公式9
YV=Yi./sum(Yi);
%公式10
Xi=Z./(1+(K-1).*V);
XL=Xi./sum(Xi);
end
第五题
%第五题 求气相混合物平均引力系数amV、平均斥力系数bmV，液相混合物平均引力系数amL、平均斥力系数bmL
function [amV,bmV,amL,bmL]=AC(T,P,YV,XL)
A=[ 1 1.297 33.2 -0.22 65 0.637503;
2 4.600 190.6 0.008 99 0.093293;
3 4.884 305.4 0.098 148 0.017343;
4 22.05 647.3 0.344 56 0.168547;
5 11.28 405.6 0.250 72.5 0.002714;
6 8.937 373.2 0.100 98.5 0.012105;
7 8.309 324.6 0.12 81 0.000000283;
8 3.384 460.4 0.227 306 0.068494;];
Bi=A(:,2);%Pci
Ci=A(:,3);%Tci
D=A(:,4);%Wi
Vi=A(:,5);%Vi
B=P./Bi;%公式2
C=T./Ci;%公式3
%公式12
%各组分引力系数
a=0.45724*82.06*82.06*C.*C./B;
%公式13
%各组分斥力系数
b=0.07780*82.06*C./B;
%公式15
m=0.37464+1.5422*D-0.26992*D.*D;
%公式14
%各组分温度可调函数
ai=power((1+m.*(1-power(C,0.5))),2);
%公式17和18
k=zeros(8,8);
for i=1:8
for j=1:8
if i==j
k(i,j)=0;
else
k(i,j)=1-power((2*power(Vi(i),1/6).*power(Vi(j),1/6))./(power(Vi(i),1/3)+power(Vi(j),1/3)),exp(1));
end
end
end
bmV=sum(YV.*b);
bmL=sum(XL.*b);
amV=0;
amL=0;
for i=1:8
for j=1:8
amV=amV+YV(i).*YV(j).*power(a(i).*a(j).*ai(i).*ai(j),0.5).*(1-k(i,j));
amL=amL+XL(i).*XL(j).*power(a(i).*a(j).*ai(i).*ai(j),0.5).*(1-k(i,j));
end
end
end
第六题
%第六题 气相混合物偏差系数ZmV，液相混合物偏差系数ZmL
function [ZmV,ZmL]=MDC(T,P,amV,bmV,amL,bmL)
%ZMV
%公式23
AmV=amV.*P./(82.06*82.06*T.*T);
%公式24
BmV=bmV.*P./(82.06*T);
ZmV0=1;
while(1)
%公式25
F_ZmV=power(ZmV0,3)-(1-BmV).*power(ZmV0,2)+(AmV-2*BmV-3*power(BmV,2)).*ZmV0-(AmV.*BmV-power(BmV,2)-power(BmV,3));
%公式26
F1_ZmV=3*power(ZmV0,2)-2*(1-BmV).*ZmV0+(AmV-2*BmV-3*power(BmV,2)).*ZmV0;
%公式27
ZmV=ZmV0-F_ZmV./F1_ZmV;
%公式28
dif1=abs((ZmV-ZmV0)/ZmV);
if (dif1<0.0005)
break
else
ZmV0=ZmV;
end
end
%ZmL
%公式29
AmL=amL.*P./(82.06*82.06*T.*T);
%公式30
BmL=bmL.*P./(82.06*T);
ZmL0=1;
while(1)
%公式31
F_ZmL=power(ZmL0,3)-(1-BmL).*power(ZmL0,2)+(AmL-2*BmL-3*power(BmL,2)).*ZmL0-(AmL.*BmL-power(BmL,2)-power(BmL,3));
%公式32
F1_ZmL=3*power(ZmL0,2)-2*(1-BmL).*ZmL0+(AmL-2*BmL-3*power(BmL,2)).*ZmL0;
%公式33
ZmL=ZmL0-F_ZmL./F1_ZmL;
%公式34
dif2=abs((ZmL-ZmL0)/ZmL);
if (dif2<0.0005)
break
else
ZmL0=ZmL;
end
end
end
第七题
%第七题 气相中各组分逸度fiV，液相中各组分逸度fiL
function [fiV,fiL]=FOEC(T,P,ZmV,ZmL,a,ai,k,YV,XL,bmV,bmL,AmV,AmL,BmV,BmL,amV,amL)
A=[ 1 1.297 33.2 -0.22 65 0.637503;
2 4.600 190.6 0.008 99 0.093293;
3 4.884 305.4 0.098 148 0.017343;
4 22.05 647.3 0.344 56 0.168547;
5 11.28 405.6 0.250 72.5 0.002714;
6 8.937 373.2 0.100 98.5 0.012105;
7 8.309 324.6 0.12 81 0.000000283;
8 3.384 460.4 0.227 306 0.068494;];
Bi=A(:,2);%Pci
Ci=A(:,3);%Tci
D=A(:,4);%Wi
Vi=A(:,5);%Vi
B=P./Bi;%公式2
C=T./Ci;%公式3
%各组分斥力系数
b=0.07780*82.06*C./B;
%fiV
%公式36
HiV=0;
for j=1:8
HiV=HiV+YV(j).*(1-k(j)).*power(a.*ai.*a(j).*ai(j),0.5);
end
fiV=YV.*P.*exp(b./bmV.*(ZmV-1)-log(ZmV-BmV)-AmV./(2.*sqrt(2).*BmV).*(2.*HiV./amV-b./bmV).*(log((ZmV+2.414.*BmV)./(ZmV-0.414.*BmV))));
%fiL
HiL=0;
for j=1:8
HiL=HiL+XL(j).*(1-k(j)).*power(a.*ai.*a(j).*ai(j),0.5);
end
fiL=XL.*P.*exp(b./bmL.*(ZmL-1)-log(ZmL-BmL)-AmL./(2.*sqrt(2).*BmL).*(2.*HiL./amL-b./bmL).*(log((ZmL+2.414.*BmL)./(ZmL-0.414.*BmL))));
end
第八题
%第八题 热力学平衡检验
function [V,YV,XL]=TMET(T,P)
A=[ 1 1.297 33.2 -0.22 65 0.637503;
2 4.600 190.6 0.008 99 0.093293;
3 4.884 305.4 0.098 148 0.017343;
4 22.05 647.3 0.344 56 0.168547;
5 11.28 405.6 0.250 72.5 0.002714;
6 8.937 373.2 0.100 98.5 0.012105;
7 8.309 324.6 0.12 81 0.000000283;
8 3.384 460.4 0.227 306 0.068494;];
Bi=A(:,2);%Pci
Ci=A(:,3);%Tci
D=A(:,4);%Wi
Vi=A(:,5);%Vi
Z=A(:,6);%Zi
B=P./Bi;%公式2
C=T./Ci;%公式3
while(1)
%公式1
K=exp(5.37*(1+D).*(1-1./C))./B;%Ki
%V的初值
V0=0.9;
while(1)
%公式4
F=sum((Z.*(K-1))./(1+((K-1)*V0)));
%公式5
F1=sum(-(Z.*(K-1).*(K-1))./((1+(K-1)*V0).*(1+(K-1)*V0)));
%公式6
V=V0-F/F1;
%公式7
dif=abs((V-V0)/V);
if (dif<0.001)
break
else
V0=V;
end
end
%公式8
Yi=Z.*K./(1+(K-1).*V);
%公式9
YV=Yi./sum(Yi);
%公式10
Xi=Z./(1+(K-1).*V);
XL=Xi./sum(Xi);
%公式12
%各组分引力系数
a=0.45724*82.06*82.06*C.*C./B;
%公式13
%各组分斥力系数
b=0.07780*82.06*C./B;
%公式15
m=0.37464+1.5422*D-0.26992*D.*D;
%公式14
%各组分温度可调函数
ai=power((1+m.*(1-power(C,0.5))),2);
%公式17和18
k=zeros(8,8);
for i=1:8
for j=1:8
if i==j
k(i,j)=0;
else
k(i,j)=1-power((2*power(Vi(i),1/6).*power(Vi(j),1/6))./(power(Vi(i),1/3)+power(Vi(j),1/3)),exp(1));
end
end
end
bmV=sum(YV.*b);
bmL=sum(XL.*b);
amV=0;
amL=0;
for i=1:8
for j=1:8
amV=amV+YV(i).*YV(j).*power(a(i).*a(j).*ai(i).*ai(j),0.5).*(1-k(i,j));
amL=amL+XL(i).*XL(j).*power(a(i).*a(j).*ai(i).*ai(j),0.5).*(1-k(i,j));
end
end
%ZMV
%公式23
AmV=amV.*P./(82.06*82.06*T.*T);
%公式24
BmV=bmV.*P./(82.06*T);
ZmV0=1;
while(1)
%公式25
F_ZmV=power(ZmV0,3)-(1-BmV).*power(ZmV0,2)+(AmV-2*BmV-3*power(BmV,2)).*ZmV0-(AmV.*BmV-power(BmV,2)-power(BmV,3));
%公式26
F1_ZmV=3*power(ZmV0,2)-2*(1-BmV).*ZmV0+(AmV-2*BmV-3*power(BmV,2)).*ZmV0;
%公式27
ZmV=ZmV0-F_ZmV./F1_ZmV;
%公式28
dif1=abs((ZmV-ZmV0)/ZmV);
if (dif1<0.0005)
break
else
ZmV0=ZmV;
end
end
%ZmL
%公式29
AmL=amL.*P./(82.06*82.06*T.*T);
%公式30
BmL=bmL.*P./(82.06*T);
ZmL0=1;
while(1)
%公式31
F_ZmL=power(ZmL0,3)-(1-BmL).*power(ZmL0,2)+(AmL-2*BmL-3*power(BmL,2)).*ZmL0-(AmL.*BmL-power(BmL,2)-power(BmL,3));
%公式32
F1_ZmL=3*power(ZmL0,2)-2*(1-BmL).*ZmL0+(AmL-2*BmL-3*power(BmL,2)).*ZmL0;
%公式33
ZmL=ZmL0-F_ZmL./F1_ZmL;
%公式34
dif2=abs((ZmL-ZmL0)/ZmL);
if (dif2<0.0005)
break
else
ZmL0=ZmL;
end
end
%fiV
%公式36
HiV=0;
for j=1:8
HiV=HiV+YV(j).*(1-k(j)).*power(a.*ai.*a(j).*ai(j),0.5);
end
fiV=YV.*P.*exp(b./bmV.*(ZmV-1)-log(ZmV-BmV)-AmV./(2.*sqrt(2).*BmV).*(2.*HiV./amV-b./bmV).*(log((ZmV+2.414.*BmV)./(ZmV-0.414.*BmV))));
HiL=0;
for j=1:8
HiL=HiL+XL(j).*(1-k(j)).*power(a.*ai.*a(j).*ai(j),0.5);
end
fiL=XL.*P.*exp(b./bmL.*(ZmL-1)-log(ZmL-BmL)-AmL./(2.*sqrt(2).*BmL).*(2.*HiL./amL-b./bmL).*(log((ZmL+2.414.*BmL)./(ZmL-0.414.*BmL))));
dif2=abs(fiV-fiL);
if (dif2<0.0001)
break;
else
K=YV.*fiL./(XL.*fiV);
end
end
end