Why do I receive incorrect result from ode45?
3 views (last 30 days)
Show older comments
Hello, I'm trying to solve this system of differential equations using ode45 function.
The code I designed is written below.
concentration.m file
%This models concentration of CH4 over FBR length
function diffeqs=concentration(~,var)
cH2 = var(1);
cCH4 = var(2); %define variable
cC2H4 = var(3); %Imporant!!
cC2H6 = var(4);
cC6H6 = var(5);
cC7H8 = var(6);
cC10H8 = var(7);
a1 = var (8);
a2 = var (9);
%Rate constants
k1 = 0.01831;
k2 = 1231;
k3 = 11340000;
k4 = 1.697;
k5 = 706300;
k6 = 2917;
%Equilibrium constants
K1 = (3.07941*10^-5);
K2 = 5.094435141; %chane k2 and k6
K3 = 721589.6019;
K4 = 142356.6252;
K5 = 0.007763912;
K6 = 9405.352961;
%Adsorption Equilibrium constants
K1H2 = 0.07581; %Site 1
K1CH4 = 0.09664;
K1C2H4 = 0.07944;
K1C2H6 = 0.08095;
K2H2 = 0.09396; %Site 2
K2CH4 = 0.1189;
K2C2H4 = 0.2269;
K2C2H6 = 0.05123;
K2C6H6 = 0.04881;
K2C7H8 = 0.05219;
K2C10H8 = 0.1175;
%Conversion between concentration, partial pressure and molar flowrate
R = 8.314;
T = 973.15; %in K
PH2 = (cH2*R*T*10^-5); %bar pressure can be drop
PCH4 = (cCH4*R*T*10^-5);
PC2H4 = (cC2H4*R*T*10^-5);
PC2H6 = (cC2H6*R*T*10^-5);
PC6H6 = (cC6H6*R*T*10^-5);
PC7H8 = (cC7H8*R*T*10^-5);
PC10H8 = (cC10H8*R*T*10^-5);
%DEN
DEN1 = (1 + K1H2*PH2 + K1CH4*PCH4 + K1C2H4*PC2H4 + K1C2H6*PC2H6);
DEN2 = (1 + K2H2*PH2 + K2CH4*PCH4 + K2C2H4*PC2H4 + K2C2H6*PC2H6 + K2C6H6*PC6H6 + K2C7H8*PC7H8 + K2C10H8*PC10H8);
%Deactivation
diffeqs(8,1) = -1.5;
diffeqs(9,1) = -0.0091;
%Reaction rate calculation
r1 = (k1*(PCH4^2 - PC2H4*PH2^2/K1)*a1/DEN1);
r2 = (k2*(PC2H4^4 - PC6H6*PC2H6*PH2^2/K2)*a2/DEN2);
r3 = (k3*(PC2H4^3 - PC6H6*PH2^3/K3)*a2/DEN2);
r4 = (k4*(PC6H6*PCH4 - PC7H8*PH2/K4)*a2/DEN2);
r5 = (k5*(PC6H6*PC2H4^2 -PC10H8*PH2^3/K5)*a2/DEN2);
r6 = (k6*(PC2H4*PH2 - PC2H6/K6)*a1/DEN1);
%Constants
us = 0.3068; % Feed velocity meter/h
pB = 500; %kg/m3
%Differential eqs
diffeqs(1,1) = pB*(2*r1 +2*r2 +3*r3 +r4 +3*r5 -r6)/us; %H2
diffeqs(2,1) = pB*(-2*r1 -r4)/us; %CH4
diffeqs(3,1) = pB*(r1 -4*r2 -3*r3 -2*r5 -r6)/us; %C2H4
diffeqs(4,1) = pB*(r2 +r6)/us; %C2H6
diffeqs(5,1) = pB*(r2 + r3 - r4 -r5)/us; %C6H6
diffeqs(6,1) = pB*r4/us; %C7H8
diffeqs(7,1) = pB*r5/us; %C8H10
end
solving.m file
%Script to solve DEs
range =[0 0.39]; %fix to real length
ICs = [0, 0.9, 0, 0, 0, 0, 0, 1, 1]; %mol/m3
[z,conc]=ode45(@concentration,range,ICs);
figure
plot(z,conc(:,1)
%xaxis('z, m')
%ylabel('conc, mol/m3')
From the rate constant values, k3 is much greater than k6 so it expects to get cC6H6 larger than cC2H6 but I got the otherwise result.
Please help me check if the incorrect result coming from the code itself or it might be the logic of differential eqn system.
Many thanks!
0 Comments
Answers (1)
Sulaymon Eshkabilov
on 24 Sep 2021
Edited: Sulaymon Eshkabilov
on 24 Sep 2021
You can try to tighten the relative and absolute error tolerances, e.g.:
OPTS = odeset('reltol', 1e-5, 'abstol', 1e-7);
[z,conc]=ode45(@concentration, Range,ICs, OPTS);
Another point (a good pratice) is not use MATLAB's built in fcn names. In your code, you have range() that should be renamed, e.g.: Range
One typo - missing ) in plot(z,conc(:,1) that was missed due to copy and paste.
See Also
Categories
Find more on Loops and Conditional Statements 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!