ODE Chemical Reaction Engineering with MATLAB

36 visualizzazioni (ultimi 30 giorni)
Chloe Chloe
Chloe Chloe il 15 Giu 2020
Commentato: Rena Berman il 12 Ott 2020
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 Commenti
Rik
Rik il 16 Giu 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 il 12 Ott 2020
(Answers Dev) Restored edit

Accedi per commentare.

Risposte (1)

Stephan
Stephan il 15 Giu 2020
Modificato: Stephan il 15 Giu 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.

Categorie

Scopri di più su Chemistry in Help Center e File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by