%{
  Lecture 13: Non-isothermal reactors
%}
clear;
home;

% adiabatic reaction, A -> B
R0 = 8.31441;
T0 = 300;

Ea = 125000;  % J/mol
k0 = 2.9e20;  % min-1
dH = -90000;  % J/mol
rhoCp = 4000; % J/l K
CA0 = 2;      % M

% T, k, and rate relations same for CSTR and PFR
T=@(X) T0-CA0*X*dH/rhoCp;        % linear
k=@(X) k0*exp(-Ea./(R0*T(X)));   % exponential
rate=@(X) CA0*(1-X).*k(X);

X=[0:0.01:0.99];
TX=T(X);
h=figure('color',[1 1 1]);
set(h,'WindowStyle','Docked')
plot(X,TX)
xlabel('Conversion')
ylabel('Temperature (K)')
title('Temperature and conversion linearly related in adiabatic reactor')

%%%%%%%%%%%%%%%%%%%%%%%%%%
% PFR
%%%%%%%%%%%%%%%%%%%%%%%%%%
% isothermal
tauX=-log(1-X)/k(0);
invrate=1./(k(0)*CA0*(1-X));

h=figure('color',[1 1 1]); set(h,'WindowStyle','Docked')
subplot(2,2,1)
plot(tauX,X)
xlabel('\tau (minutes)');ylabel('Conversion')
title('Isothermal PFR, A \rightarrow B')

subplot(2,2,2)
plot(invrate,X)
ylabel('Conversion'); xlabel('1/rate');
title('Isothermal PFR, A \rightarrow B')

% adiabatic
dXdt = @(t,X) rate(X)/CA0;
tau_range=[0 4];
init=0;
[tau,Xpfr] = ode45(dXdt,tau_range,init);

subplot(2,2,3)
plot(tau,Xpfr); axis([0 4 0 1]);
xlabel('\tau (minutes)'); ylabel('Conversion')
title('Adiabatic PFR, A \rightarrow B')

subplot(2,2,4)
plot(1./rate(X),X)
ylabel('Conversion'); xlabel('1/rate')
title('Adiabatic PFR, A \rightarrow B')


%%%%%%%%%%%%%%%%%%%%%%%%%%
% CSTR
%%%%%%%%%%%%%%%%%%%%%%%%%%
X=[0:0.01:0.99];
h=figure('color',[1 1 1]);
set(h,'WindowStyle','Docked')

% isothermal
tauX = (1/k(0))*X./(1-X);
subplot(2,2,1)
plot(tauX,X)
xlabel('\tau (minutes)');ylabel('Conversion')
title('Isothermal CSTR, A \rightarrow B')

subplot(2,2,2)
invrate=1./(k(0)*CA0*(1-X));
plot(invrate,X)
ylabel('Conversion');xlabel('1/rate')
title('Isothermal CSTR, A \rightarrow B')

% adiabatic
tau=@(X) (X./(1-X))./k(X);

tauX=tau(X);
rateX=rate(X);
invrate=1./rate(X);
kX=k(X);

subplot(2,2,3)
plot(tauX,X)
xlabel('\tau (minutes)');ylabel('Conversion')
title('Adiabatic CSTR, A \rightarrow B')

subplot(2,2,4)
plot(invrate,X)
ylabel('Conversion');xlabel('1/rate')
title('Adiabatic CSTR, A \rightarrow B')

%
% Decompose adiabatic rate
%
h=figure('color',[1 1 1]);
set(h,'WindowStyle','Docked')
subplot(3,1,1)
plot(X,invrate)
ylabel('1/rate (min liter/mol)')
subplot(3,1,2)
plot(X,1./kX)
ylabel('1/k')
subplot(3,1,3)
plot(X,1./(CA0*(1-X)))
ylabel('1/CA')
xlabel('Conversion')

%
%  PFR with heat exchange, non-adiabatic
%
R0 = 8.31441;
T0 = 300;

Ea = 125000;  % J/mol
k0 = 2.9e20;  % min-1
dH = -90000;  % J/mol
rhoCp = 4000; % J/l K
CA0 = 2;      % M
Ta = T0-5;

init=[0 T0];
time_range=[0 80];

Ua0 = 250; % J/K s
h=figure('color',[1 1 1]);
set(h,'WindowStyle','Docked')
hold all % prevents color from being reset to blue each time
for i=0:10
    Ua=Ua0*i;
    dAdt = @(t,A) [k0*exp(-Ea/(R0*A(2)))*(1-A(1));
    (Ua*(Ta-A(2))-k0*exp(-Ea/(R0*A(2)))*CA0*(1-A(1))*dH)/rhoCp];
    [tau,A]=ode45(dAdt,time_range,init);

    plot(A(:,2),A(:,1))
end
% A(:,1) is conversion, A(:,2) is temperature
xlabel('Temperature (K)'); ylabel('Conversion')
axis([295 350 0 1])
title({'PFR with various degrees of heat exchange';'T_0 = 300 K, T_a = 295 K'})

h=figure('color',[1 1 1]);
set(h,'WindowStyle','Docked')
Ua=750;
init=[0 T0];
time_range=[0 20];
Ta = T0-5;

dAdt = @(t,A) [k0*exp(-Ea/(R0*A(2)))*(1-A(1));
(Ua*(Ta-A(2))-k0*exp(-Ea/(R0*A(2)))*CA0*(1-A(1))*dH)/rhoCp];
[tau,A]=ode45(dAdt,time_range,init);
subplot(2,2,1)
plot(tau,A(:,1))
xlabel('\tau');ylabel('Conversion');
title({'A \rightarrow B, PFR with heat exchange';'Ua=750 J/K s'})

subplot(2,2,3)
plot(tau,A(:,2))
xlabel('\tau');ylabel('Temperature');

subplot(2,2,2)
rate=k0*CA0*exp(-(Ea/R0)./A(:,2)).*(1-A(:,1));
invrate=1./rate;
plot(invrate,A(:,1))
axis([0 20 0 1])
xlabel('1/rate');ylabel('Conversion')
title('Levenspiel on his side')

Published with MATLAB® 7.0.1