%{
  Lecture 10
  Multiple Reaction Networks
%}

% Reactions in series:  A-->B-->C
clear
close all
home
sprintf('Reactions in series, first-order A-->B-->C')

k1 = 5;
k2 = 1;

dCdt=@(t,C) [-k1*C(1);
    k1*C(1)-k2*C(2);
    k2*C(2)];

time_range=[0 5];
init=[1 0 0];

[t,C]=ode45(dCdt,time_range,init);
h=figure('color',[1 1 1]);
set(h,'WindowStyle','Docked')
subplot(2,1,1)
plot(t,C(:,1),t,C(:,2),t,C(:,3))
legend('A','B','C')
xlabel('time')
ylabel('concentration')
title('A \rightarrow B \rightarrow C, k_1 = 5k_2')

k1=1;
k2=5;
dCdt=@(t,C) [-k1*C(1);
    k1*C(1)-k2*C(2);
    k2*C(2)];

[t,C]=ode45(dCdt,time_range,init);
subplot(2,1,2)
plot(t,C(:,1),t,C(:,2),t,C(:,3))
legend('A','B','C')
xlabel('time')
ylabel('concentration')
title('A \rightarrow B \rightarrow C, 5k_1 = k_2')

sprintf('First-order reactions in parallel, A-->B, A-->C')
k1=5;
k2=1;
dCdt = @(t,C) [-(k1+k2)*C(1);
    k1*C(1);
    k2*C(1)];

[t,C]=ode45(dCdt,time_range,init);
h=figure('color',[1 1 1]);
set(h,'WindowStyle','Docked')
subplot(2,1,1)
plot(t,C(:,1),t,C(:,2),t,C(:,3))
legend('A','B','C')
xlabel('time')
ylabel('concentration')
title('A \rightarrow B, A \rightarrow C, k_1 = 5k_2')

k1=1;
k2=5;
dCdt = @(t,C) [-(k1+k2)*C(1);
    k1*C(1);
    k2*C(1)];

[t,C]=ode45(dCdt,time_range,init);
subplot(2,1,2)
plot(t,C(:,1),t,C(:,2),t,C(:,3))
legend('A','B','C')
xlabel('time')
ylabel('concentration')
title('A \rightarrow B, A \rightarrow C, 5k_1 = k_2')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sprintf('Competing first and second order Reactions')

k1=1;
k2=1;
t=[0:0.01:10];
for CA0=0.5:0.5:1.5
% CSTR
    fCA=@(t) (-(1+k2*t)+((1+k2*t).^2+4*CA0*k1*t).^(0.5))./(2*k1*t);
    CA=fCA(t);
    CB= k1*t.*CA.*CA;
    CC= k2*t.*CA;
    h=figure('color',[1 1 1]);
    set(h,'WindowStyle','Docked')
    subplot(2,1,1)
    plot(t,CA,t,CB,t,CC,t,CB./CC); title(['Competing first and second order reactions in a CSTR, C_{A0}=',num2str(CA0)])
%plot(t,CB,t,CB./CC)
%plot(t,CC,t,CC./CB)
    legend('C_A','C_B','C_C','C_B/C_C')
    xlabel('residence time')
    ylabel('concentration (M)')
    hold

% PFR
    subplot(2,1,2)
    dCdt=@(t,C) [-k1*C(1).*C(1)-k2*C(1);
    k1*C(1).*C(1);
    k2*C(1)];
    time_range=[0 10];
    init=[CA0;0;0];
    [t,C]=ode45(dCdt,time_range,init);
%plot(t,C(:,2),'bo',t,C(:,2)./C(:,3),'go')
    plot(t,C,t,C(:,2)./C(:,3))
    legend('C_A','C_B','C_C','C_B/C_C')
%plot(t,C(:,3),'b+',t,C(:,3)./C(:,2),'g+')
    title(['Competing first and second order reactions in a PFR, C_{A0}=',num2str(CA0)])
end
%title('CSTR vs. PFR for competing first and second order reactions (CSTR solid lines)')
%legend('CSTR C_B','CSTR C_B/C_C','PFR C_B','PFR C_B/C_C')
%{
 Benzene condensation

  2 B --> D + H
  B + D --> T + H

  PFR
%}
sprintf('Benzene condensation in terms of species rates')
nu = [ -2 1 1 0;
    -1 -1 1 1];
R0=0.0821; % l atm/mol K
P0 = 1; % atm
T0 = 1033; % K
FB0 = 60000; % mol/hr
CB0 = P0/(R0*T0); % mol/l
v0 = FB0/CB0;

k1 = 7e5;
k2 = 4e5;
K1 = 0.31;
K2 = 0.48;

Q1 = @(F) (F(1)^nu(1,1))*(F(2)^nu(1,2))*(F(3)^nu(1,3))*(F(4)^nu(1,4));
Q2 = @(F) (F(1)^nu(2,1))*(F(2)^nu(2,2))*(F(3)^nu(2,3))*(F(4)^nu(2,4));

R = @(F) [(k1/v0^2)*(F(1)^2-F(2)*F(3)/K1);
    (k2/v0^2)*(F(1)*F(2)-F(3)*F(4)/K2)];

dFdV = @(t,F) nu'*R(F);

volume_range=[0 1000];
init=[FB0 0 0 0];
[V,F]=ode45(dFdV,volume_range,init);

XB=(F(1,1)-F(:,1))/F(1,1);

h=figure('color',[1 1 1]);
set(h,'WindowStyle','Docked')
subplot(2,1,1)
plot(V,XB)
xlabel('Volume (l)')
ylabel('B conversion')
title('Benzene condensation reaction network, PFR')

subplot(2,1,2)
plot(V,F)
legend('B','D','H','T')
xlabel('Volume (l)')
ylabel('Flow rates (mol/hr)')

index=find(abs(XB-0.5)<0.001);

sprintf('Volume for %f%% conversion %f liters',XB(index)*100,V(index))

h=figure('color',[1 1 1]);
set(h,'WindowStyle','Docked')
SC=F(:,2)./F(:,4);
plot(V,SC,V(index),SC(index),'r>')
xlabel('Volume')
ylabel('Selectivity D/T')
title('Benzene Condensation Network, Selectivity Plot')
axis([0 1000 0 50])
ans =
Reactions in series, first-order A-->B-->C
ans =
First-order reactions in parallel, A-->B, A-->C
ans =
Competing first and second order Reactions
Warning: Divide by zero.
Current plot held
Warning: Divide by zero.
Warning: Divide by zero.
Current plot held
Warning: Divide by zero.
Warning: Divide by zero.
Current plot held
Warning: Divide by zero.
ans =
Benzene condensation in terms of species rates
ans =
Volume for 49.924484% conversion 401.440322 liters
Warning: Divide by zero.