ODE Chemical Reaction Engineering with MATLAB

Question

Chloe Chloe 2020-6-15

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/548430-ode-chemical-reaction-engineering-with-matlab

评论： Rena Berman 2020-10-12

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 个评论
显示 3更早的评论隐藏 3更早的评论

Rik 2020-6-15

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 2020-6-16

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 2020-10-12

(Answers Dev) Restored edit

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Stephan 2020-6-15

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/548430-ode-chemical-reaction-engineering-with-matlab#answer_451401

编辑：Stephan 2020-6-15

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