Way to visualise system of ODEs

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');

Answer 1

[Star Strider]

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.

[Luke Benham]

For whatever reason this only returns ys as a NaN vector?

[Star Strider]

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.