% Simulates volcano plot for A-to-B reaction
clear
clc

% Constants:
% Gas constant
R = 8.31446261815324; % J/mol-K
Rg = R*0.01; % L-bar/mol-K

% Conditions:
% Reactor temperature
Tc = 200; % deg C
T = Tc + 273.15; % K

% Feed pressure (bar)
Pf(1) = 100; % A(g)
Pf(2) = 0; % B(g)

% Feed concentration
Cf = Pf/Rg/T; % M

% Initial conditions
x0 = [Cf';0;0];

% Volumetric flowrate
qsdot = 5; % mL/min
qdot = qsdot/60000; % L/s

% Catalyst weight
w = 138; % mg
% Catalyst density
rho = 0.53; % g/mL
% Reactor bed volume (L)
V = w*0.001/rho*0.001;

% Space velocity
SV = qdot/V; % 1/s

% Parameters:
% Vapor species enthalpy (kJ/mol)
H(1) = 0; % A(g)
H(2) = 0; % B(g)

% Linear scaling relationship
gamma = -2.0; % unitless
delta = -1.4; % eV

% Bronsted-Evans-Polanyi relationship
alpha = 0.6; % unitless
beta = 102; % kJ/mol

% Initial binding energy
BEa0 = 1.0; % eV

% Binding energy shift bounds (eV)
UL = 0.3; % Weak binding
UR = 0.5; % Strong binding
int = 0.005; % Sampling interval

% Loop indices
is = UL/int;
ie = UR/int;

% Preallocate matrices
TOFb = zeros(ie - is + 1,1);
Theta = zeros(ie - is + 1,3);

% Time-on-stream
tos = 5e0; % s

% Target conversion
C = 1; % percent

% Iterate over binding energy
for i = 1:ie - is + 1
    
    % Binding energy shift (eV)
    delBEa = int*(i - 1 + is);
    
    % Compute rate constants
    k = cstr_constants(BEa0,delBEa,gamma,H,delta,alpha,beta,T);
    
    % Compute CSTR performance
    [TOFb(i),Theta(i,:)] = volcano(k,T,SV,Cf,tos,x0,qdot,C);
end

% Compute binding energy matrix
BEa(:,1) = BEa0 + int*(is:ie); % eV

% Collect results
results = [BEa,TOFb,Theta];

% Plot Sabatier volcano
semilogy(BEa,TOFb)
xlabel('Binding energy of A (eV)')
ylabel('Production rate of B (1/s)')

% Plot surface coverage
figure(2)
plot(BEa,Theta)
xlabel('Binding energy of A (eV)')
ylabel('Surface coverage (-)')
legend('Theta A*','Theta B*','Theta *')

%% Functions
function k = cstr_constants(BEa0,delBEa,gamma,H,delta,alpha,beta,T)

% Constants:
% Gas constant
R = 8.31446261815324; % J/mol-K
R = R*0.001; % kJ/mol-K

% Boltzmann constant
kB = 1.380649e-23; % J/K

% Planck constant
h = 6.62607015e-34; % J-s

% Parameters:
% Binding energy (eV)
BE(1) = BEa0 + delBEa; % A*
BE(2) = gamma*BE(1) + H(2)/96.485 - gamma*H(1)/96.485 - (1 - gamma)*delta; % B*
% Restrict to positive values
BE = max(BE,0);

% Heat of reaction (kJ/mol)
delH(1) = -BE(1)*96.485; % A(g) + * <--> A*
delH(2) = (H(2) - BE(2)*96.485) - (H(1) - BE(1)*96.485); % A* <--> B*
delH(3) = BE(2)*96.485; % B* <--> B(g) + *

% Activation energy (kJ/mol)
Ea = alpha*delH(2) + beta;
% Restrict to positive values
Ea = max(Ea,delH(2));
Ea = max(Ea,0);

% Entropy of reaction (J/mol-K)
delS(1) = -135; % A(g) + * <--> A*
delS(2) = 0; % A* <--> B*
delS(3) = 135; % B* <--> B(g) + *

% Gibbs free energy of reaction 
delG = delH - T*delS*0.001; % 

% Equilibrium constants
K = exp(-delG/R/T); % unitless

% Compute rate constants:
% Forward (1/s)
k(1) = kB*T/h*exp(delS(1)*0.001/R); % A(g) + * --> A*
k(3) = kB*T/h*exp(-Ea/R/T); % A* --> B*
k(6) = kB*T/h*exp(-delS(3)*0.001/R); % B(g) + * --> B*

% Reverse (1/s)
k(2:2:4) = k(1:2:3)./K(1:2); 
k(5) = k(6)*K(3); % B* --> B(g) + *
end

function xdot = cstr(t,x,k,T,SV,Cf)

% Constants:
% Gas constant
R = 8.31446261815324; % J/mol-K
Rg = R*0.01; % L-bar/mol-K

% Reference pressure
Po = 1.01325; % bar

% Parameters:
% Catalyst weight
w = 138; % mg
% Active site density
srho = 20; % micromol/g
% Number of active sites (gmol)
Nsites = w*0.001*srho*0.000001;

% Catalyst density
rho = 0.53; % g/mL
% Reactor bed volume (L)
V = w*0.001/rho*0.001;

% Stoichiometry matrix (unitless)
nu(1,:) = [-1 0 0]; % A(g)
nu(2,:) = [0 0 1]; % B(g)
nu(3,:) = [1 -1 0]; % A*
nu(4,:) = [0 1 -1]; % B*

% Site balance (unitless)
theta_star = 1 - sum(x(3:4));

% Rate expressions (1/s)
r(1,1) = k(1)*x(1)*Rg*T/Po*theta_star - k(2)*x(3); % A(g) + * <--> A*
r(2,1) = k(3)*x(3) - k(4)*x(4); % A* <--> B*
r(3,1) = k(5)*x(4) - k(6)*x(2)*Rg*T/Po*theta_star; % B* <--> B(g) + *

% Compute xdot:
% Vapor species (M/s)
xdot(1:2,1) = nu(1:2,:)*r*Nsites/V + SV*(Cf' - x(1:2));

% Surface species (1/s)
xdot(3:4,1) = nu(3:4,:)*r;
end

function [TOFb,Theta] = volcano(k,T,SV,Cf,tos,x0,qdot,C)

% ODE solver options
options = odeset('RelTol',1e-9,'AbsTol',1e-12);

% Iterate until convergence
for n = 1:9e18
    
    % Obtain CSTR data
    [~,x] = ode15s(@(t,x) cstr(t,x,k,T,SV,Cf),[0 tos],x0,options);

    % Data analysis:
    % Catalyst weight
    w = 138; % mg
    % Active site density
    srho = 20; % micromol/g
    % Number of active sites (gmol)
    Nsites = w*0.001*srho*0.000001;

    % Compute conversion (%):
    Xa = (Cf(1) - x(end,1))/Cf(1)*100;

    % Compute consumption rate (1/s):
    TOFa = (Cf(1) - x(end,1))*qdot/Nsites;
    % Compute production rate (1/s):
    TOFb = (x(end,2) - Cf(2))*qdot/Nsites;

    % Compute mass balance
    MB = TOFb/TOFa*100; % mol %
   
    % Check steady state
    if abs(MB - 100) > 1
        % Modify time span
        tos = tos*10; % s
    else
        % Check conversion
        if abs(Xa - C) > 0.01
            % Modify parameters
            SV = SV*abs(Xa)/C; % 1/s
            qdot = qdot*abs(Xa)/C; % L/s
        else
            % Exit infinite loop
            break
        end
    end
end

% Retrieve surface coverage (unitless)
Theta(1) = x(end,3); % A*
Theta(2) = x(end,4); % B*
Theta(3) = 1 - sum(x(end,3:4)); % *
end