Way to visualise system of ODEs
3 views (last 30 days)
Show older comments
Luke Benham
on 18 Feb 2021
Commented: Star Strider
on 18 Feb 2021
Hi,
I have a problem involving a Plug-Flow reactor where I have a system of coupled differential equations both with respect to the length travelled through the reactor.
A link to the theory; Here
I essetially have two ODEs that need to be solved;
dX/dL = n*R*A/(2*Fn)
dT/dL = n*H*A*R/(Ft*Cp)
where; H=f(T), R = f(T), Cp = f(T), n = f(T,X) etc. My code is included below.
I believe I cannot solve this analytically (dsolve) so was thinking down the lines of ode45 etc, however I'm not sure why my results so far haven't been successful.
Any help with where to begin would be amazing!
clear; close all; format shortg; opengl software; clc;
%% Input specifications
R = 8.314; %Gas Const [J/molK]
alpha = 0.5; %Constant for use in Rate Eqn
D = 1; %Diameter of Reactor
A = pi*D^2/4; %Area of Reactor
% Temp & Pressure
T_inC = 306; %[DegC]
T_inK = T_inC+273.15; %[K]
P_in = 208; %[bar] ----- subject to change
P_inA = P_in*1.01325; %[atm] -------||--------------
% Initial Mole Fractions
F_H2_in = 73.95; F_N2_in = 20.86; F_NH3_in = 1.24;
F_T_in = F_N2_in+F_H2_in+F_NH3_in;
Y_N2 = F_N2_in/F_T_in; Y_H2 = F_H2_in/F_T_in;
Y_NH3 = F_NH3_in/F_T_in;
sumY = Y_N2+Y_H2+Y_NH3; %Make sure sum(Yi)=1
%% Equilibrium
disp('------------------------------------------')
disp('Starting Simulation');
syms T(L) X(L)
P = P_inA;
log10Ka = -2.691122.*log(T)-5.519265e-5.*T...
+1.848863e-6.*T.^2+2001.6./T+2.689;
Ka = 10.^(log10Ka);
% Arrhenius
A0 = 8.849e14; %Arrhenius Constant
Ea = 170560; %Activation Energy [J/mol]
k = A0.*exp(Ea./(R.*T)); %Rate Constant
% Fugacities & Activities
phi_N2 = 0.934+0.203e-3.*T+0.296e-3*P-...
0.271e-6.*T.^2+0.4775e-6*P^2; %Fugacity Coeff - N2 (Reactant)
y_N2 = Y_N2*(1-X)./(1-2.*X.*Y_N2); %Mol of N2
a_N2 = phi_N2.*y_N2.*P; %Activity of N2
phi_H2 = exp(exp(-3.84.*T.^0.125+0.541)*P...
-exp(-0.126.*T.^0.5-15.98)*P^2+...
300*exp(-0.0119.*T-5.941)*exp(+P/300)); %Fugacity Coeff - H2 (Reactant)
y_H2 = (Y_H2-3.*X.*Y_N2)./(1-2.*X.*Y_N2); %Mol of H2
a_H2 = phi_H2.*y_H2.*P; %Activity of N2
phi_NH3 = 0.144+0.203e-2.*T-0.449e-3*P-...
0.114e-5.*T.^2+0.276e-6*P^2; %Fugactiy Coeff - NH3 (Product)
y_NH3 = (Y_NH3+2.*X.*Y_N2)./(1-2.*X.*Y_N2); %Mol of NH3
a_NH3 = phi_NH3.*y_NH3.*P; %Activity of N2
% Other Vars
C_pT = 35.31+0.02*T+0.00000694*T^2-0.0056*P+0.000014*P*T; %Specific Heat Cap. of Mixture
F_N2 = F_N2_in*(1-X);
F_H2 = F_H2_in-3*X*F_N2_in;
F_NH3 = F_NH3_in+2*X*F_N2_in;
F_T = F_N2+F_H2+F_NH3;
% Final Vars
dH_R = 4.184*(-(0.54526+840.609./T+459.734e7./(T.^3))*P-...
5.34685.*T-0.2525e-3.*T.^2+1.69167e-6.*T.^3-9157.09); %Heat of Reaction
R_NH3 = 2.*k.*(Ka.^2.*a_N2.*(a_H2.^3./a_NH3.^2).^alpha-...
(a_NH3.^2./a_H2.^3).^(1-alpha)); %Rate of Reaction
eps = -17.539096+0.07697849.*T+6.900548.*T.^2-1.082e-4.*T.^3-...
26.4247.*T.^4+4.9276e-8.*T.^5+38.937.*X.^3; %Catalyst Effectiveness
disp('Variables Simulated');
%% Solve Differential Equations
disp('Starting ODE solving');
ODE1 = diff(X) == eps.*R_NH3.*A./(2.*F_N2);
ODE2 = diff(T) == eps.*(-dH_R).*A.*R_NH3./(F_T.*C_pT);
ODEs = [ODE1; ODE2];
V1 = odeToVectorField(ODE1,ODE2);
F1 = matlabFunction(V1,'vars',{'L','Y','X','T'});
cond1 = T(0) == T_inK;
cond2 = X(0) == 0;
conds = [cond1; cond2];
Lspan = [0, 30];
y0 = [T_inK, 0];
s = ode45(F1)
disp('ODEs Solved');
0 Comments
Accepted Answer
Star Strider
on 18 Feb 2021
Try your code with these changes:
[V1,Sbs] = odeToVectorField(ODE1,ODE2);
F1 = matlabFunction(V1,'vars',{'t','Y'});
cond1 = T(0) == T_inK;
cond2 = X(0) == 0;
conds = [cond1; cond2];
Lspan = [0, 30];
y0 = [T_inK, 0];
[t,ys] = ode45(F1, Lspan, y0);
disp('ODEs Solved');
figure
yyaxis left
plot(t,ys(:,1))
yyaxis right
plot(t, ys(:,2))
legend(string(Sbs), 'Location','E')
grid
The code prior to that is unchanged. The first argument to ‘F1’ must be the independent variable vector, even if you don’t use it. Since the only variable in ‘F1’ is ‘Y’ otherwise (as the default dependent variable), that’s all matlabFunction needs to know.
2 Comments
Star Strider
on 18 Feb 2021
Not when I ran it.
The beginning of the integrated result is a (Nx2) complex matrix, so it plots up to about
seconds (or whatever the independent variable is). After that everything is NaN. The NaN values are the result of 0/0, Inf/Inf or some such. I have no idea what you’re doing (and I have no significant experience with Plug-Flow reactors), so I can’t help you with that, and defer to you to troubleshoot it.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/524117/image.png)
I got your code to work, and integrate your ODEs, which is what you requested help with. That’s the best I can do.
Also, I noticed that you’re using:
opengl software
If you have an AMD graphics card and are having problems with it, update the driver. There were known issues with AMD drivers several months ago that have been corrected with the latest driver updates.
More Answers (0)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!