En este video se resuelve el problema P11-6 del libro H. Scott Fogler, Essentials of Chemical Reaction Engineering, Prentice Hall, 2011.
Si tienen dudas o requieren asesorías sobre algunos temas relacionados no duden en contactarme:
WhatsApp: (+52) 2461474018

CODIGO
NOMBRE ARCHIVO: PBR_NI
% PROGRAMA PARA RESOLVER EDO
% USANDO EL METODO DE EULER

clear
clc
format long

%CONDICIONES DE ALIMENTACI?N
T0 = 300; %K
Fa0 = 5; %mol/min
Ca0 = 2; %mol/dm3
Ci = 2*Ca0; %mol/dm3

%DATOS CINETICOS
k1 = 0.1; %min-1 a 300 K
Tr1 = 300; %K
Ke1 = 1000; %
Tr2 = 300; %K

E = 10000; %cal/mol
dHrxn = -20000; %cal/mol

%DATOS DEL CATALIZADOR
dc = 1.2; %kg/dm3
alphadc = 0.02; %dm-3

%CAPACIDADES CALORIFICAS
cpa = 160; %cal/mol-k
cpb = 160; %cal/mol-k
cpi = 18; %cal/mol-k

%CONSTANTE DE LOS GASES
R = 1.987; %cal/mol-K

%Relacion de concentraciones
tethaI = Ci/Ca0; %para el inerte

% CONDICIONES INICIALES
CONV = 0;
DP = 1;

% RANGO DE INTEGRACION
W0 = 0; % Valor inicial de la variable dependiente
Wf = 48; % Valor final de variable dependiente

dW = 1e-3; % Tama?o de paso en el tiempo
n = (Wf-W0)/dW; % Numero de iteraciones

% INICIO DE VARIABLES PARA GRAFICAR
W(1) = W0;
x1(1) = CONV;
x2(1) = DP;
x3(1) = Ke1/(Ke1+1);
x4(1) = T0;
x5(1) = dc*k1*Ca0;

for i = 1:n
%%SISTEMA DE ECUACIONES
%BALANCE DE ENERG?A (SISTEMA ADIABATICO)
T = T0+CONV*(-dHrxn)/(cpa+tethaI*cpi);

%ECC DE Arrhenius
k = k1*exp((E/R)*((1/Tr1)-(1/T)));

%ECC DE van't Hoff
Kc = Ke1*exp((dHrxn/R)*((1/Tr2)-(1/T)));

%CONVERSION EN EL EQUILIBRIO
Xe = Kc/(1+Kc);

%CINETICA Y ESTEQUIOMETRIA (COMINACION)
ra = k*Ca0*((1-CONV)-CONV/Kc)*(T0/T)*DP;

%BLOQUE 1 y ECC DE ERGUN
dX = ra/(dc*Fa0);
dY = -alphadc/(2*DP*dc)*(T/T0);

%INTEGRA LA VARIABLE DEPENDIENTE CON EL METODO DE EULER
CONV = INT_EULER(dX,CONV);
DP = INT_EULER(dY,DP);

%CAMBIO A VOLUMEN

% ACTUALIZACIONES DE VARIABLES PARA GRAFICAR
W(i+1) = W(i) + dW;
x1(i+1) = CONV;
x2(i+1) = DP;
x3(i+1) = Xe;
x4(i+1) = T;
x5(i+1) = ra;

end

V = W./dc;


% GRAFICA Volumen vs variables de interes
subplot(4,1,1),plot(V,x1,V,x3),title('X y Xe')
xlabel('V'),ylabel('X')
subplot(4,1,2),plot(V,x2),title('caida de presion')
xlabel('V'),ylabel('y')
subplot(4,1,3),plot(V,x4),title('Temperatura')
xlabel('V'),ylabel('T')
subplot(4,1,4),plot(V,x5),title('ra')
xlabel('V'),ylabel('ra')

FUNCIÓN PARA RESOLVER POR METODO DE EULER
NOMBRE ARCHIVO: INT_EULER

function x = INT_EULER(fx0,x0)

dt = 1e-3;
x = x0 + fx0*dt;

end