小弟初學(xué)氣液相平衡

根據(jù)已知的一組相平衡數(shù)據(jù)，學(xué)習(xí)用PR方程計(jì)算組分逸度，但計(jì)算結(jié)果總是不對(duì)!

現(xiàn)在有一個(gè)問(wèn)題迷惑，在用PR方程計(jì)算壓縮因子時(shí)，求得三個(gè)結(jié)果

不知哪個(gè)是氣態(tài)的壓縮因子，哪個(gè)是液態(tài)的壓縮因子？

先貼上此問(wèn)題的m文件，文件中有詳細(xì)說(shuō)明

望各位前輩指教一二，不勝感激！

clear
clc
format long
%在以下壓力、溫度以及組分條件下，氮氧兩組份滿足相平衡，相平衡常數(shù)KN=1.14，這段程序是為了學(xué)習(xí)并驗(yàn)證PR方程

p=5.23e5;    %壓力
T=96.6;    %溫度
yN=0.904;  %氣相氮組分
yO=1-yN;    %氣態(tài)氧組分
xN=0.793;     %液態(tài)氮組分
xO=1-xN;    %液態(tài)氧組分

%物性基本常數(shù)
%N2
TcN=126.2;  %臨界溫度
pcN=34e5;   %臨界壓力
wN=0.03726; %偏心因子
RN=296.8;   %氣體常數(shù)
mN=0.37464+1.54226*wN-0.26992*wN^2; %用來(lái)計(jì)算PR方程中的a(T)

%O2
TcO=154.58; %同氮?dú)?br /> pcO=50.43e5;
wO=0.022;
RO=259.84;
mO=0.37464+1.54226*wO-0.26992*wO^2;

%PR方程參數(shù)求取
%N2

TrN=T/TcN;  %對(duì)比溫度
alphaN=(1+mN*(1-TrN^0.5))^2;    %用來(lái)求a(T)
acN=0.45724*RN^2*TcN^2/pcN;     %用來(lái)求a(T)
aN=acN*alphaN;                  %a(T)計(jì)算式
bN=0.0778*RN*TcN/pcN;           %b計(jì)算式
AN=aN*p/(RN^2*T^2);             %A計(jì)算式
BN=bN*p/RN/T;                   %B計(jì)算式

aaN=BN-1;
bbN=AN-2*BN-3*BN^2;
ccN=BN^3+BN^2-AN*BN;            %以上三個(gè)參數(shù)用做求體積壓縮率的多項(xiàng)式系數(shù)

%O2  注釋同氮?dú)?br /> TrO=T/TcO;
alphaO=(1+mO*(1-TrO^0.5))^2;
acO=0.45724*RO^2*TcO^2/pcO;
aO=acO*alphaO;
bO=0.0778*RO*TcO/pcO;
AO=aO*p/(RO^2*T^2);
BO=bO*p/RO/T;

aaO=BO-1;
bbO=AO-2*BO-3*BO^2;
ccO=BO^3+BO^2-AO*BO;

%混合規(guī)則

%氣相
kNO=-0.023644+0.15427e-3;   %兩元交互因子
am=yN*yN*aN+yO*yO*aO+2*yN*yO*sqrt(aN*aO)*(1-kNO);   %混合a
% am=(yN*sqrt(aN)+yO*sqrt(aO))^2;
bm=yN*bN+yO*bO;                                     %混合物的b
Am=am*p/(RN^2*T^2);
Bm=bm*p/RN/T;

aam=Bm-1;
bbm=Am-2*Bm-3*Bm^2;
ccm=Bm^3+Bm^2-Am*Bm;     %以上三個(gè)參數(shù)用做求體積壓縮率的多項(xiàng)式系數(shù)


%液相
am1=xN*xN*aN+xO*xO*aO+2*xN*xO*sqrt(aN*aO)*(1-kNO);
bm1=xN*bN+(1-xN)*bO;
Am1=am1*p/(RN^2*T^2);
Bm1=bm1*p/RN/T;

aam1=Bm1-1;
bbm1=Am1-2*Bm1-3*Bm1^2;
ccm1=Bm1^3+Bm1^2-Am1*Bm1;



%求體積壓縮因子
n=[1,aam,bbm,ccm];
z=roots(n)
n1=[1,aam1,bbm1,ccm1];
z1=roots(n1)


%求兩相N2的逸度


phi_yN=exp(bN*(z-1)/bm-log(z-Bm)-Am/(2*sqrt(2))*((2*(yO*sqrt(aN*aO)+yN*sqrt(aN*aN)))/am-bN/bm)*(log((z+2.414*Bm)./(z+0.414*Bm))))
phi_xN=exp(bN*(z1-1)/bm1-log(z1-Bm1)-Am1/(2*sqrt(2))*((2*((1-xN)*sqrt(aN*aO)+xN*sqrt(aN*aN)))/am1-bN/bm1)*(log((z1+2.414*Bm1)./(z1+0.414*Bm1))))


%      無(wú)論氣相還是液相都取得了三個(gè)體積壓縮因子，根據(jù)相平衡的原則，KN=1.14應(yīng)該等于phi_xN/phi_yN，分別按照三個(gè)壓縮因子求得逸度，但無(wú)
%      論怎么組合，都不能滿足以上式子。
%問(wèn)題：可能是因?yàn)閴嚎s因子求取有問(wèn)題?另問(wèn)：當(dāng)求得三個(gè)壓縮因子時(shí)，到底哪個(gè)用作液相的壓縮因子，哪個(gè)用作氣相壓縮因子? 返回小木蟲查看更多