CSTR Reaction Model: Yield vs Conversion graph
5 views (last 30 days)
Show older comments
I am currently modelling a CSTR in which a third order reaction occurs to produce C and second order for D. A parallel reaction occurs. A + B ---> C and A + B ---> D. I'm trying to produce a yield profile graph that looks something like this by testing different values for k2. The code that I used is shown below. I used this approach for a fed-batch and pfr and produced succesful results. My CSTR doesn't work though. Any ideas how to fix it?

clear all; close all
clc
%CONTINUOUS STIRRED TANK REACTOR CODE CONCENTRATION AS A FUNCTION OF CSTR
%TIME
%time span specification
timeStart = 0; timeEnd = 600;%s
timeSpan = [timeStart:0.1:timeEnd];
cao = 0.13;
%initial condition specification
initCA = 0.13; initCB = 0.1; initCC = 0; initCD = 0; %mol/m^3
initCon= [initCA, initCB, initCC, initCD];
v0A = 2.5*10^-5;
%v0A = 1.93*10^-5; %L/s
k2 = 0.86;
%ODE solver CSTR
[time, Ccstr] = ode45(@diffcstr11, timeSpan, initCon,[], v0A, k2);
%plot CSTR
figure (1)
plot(time, Ccstr, 'Linewidth',1.5)
grid on
xlabel('time [s]')
ylabel('concentration [mol/m^3]')
legend('c_A','c_B','c_C','c_D')
% Yield Profile
k2s = [0.01, 0.05 , 0.1, 0.86];
Y = zeros(numel(timeSpan), numel(k2s));
X = zeros(numel(timeSpan), numel(k2s));
Z = zeros(numel(timeSpan), numel(k2s));
for i = 1:numel(k2s)
[time2, C2] = ode45(@diffcstr11, timeSpan, initCon,[],v0A, k2s(i));
Z(:, i) = C2(:,4)./cao;
Y(:, i) = C2(:,3)./cao;
X(:, i) = 1 - C2(:,1)./cao;
figure (3)
plot(X(:,1), Y(:,1),'-', X(:,2), Y(:,2), '-',X(:,3), Y(:,3), '-', X(:,4), Y(:,4), '-','Linewidth',1)
grid on
hold on
ax = gca;
ax.ColorOrderIndex = 1;
plot(X(:,1), Z(:,1),'--', X(:,2), Z(:,2), '--',X(:,3), Z(:,3), '--', X(:,4), Z(:,4), '--','Linewidth',1)
hold off
xlabel('xa')
ylabel('C_C/C_A_0 C_D/C_A_0' )
legend('0.01', '0.05', '0.1', '0.86', '0.01', '0.05', '0.1', '0.86')
end
function dcdt = diffcstr11(t, C, v0A, k2)
% C(1) is cA, C(2) is cB, C(3) is cC, C(4) is cD
% parmeters
cA0 = 0.13; %mol/L
cB0 = 0.1; %mol/L
V = 30*10^-3; %L
k1 = 1; %L/(mol*s)
%k2 = 0.15; %L/(mol*s)
%v0A = 0.03; %L/s
v0B = (5*10^-5); %L/s
v0 = v0A+ v0B;%L/s
tau = V/v0;
%Rate equations
rA = -(k1*C(1)*C(2))-(k2*C(1)*C(2));
rB = rA;
rC = (k1*C(1)*C(2));
rD = (k2*C(1)*C(2));
% differential equations
dcdt = [ v0A*(cA0/V) - (C(1)/tau) - (-rA);
v0B*(cB0/V) - (C(2)/tau) - (-rB);
-(C(3)/tau) + rC;
-(C(4)/tau) + rD];
end
0 Comments
Answers (0)
See Also
Categories
Find more on Graph and Network Algorithms in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!