# question ..question about engine matlab

3 views (last 30 days)
Chloe Chloe on 15 Jun 2020
Commented: Rik on 16 Jun 2020
This question was flagged by Rik

Show 1 older comment
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?
Stephan on 16 Jun 2020
Nice work 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.

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.