clear all
global X;
global CA0;
global i;
global k;
global Ke;
global KA;
global KB;
global KC;
global KD;
global tetb;
global rop;
global npontos;
global eta;
% Constantes
%KA=0.875258325;%(dm3/mol)
%KB=1.286220417;%(dm3/mol)
%KC=0.039345976;%(dm3/mol)
%KD=0.075259805;%(dm3/mol)
k(1) = 0.02907; % k'ap do trabalho 1 (dm^6/g.mol.min)
KA = 0.875258325;
KB = 1.286220417;
KC = 0.039345976;
KD = 0.075259805;
CA0=4.850; % Concentração do reagente limitante à entrada (M)
rop=322; % Catalyst density (g/l)
v0=0.005; % Volumetric flow- Feed (L/min)
Dc=0.1; % Reactor cross section diameter (dm)
tetb=1.66804;% (CB0/CA0)
eb=0.45; % Bulk porosity
Cp=8; % Specific heat- Feed (cal/mol.K)
HR=-5000; % Het of reaction (cal/mol)
Ea=15000; % Activation energy (cal/mol)
R=1.987; % Perfect Gas constant (cal/mol.K)
T0=335; % Feed temperature (K)
Lmax=0.8; % Reactor length (dm). 1st we have to give a huge
value t achieve the new Xe!
% Then we change Lmax to chieve X= 0.995Xe
npontos = 500; % Total number of points
z(1)=0; % Linear variable z along the rector (dm)
k(1)=0.02907; % intrinsic kinetic constant at the entry- kap from
work#1
Ke(1)=3.59855; % equilíbrium constant at the entry of the reactor
Ac = pi*Dc^2/4 ; % Cross section area
dz = Lmax/npontos; % Increment of the independent variable;
X(1)=0; % Conversion at the entry
TR=T0;
T(1)=TR;
for i=1:npontos
z(i+1)=z(i)+dz;
k(i)=k(1)*exp((-Ea/R)*((1/T(i))-(1/TR))); % Arrhenius Law
Ke(i)=Ke(1)*exp((-HR/R)*((1/T(i))-(1/TR))); % VAn´t Hoff Law
% Integration of the equation of the molar balance to the pellet
esferabvp4exemplotc;
disp('z');
disp(z(i+1));
% Velocidade da reação observada
robs = eta*((k(i)*CA0*(1-X(i))*(((tetb-X(i))/(1-X(i)))-
(((X(i)/(1-X(i)))^2)/Ke(i))))/(1+(KA*CA0*(1-
X(i)))+((KB*CA0*(X(i)^2))/(Ke(i)*(1-X(i))))+((KC+KD)*CA0*(1-
X(i))*((X(i)/(1-X(i)))^2))));
% Molar balance to the reactor
X(i+1)=X(i)+(rop*(1-eb)*Ac/(CA0*v0))*robs*dz;
% Energy balance
T(i+1)=T(i)+(-HR/((1+tetb)*Cp))*((rop*(1-
eb)*Ac*robs)/(CA0*v0))*dz;
end
Tmax=T(i+1);
disp('Conversion at the end:');
disp('X=');
disp(X(i+1));
disp('Temperatura à saída:');
disp('T=');
disp(T(i+1));
disp('L (dm):');
disp(Lmax); % Comprimento do reator
%-----------------------------------------------------------------------
% Grafico da conversão
figure
plot(z, X);
axis([0 Lmax 0 1]);
xlabel('z - Comprimento do reator');
ylabel('X - Conversão');
% Gráfico da temperatura
figure;
plot(z,T);
axis([0 Lmax T0 Tmax]);
xlabel('z - Comprimento do reator');
ylabel('T - Temperatura');
%save profileXpb.txt X -ASCII;
%save profileXpb.txt T -ASCII;
