srk.m

化工热力学计算

function fai=SRK(T,P,V)
% SRK 方程
% 计算逸度系数
% 计算液相逸度系数时，z为液相摩尔分数；z为nx1 列向量
% 计算气相逸度系数时，z为气相摩尔分数；z为nx1 列向量% P=R*T/(Vm-b)-a/(Vm*(Vm+b));
R=8.3145;%R 摩尔气体常数，J/(mol*K)

% 顺序依次是 二甲醚DME,甲醇CH3OH,一氧化碳CO,二氧化碳CO2,氢气H2,水H2O
Tc=[400 512.6 132.9 304.2 33.2 647.3]'; %Tc是临界温度 K
Pc=[5.37 8.096 3.496 7.376 1.297 22.048]';%Pc是临界压力 Mpa
w=[0.192 0.559 0.049 0.225 -0.22 0.344]';%w是偏心因子
z=[0.0008 0.003 0.4355 0.0011 0.3556 0.0012]';
bi=0.086664*R*Tc./Pc/1e6;
b=sum(z.*bi);
m=0.48+1.574.*w-0.176*w.^2;
ai=(0.42747*R^2*Tc.^2./Pc/1e6).*(ones(size(m)).+m.*(ones(size(m)).-(T./Tc).^0.5)).^2;
ai(5)=0.42747*R^2*Tc(5)^2/Pc(5)/1e6*1.202*exp(-0.30288*T/Tc(5));  %修正氢气的系数
a=sum(sum((ai*ai').^0.5.*(z*z')));%k是二元交互作用系数，nxn矩阵
Z=P*1e6*V/R/T;%Z是压缩因子
A=a*P*1e6/(R*T)^2;
B=b*P*1e6/(R*T);
fai=exp(bi./b*(Z-1)-log(Z-B)-A/B*(2*ai.^0.5/a^0.5-bi/b)*log(1+B/Z));%fai是逸度系数，nx1列向量