||
单组分
function [rho,f] =PengRobinson(P1,T,N)
%+++++++++++++++++++++++++++++++++++++++++++++
%PengRobinson is used 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), O2(8).
%Where P1 means input pressure(kPa), T is temperature(K), N means the serial 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).
%+++++++++++++++++++++++++++++++++++++++++++++
%set physical parameters(分子量,偏心因子,临界温度K,临界压力bar)
%+++++++++++++++++++++++++++++++++++++++++++++
P=P1./100;
t_M=[16.043 30.070 44.097 58.123 72.150 44.01 28.01 31.999];
t_omiga=[0.012 0.100 0.152 0.2 0.252 0.224 0.048 0.022];
t_Tc=[190.6 305.3 369.8 425.1 469.7 304.2 132.9 154.58];
t_Pc=[45.99 48.72 42.48 37.96 33.70 73.83 34.99 50.43];
%+++++++++++++++++++++++++++++++++++++++++++++
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);
disp(f);
disp(rho);
双组份
function [rho1,rho2,f1,f2] =PengRobinson_Binary(P1,T,N,y)
%+++++++++++++++++++++++++++++++++++++++++++++
%PengRobinson is used 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) CO(11)
%Where P1 means input pressure(kPa), T is temperature(K), N means the serial 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]).
%+++++++++++++++++++++++++++++++++++++++++++++
%set 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 28.01];
t_omiga=[0.012 0.100 0.152 0.2 0.252 0.181 0.229 0.197 -0.216 0.224 0.048];
t_Tc=[190.6 305.3 369.8 425.1 469.7 408.1 460.39 433.75 33.19 304.2 132.9];
t_Pc=[45.99 48.72 42.48 37.96 33.70 36.48 33.81 31.99 13.13 73.83 34.99];
%+++++++++++++++++++++++++++++++++++++++++++++
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));
%+++++++++++++++++++++++++++++++++++++++++++++
%gain the density of gas
%+++++++++++++++++++++++++++++++++++++++++++++
%rho=(P./(Z1.*R.*T)).*M.*1e6;
%rho=vpa(rho,6);
rho=(P./(Z1.*R.*T)).*1e6;
rho1=vpa(y(1).*rho.*M(1),6);
rho2=vpa(y(2).*rho.*M(2),6);
phi1=zeros(length(P),1);
phi2=zeros(length(P),1);
f1=zeros(length(P),1);
f2=zeros(length(P),1);
phi1=exp((b(1)./bm).*(Z1-1)-log(Z1-beta)-q.*I.*(2.*(y(1).*a(1)+y(2).*a12)./am-b(1)./bm));
phi2=exp((b(2)./bm).*(Z1-1)-log(Z1-beta)-q.*I.*(2.*(y(2).*a(2)+y(1).*a12)./am-b(2)./bm));
f1=phi1.*P1.*y(1);
f2=phi2.*P1.*y(2);
f1=vpa(f1,5);
f2=vpa(f2,5);
disp(rho);
disp(f1);
disp(f2);
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-27 15:19
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社