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!
