clc;clear;close all
%% Datos
P0=134000; % [Pa] Presión a la entrada del reactor
T0=431.1; % [K] Temperatura a la entrada del reactor
Ta=683.15; % [K] Temperatura de servicio
e=0.45; % Porosidad del catalizador
du=0.021; % [m] Diametro interno del tubo del reactor
L=3.5; % [m] Longitud del tubo del reactor
Ap=9.33e-5; % [m^2] Área de la superficie externa del catalizador
Ac=(du^2*pi)/4; % [m^2] Sección transversal del tubo del reactor
W=0.966; % [kg] masa del catalizador
rhob=W/(Ac* L); % [Kg/m^3] Densidad del lecho catalitico
Vp=W/rhob; % [m^3] Volumen del catalizador
rhoc=rhob/(1-e);% [Kg/m^3] Densidad del catalizador
Dp=6*Vp/Ap; % [m^3/m^2] Diametro de particula
m=(20.98*58.12+262.1*32+0.4*44.01+6.72*18); % [kg/h] Flujo másico
FT0=((20.98+262.1+0.4+6.72+985.995)); %[Kmol/h] Flujo total inicial a la entrada
del reactor
R=8.314; % [J/(mol*K)] Constante de los gases
Cto=(P0/(R*T0)/1000); % [kmol/m3] Concentración total inicial
U=107; % [W/(m2*K)] Coeficiente general de transferencia de calor
a=4/du;% 0.26923; % [m^-1]
G=(m/Ac); % [Kg/(m^2*s)] Velocidad de masa superficial
Ci=[20.980/11600 262.100/11600 0/11600 0.400/11600 6.720/11600 985.995/11600 P0
T0]; % Condiciones iniciales [Fa0 Fb0 Fc0 Fd0 Fe0 Fi0 P0 T0] [kmol/h];
L1=linspace(0,3.5);
[~,Flujos]=ode45(@Balances,L1,Ci,[],Cto,P0,T0,U,a,Ta,rhob,G,FT0,Dp,Ac,rhoc);
Xa=1-(Ci(1,1)./Flujos(:,1));
function [F]=Balances(~,Ci,Cto,P0,T0,U,a,Ta,rhob,G,FT0,Dp,Ac,rhoc)
Fa=Ci(1); Fb=Ci(2); Fc=Ci(3); Fd=Ci(4); Fe=Ci(5); Fi=Ci(6); P=Ci(7); T=Ci(8);
Ft=(Fa+Fb+Fc+Fd+Fe+Fi)/11600; % [kmol/h] Flujo total
%% Velocidades de rxn
r1=((2.191e-4*2616*(Cto*(Fa/Ft)*(P/P0)*(T0/T))*(Cto*(Fb/Ft)*(P/P0)*(T0/T))^...
0.2298)/(1+2616*(Cto*(Fa/Ft)*(P/P0)*(T0/T)))); % Reacción 1
r2=(7.028e-5*(Cto*(Fb/Ft)*(P/P0)*(T0/T))^0.2298); % Reacción 2
r3=(4.989e-6*(Cto*(Fc/Ft)*(P/P0)*(T0/T))*((Cto*(Fb/Ft)*(P/P0)*(T0/T))^0.6345...
/(Cto*(Fa/Ft)*(P/P0)*(T0/T))^1.151)); % Reacción 3
%% Balances molares
% Fa
dFa_dz=(rhoc*(1-0.45)*Ac)*(-r1-1/4*r2);
% Fb
dFb_dz=(rhoc*(1-0.45)*Ac)*(-3.5*r1-1.625*r2+3*r3);
%Fc
dFc_dz=(rhoc*(1-0.45)*Ac)*(r1+r3);
%Fd
dFd_dz=(rhoc*(1-0.45)*Ac)*(r2-4*r3);
%Fe
dFe_dz=(rhoc*(1-0.45)*Ac)*(4*r1+1/5*r2-r3);
% F inerte
DFi_dz=0;
%% Balance de energía
% Cp de cada componente
Cp_c=Cpc(T);
% Entalpia de cada rxn
DH_rxn=DHrxn(T);
% Ecuación diferencial
dT_dz=(rhoc*(1-0.45)*Ac)*((U*a*(Ta-T)
+r1*DH_rxn(1,1)+r2*DH_rxn(1,2)+r3*DH_rxn(1,3))...
/(rhob*(Fa*Cp_c(1,1)+Fb*Cp_c(1,2)*Fc*Cp_c(1,3)+Fd*Cp_c(1,4)+Fe*Cp_c(1,5))));
%% Caída de presión
% Fraciones molares
ya=Fa/Ft; yb=Fb/Ft; yc=Fc/Ft; yd=Fd/Ft; ye=Fe/Ft;
% Masa molar de la mezcla
Mmix=ya*58.12+yb*32+yc*98.06+yd*44.01+ye*18; % [kg/kmol]
% Densidad del lecho catalitico
rhoo=(Mmix/22.4*273.15/T0*P0/101325); %[Kg/m^3]
% viscosidad de la mezcla
u=viscosidad(Fa,Fi,Fd,Fe,Fb,P,T);
% Alfa
alfa=(2*G*(1-0.45))/(rhoo*Dp*0.45^3*rhob*Ac*P0)*(((150*(1-0.45)*u)/Dp+1.75*G));
% Ecuación diferencial de la caida de presión
dP_dz=(rhoc*(1-0.45)*Ac)*(-alfa/2*T/T0*P0/(P/P0)*Ft/FT0);
F=[dFa_dz ;dFb_dz; dFc_dz; dFd_dz; dFe_dz; DFi_dz; dP_dz;dT_dz ];
%disp(F)
end
%% Cp de cada compuesto
function Cp_c=Cpc(T)
% [A B C D E] Orden de los componentes
a=[9.487 28.106 -13.075 19.975 32.243];
b=[0.3313 -3.68e-6 0.3484 7.343e-2 1.923e-3];
c=[-1.108e-4 1.475e-5 -2.184e-4 -5.601e-5 1.055e-5];
d=[-2.2821e-9 -1.065e-8 4.839e-8 1.715e-8 -3.596e-9];
%Cp=a + b*T + c*T.^2 + d*T.^3
Cp_c=zeros(1,5);
for i=1:5
Cp=@(T) (a(1,i) + (b(1,i))*T + (c(1,i))*T.^2 + (d(1,i))*T.^3);
Cp_c(i)=Cp(T);
end
end
%% Entalpia de cada rxn
function DH_rxn_c=DHrxn(T)
% Hrxn= a+b*T+c*T^2+d*T^3+e*T^4
%[ rxn_1 rxn_2 rxn_3]
a=[-1242655.6 -2646628.4 -1428408];
b=[8.039 48.939 40.9];
c=[0.0124025 0.013971 -0.026385];
d=[-0.00004 -0.0000517 -0.0000119333];
e=[1.85255e-8 3.0700e-8 1.2141e-8];
DH_rxn_c=zeros(1,3);
for i=1:3
DH_rxn=@(T) (a(1,i)+b(1,i)*T+c(1,i)*T.^2+d(1,i)*T.^3+e(1,i)*T.^4);
DH_rxn_c(i)=DH_rxn(T);
end
end
% Viscosidad de la mezcla
function u=viscosidad(Fa,Fi,Fd,Fe,Fb,P,T)
% [n-butano N CO2 H2O O]
c1=[261.5278 17.8768 28.3005 730.3278 -5.3873];
c2=[12.5519 4.2383 5.5417 11.1122 2.6834];
c3=[-0.55997 0.0169 -0.0036 -1.1165 0.8772];
ui=zeros(1,5);
for i=1:5
ui(i)=(c1(1,i)+c2(1,i)*P/100+c3(1,i)*(T-273.15))*3600/1e6;
end
Ft=Fa+Fi+Fd+Fe+Fb;
u=(ui(1)*(Fa/Ft)+ui(2)*(Fi/Ft)+ui(3)*(Fd/Ft)+ui(4)*(Fe/Ft)+ui(5)*(Fb/
Ft))*(58.12+28+44.01+18+32)...
/(((Fa+Fi+Fd+Fe+Fb)/Ft*(58.12+28+44.01+18+32)));
end
