ODE Chemical Reaction Engineering with MATLAB

36 views (last 30 days)
I am solving a chemical reaction engineering Problem like the following.
function [W,Fa] = PBR_Isothermal
clc; clear;
%A = Dammitol
%B = Valualdehyde
%C = Oxygen
%D = Carbon Dioxide
%E = Water
%I = Nitrogen
%Feed based on inlet as 100 kmol/h
Fa0 = 0.1*100;
Fb0 = 0;
Fc0 = 0.07*100;
Fd0 = 0;
Fe0 = 0.02*100;
Fi0 = 0.81*100;
[W, Fa] = ode45(@FunA,[0 1000],[Fa0 Fb0 Fc0 Fd0 Fe0 Fi0]);
% [Fb W] = ode45(FunA,[0 1000],Fb0);
% [Fc W] = ode45(FunA,[0 1000],Fc0);
% [Fd W] = ode45(FunA,[0 1000],Fd0);
% [Fe W] = ode45(FunA,[0 1000],Fe0);
end
function A = FunA(W,F)
%Total Pressure
Ptot = 114.3;
%Defining Constant
FT = F(1)+F(2)+F(3)+F(4)+F(5)+F(6);
RT = (1/353-1/373)/1.987;
%Partial Pressure of Each
PD = (F(1)/FT)*Ptot;
PV = (F(2)/FT)*Ptot;
PO2 = (F(3)/FT)*Ptot;
PCO = (F(4)/FT)*Ptot;
PWA = (F(5)/FT)*Ptot;
PN2 = (F(6)/FT)*Ptot;
%Constants
k1 = 1.771*10^-3;
k2 = 23295;
k3 = 0.5;
k4 = 1.0;
k5 = 0.8184;
k6 = 0.0;
k7 = 0.5;
k8 = 0.2314;
k9 = 0.0;
k10 = 1.0;
k11 = 1.25;
k12 = 0.0;
k13 = 2.795*10^-4;
k14 = 33000;
k15 = 0.5;
k16 = 2.0;
k17 = 2.0;
%Rate of Reaction
r1 = (k1*exp(-k2*RT)*PO2^k3*PD^k4)/(1+k5*exp(k6*RT)*PO2^k7+k8*exp(k9*RT)*PD^k10+k11*PV^k17*exp(k12*RT));
r2 = (k13*exp(-k14*RT)*PO2^k15*PV^k16)/(1+k5*exp(k6*RT)*PO2^k7+k8*exp(k9*RT)*PD^k10+k11*PV^k17*exp(k12*RT));
%Molar Flow Rate of Each Component
A(1) = -r1;
A(2) = r1-r2;
A(3) = -0.5*r1-2.5*r2;
A(4) = 2*r2;
A(5) = r1+2*r2;
Fi0 = 0.81*100;
A(6) = Fi0;
A=A';
end
I got till here but I'm having trouble working with the rest of the code. Can anyone suggest how I should implement the rest of the questions?
(Please just ignore the complicated coefficients. I've already made the library.)
  5 Comments
Rik
Rik on 16 Jun 2020
Regarding your flag ("i wish to remove this question due to plagirism"): how is this plagiarism? Your code doesn't claim that you are the sole creator (although it doesn't provide a reference). I could imagine those screenshots being copyright infringment, but that is a completely different thing.

Sign in to comment.

Answers (1)

Stephan
Stephan on 15 Jun 2020
Edited: Stephan on 15 Jun 2020
Your code works - you only need to call your function and plot the results:
% Call the function and save results in W and Fa
[W,Fa] = PBR_Isothermal;
% plot results
subplot(3,2,1)
plot(W,Fa(:,1))
title('Fa(1)')
subplot(3,2,2)
plot(W,Fa(:,2))
title('Fa(2)')
subplot(3,2,3)
plot(W,Fa(:,3))
title('Fa(3)')
subplot(3,2,4)
plot(W,Fa(:,4))
title('Fa(4)')
subplot(3,2,5)
plot(W,Fa(:,5))
title('Fa(5)')
subplot(3,2,6)
plot(W,Fa(:,6))
title('Fa(6)')
function [W,Fa] = PBR_Isothermal
clc; clear;
%A = Dammitol
%B = Valualdehyde
%C = Oxygen
%D = Carbon Dioxide
%E = Water
%I = Nitrogen
%Feed based on inlet as 100 kmol/h
Fa0 = 0.1*100;
Fb0 = 0;
Fc0 = 0.07*100;
Fd0 = 0;
Fe0 = 0.02*100;
Fi0 = 0.81*100;
[W, Fa] = ode45(@FunA,[0 1000],[Fa0 Fb0 Fc0 Fd0 Fe0 Fi0]);
% [Fb W] = ode45(FunA,[0 1000],Fb0);
% [Fc W] = ode45(FunA,[0 1000],Fc0);
% [Fd W] = ode45(FunA,[0 1000],Fd0);
% [Fe W] = ode45(FunA,[0 1000],Fe0);
end
function A = FunA(~,F)
%Total Pressure
Ptot = 114.3;
%Defining Constant
FT = F(1)+F(2)+F(3)+F(4)+F(5)+F(6);
RT = (1/353-1/373)/1.987;
%Partial Pressure of Each
PD = (F(1)/FT)*Ptot;
PV = (F(2)/FT)*Ptot;
PO2 = (F(3)/FT)*Ptot;
PCO = (F(4)/FT)*Ptot;
PWA = (F(5)/FT)*Ptot;
PN2 = (F(6)/FT)*Ptot;
%Constants
k1 = 1.771*10^-3;
k2 = 23295;
k3 = 0.5;
k4 = 1.0;
k5 = 0.8184;
k6 = 0.0;
k7 = 0.5;
k8 = 0.2314;
k9 = 0.0;
k10 = 1.0;
k11 = 1.25;
k12 = 0.0;
k13 = 2.795*10^-4;
k14 = 33000;
k15 = 0.5;
k16 = 2.0;
k17 = 2.0;
%Rate of Reaction
r1 = (k1*exp(-k2*RT)*PO2^k3*PD^k4)/(1+k5*exp(k6*RT)*PO2^k7+k8*exp(k9*RT)*PD^k10+k11*PV^k17*exp(k12*RT));
r2 = (k13*exp(-k14*RT)*PO2^k15*PV^k16)/(1+k5*exp(k6*RT)*PO2^k7+k8*exp(k9*RT)*PD^k10+k11*PV^k17*exp(k12*RT));
%Molar Flow Rate of Each Component
A(1) = -r1;
A(2) = r1-r2;
A(3) = -0.5*r1-2.5*r2;
A(4) = 2*r2;
A(5) = r1+2*r2;
Fi0 = 0.81*100;
A(6) = Fi0;
A=A';
end
Note that PCO, PWA and PN2 are not used inside your function.

Categories

Find more on Chemistry 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!