单组分

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);