ODE Chemical Reaction Engineering with MATLAB

15 Jun 2020
1 Answer
Updated 12 Oct 2020
22 Views (30 days)

1 vote

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
Show 3 older comments Hide 3 older comments

Rik
Rik on 15 Jun 2020
The funny thing is: if you change the question title there is actually a better chance that Google will keep a cache of your original question. So below you will find what I could recover from there.
[content posted by Chloe Chloe as cached on 2020/06/15 09:32 GMT]
ODE Chemical Reaction Engineering with MATLAB
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.)
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.
Rena Berman
Rena Berman on 12 Oct 2020
(Answers Dev) Restored edit

Sign in to comment.

Sign in to answer this question.

Answers (1)

Stephan
Stephan on 15 Jun 2020
Edited: Stephan on 15 Jun 2020

0 votes

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

Tags

Asked:

Chloe Chloe
Chloe Chloe
on 15 Jun 2020

Commented:

Rena Berman
Rena Berman
on 12 Oct 2020

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!