Way to visualise system of ODEs

4 次查看（过去 30 天）
Luke Benham
Luke Benham 2021-2-18
评论： Star Strider 2021-2-18
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');

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

采纳的回答

Star Strider
Star Strider 2021-2-18
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 个评论
Luke Benham
Luke Benham 2021-2-18
For whatever reason this only returns ys as a NaN vector?
Star Strider
Star Strider 2021-2-18
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.
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.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by