Optimization of Plug Flow Reactor by Varying Temperature and Pressure

Asked by Sean Powers

Sean Powers (view profile)

on 16 Apr 2019
Latest activity Commented on by Bosnian Kingdom

Bosnian Kingdom (view profile)

on 24 Apr 2019
Accepted Answer by Torsten

Torsten (view profile)

I'm trying to optimize a complex reaction system to have the largest output of desired product C by varying the temperature and pressure of the reaction in an adiabatic PFR. I've created the design equation of the reaction:
Function dYdV = Designequation(v,y,T,P)
Fa = y(1);
Fb = y(2);
Fc = y(3);
Fd = y(4);
% Explicit equations
T = T(1);
P = P(1);
Cto = P / (8.314 * 10^-5) / T;
E1 = 15000;
E2 = 17500;
Ft = Fa + Fb + Fc + Fd;
Ca = Cto * Fa / Ft;
Cb = Cto * Fb/Ft;
Cc = Cto * Fc/Ft;
Cd = Cto * Fd/Ft;
k1 = 0.075 * exp(E1 /1.987 * (1/300 - (1/T)));
k2 = 0.0015 * exp(E2 / 1.987 * (1 / 300 - (1 / T)));
Fao = 50;
ra = -k1*Ca*Cb;
rb = -k2*Cb*Cc;
% Differential equations
dFbdV = (2*ra)+rb;
dFcdV = rb-ra;
dFddV = -rb;
dYdV = [dFadV; dFbdV; dFcdV; dFddV];
The design equation is solved for volume and the initial molar flow rates of A and B using ode113 with an initial temperature of 300K and 1 atm:
Vspan = [0 2]; % Range for the volume of the reactor
y0 = [50; 25; 0; 0]; % Initial values for the dependent variables i.e Fa,Fb,and Fc, Fd
T = 300;
P = 1;
[v, y] =ode113(@(v,y) Designequation(v,y,T,P),Vspan,y0);
dFbdT =y(:,2);
dFcdT =y(:,3);
dFddT = y(:,4);
plot(v,dFcdT); hold on
Now I want to find the largest value of the molar flow rate of C (dFcdT), by varying the temperature and pressure, but I haven't found a way that I can implement that based on searches of the MatLab forum.
I apologize if this has been answered in other forum posts, but I'm desperate.

Torsten

Torsten (view profile)

on 16 Apr 2019
So you want to maximize dFcdT(end) ? Or integral_{v=0}^{v=2} dFcdT dv ? Or something else ?
Sean Powers

Sean Powers (view profile)

on 16 Apr 2019
I want to maximice dFcdT(end) while keeping the vspan 0:2 and y0 the same.

Torsten (view profile)

on 16 Apr 2019
Edited by Torsten

Torsten (view profile)

on 16 Apr 2019

Although I guess that there are bounds on P and T, here is a code that will give you a start.
In case you want to constrain P and T, use "fmincon" instead of "fminsearch".
function main
format long
T = 300;
P = 1;
x0 = [T P];
options = optimset('MaxFunEvals',100000,'MaxIter',100000)
x = fminsearch(@fun,x0,options)
end
function obj = fun(x)
T = x(1);
P = x(2);
Vspan = [0 2]; % Range for the volume of the reactor
y0 = [50; 25; 0; 0]; % Initial values for the dependent variables i.e Fa,Fb,and Fc, Fd
[v, y] =ode15s(@(v,y) Designequation(v,y,T,P),Vspan,y0);
obj = -y(end,3)
end
function dYdV = Designequation(v,y,T,P)
Fa = y(1);
Fb = y(2);
Fc = y(3);
Fd = y(4);
% Explicit equations
Cto = P / (8.314 * 10^-5) / T;
E1 = 15000;
E2 = 17500;
Ft = Fa + Fb + Fc + Fd;
Ca = Cto * Fa / Ft;
Cb = Cto * Fb/Ft;
Cc = Cto * Fc/Ft;
Cd = Cto * Fd/Ft;
k1 = 0.075 * exp(E1 /1.987 * (1/300 - (1/T)));
k2 = 0.0015 * exp(E2 / 1.987 * (1 / 300 - (1 / T)));
Fao = 50;
ra = -k1*Ca*Cb;
rb = -k2*Cb*Cc;
% Differential equations
dFbdV = (2*ra)+rb;
dFcdV = rb-ra;
dFddV = -rb;
dYdV = [dFadV; dFbdV; dFcdV; dFddV];
end

Sean Powers

Sean Powers (view profile)

on 22 Apr 2019
Modify the function Designequation given your differential equations and reaction rates and then run iterate given a temperature and pressure range and the amount of steps you want to take in those ranges. i.e.
function dYdV = Designequation(v,y,T,P)
Fa = y(1);
Fb = y(2);
Fc = y(3);
Fd = y(4);
% Explicit equations
T = T(1);
P = P(1);
Cto = P / (8.314 * 10^-5) / T;
E1 = 15000;
E2 = 17500;
Ft = Fa + Fb + Fc + Fd;
Ca = Cto * Fa / Ft;
Cb = Cto * Fb/Ft;
Cc = Cto * Fc/Ft;
Cd = Cto * Fd/Ft;
k1 = 0.075 * exp(E1 /1.987 * (1/300 - (1/T)));
k2 = 0.0015 * exp(E2 / 1.987 * (1 / 300 - (1 / T)));
Fao = 50;
rA = -(rA1 + rA2)
rB = -(rB1 + rB2 + rB3)
rC = rC1 + rC3
rD = rD1 + rD2 + rD3
rE = rE2 + rE3
% Differential equations
d(FA)/d(V) = rA
d(FB)/d(V) = rB
d(FC)/d(V) = rC
d(FD)/d(V) = rD
d(FE)/d(V) = rE
d(FF)/d(V) = rF
dYdV = [d(FC)/d(V)];
end
I didnt fully modify the script for your problem, but its something like that
Torsten

Torsten (view profile)

on 23 Apr 2019
@Sean Powers:
And did the exaustive search lead to a different result than the usual optimization with "fminsearch" ?
Bosnian Kingdom

Bosnian Kingdom (view profile)

on 24 Apr 2019
Thank you, I'll try.