模拟心得
MATERIAL STUDIO 中SORPTION
第一个课题是模拟金属有机框架和共价有机框架吸附CO2以及分离CO2/CH4,使用的软件是Material studio，使用的是Sorption模块，输入的是逸度。
单组份求逸度的MATLAB程序，只需要在主程序窗 口输入function [rho,f] =PengRobinson(P1,T,N)（P1,T,N是具体的数值）就可以得到不同的压力和温度下的逸度。
function [rho,f] =PengRobinson(P1,T,N)
%+++++++++++++++++++++++++++++++++++++++++++++
%PengRobinson is ud to calculate the density and fugacity of single
%component gas at given pressure with Peng-Robinson equation.
%PengRobinson v1.00 beta include the parameter of n-alkanes(1-5), CO2(6)
%and CO(7).
%Where P1 means input pressure(kPa), T is temperature(K), N means the rial number of gas. rho
%is density, f is fugacity.
%e.g. If you wanna calculate density and fugacity of methane at 300kPa, 298k,you
%need input [rho,f] =PengRobinson(300,298,1).
%+++++++++++++++++++++++++++++++++++++++++++++
%t physical parameters
%+++++++++++++++++++++++++++++++++++++++++++++
P=P1./100;
t_M=[16.043 30.070 44.097 58.123 72.150 44.01 28.01];
t_omiga=[0.012 0.100 0.152 0.2 0.252 0.224 0.048];
t_Tc=[190.6 305.3 369.8 425.1 469.7 304.2 132.9 ];
t_Pc=[45.99 48.72 42.48 37.96 33.70 73.83 34.99];
%+++++++++++++++++++++++++++++++++++++++++++++
Tc=t_Tc(N);
Pc=t_Pc(N);
omiga=t_omiga(N);
M=t_M(N);
%+++++++++++++++++++++++++++++++++++++++++++++
R=83.14;
epsilon=1-2.^(0.5);
sigma=1+2.^(0.5);
%+++++++++++++++++++++++++++++++++++++++++++++
%calculate the Z of PR equation
%+++++++++++++++++++++++++++++++++++++++++++++
alpha=(1+(0.37464+1.54226*omiga-0.26992*omiga.^2)*(1-(T/Tc)^(0.5))).^2;
a=((0.45724*R^2*Tc^2)/Pc)*alpha;
b=0.07779.*R.*Tc./Pc;
beta=b.*P./(R.*T);
q=a./(b.*R.*T);
Z0=zeros(length(P),1);
Z1=ones(length(P),1);
for k=1:length(P)
while abs(Z1(k)-Z0(k))>1e-6
Z0(k)=Z1(k);
Z1(k)=1+beta(k)-q.*beta(k).*(Z0(k)-beta(k))./((Z0(k)+epsilon.*beta(k)).*(Z0(k)+sigma.*beta(k)));
end
end
I=(1./(sigma-epsilon)).*log((Z1+sigma.*beta)./(Z1+epsilon.*beta));
%+++++++++++++++++++++++++++++++++++++++++++++
%gain the density of gas
%+++++++++++++++++++++++++++++++++++++++++++++
rho=(P./(Z1.*R.*T)).*M.*1e6;
rho=vpa(rho,6);
phi=exp(Z1-1-log(Z1-beta)-q.*I);
f=phi.*P1;
f=vpa(f,5);
双组份的求逸度的程序：求逸度的过程和单组份的一样。双组份的逸度求解程序：
function [rho1,rho2,f1,f2] =PengRobinson_Binary(P1,T,N,y)
%+++++++++++++++++++++++++++++++++++++++++++++
%PengRobinson is ud to calculate the density and fugacity of binary
%component gas at given pressure with Peng-Robinson equation.
%PengRobinson v1.00 beta include the parameter of n-alkanes(1-5),
%isoButane(6) isoPentane(7), neoPentane(8) hydrogen(9) carbon dioxide(10)
%Where P1 means input pressure(kPa), T is temperature(K), N means the rial number of gas. rho
%is density(g/m3), f is fugacity.
%e.g. If you wanna calculate density and fugacity of mixture of methane and ethane at 300kPa,298k, you
%need input [rho,f] =PengRobinson(300,298,[1,2],[0.5,0.5]).
%+++++++++++++++++++++++++++++++++++++++++++++
%t physical parameters
%+++++++++++++++++++++++++++++++++++++++++++++
P=P1./100;
t_M=[16.043 30.070 44.097 58.123 72.150 58.123 72.150 72.150 2.016 44.01];
t_omiga=[0.012 0.100 0.152 0.2 0.252 0.181 0.229 0.197 -0.216 0.224];
t_Tc=[190.6 305.3 369.8 425.1 469.7 408.1 460.39 433.75 33.19 304.2];
t_Pc=[45.99 48.72 42.48 37.96 33.70 36.48 33.81 31.99 13.13 73.83];
%+++++++++++++++++++++++++++++++++++++++++++++
Tc=[t_Tc(N(1));t_Tc(N(2))];
Pc=[t_Pc(N(1));t_Pc(N(2))];
omiga=[t_omiga(N(1));t_omiga(N(2))];
M=[t_M(N(1));t_M(N(2))];
%+++++++++++++++++++++++++++++++++++++++++++++
R=83.14;
epsilon=1-2.^(0.5);
sigma=1+2.^(0.5);
%+++++++++++++++++++++++++++++++++++++++++++++
%calculate the Z of PR equation
%+++++++++++++++++++++++++++++++++++++++++++++
alpha=(1+(0.37464+1.54226.*omiga-0.26992.*omiga.^2).*(1-(T./Tc).^(0.5))).^2;
a=((0.45724.*R.^2.*Tc.^2)./Pc).*alpha;
b=0.07779.*R.*Tc./Pc;
a12=(a(1).*a(2)).^0.5;
am=(y(1).^2).*a(1)+2.*y(1).*y(2).*a12+(y(2).^2)*a(2);
bm=y(1).*b(1)+y(2).*b(2);
beta=bm.*P./(R.*T);
q=am./(bm.*R.*T);
Z0=zeros(length(P),1);
Z1=ones(length(P),1);
for k=1:length(P)
while abs(Z1(k)-Z0(k))>1e-6
Z0(k)=Z1(k);
Z1(k)=1+beta(k)-q.*beta(k).*(Z0(k)-beta(k))./((Z0(k)+epsilon.*beta(k)).*(Z0(k)+sigma.*beta(k)));
end
end
I=(1./(sigma-epsilon)).*log((Z1+sigma.*beta)./(Z1+epsilon.*beta));
%+++++++++++++++++++++++++++++++++++++++++++++
