Differential Equation Problem using ODE45

In below script, I was using the equation of dF/dz [F=CO2, CO, H2O, H2, DME, CH3OH] which written as A(1) until A(6) to solve the molar flowrate of each component. It was assumed that the feed only contained CO2 and H2 in the inlet flowrate. The mole fraction of CO2 and H2 was assumed to 0.25 and 0.75 with a total flowrate of 0.2667 kmol/s [960 kmol/hr].

In this case, the script was successfully run. However, the mole flowrate obtained for each component was NaN that results to no reaction occured/error. Appreciate if you could help me and give suggestion to improve the script. Thank you in advance.

The script:

function A = test_A(Z,F)

%Fixed parameters

Ptot=50;

T=523;

pc=1900;

%Defining constant

FT=F(1)+F(2)+F(3)+F(4)+F(5)+F(6)+F(7);

RT=523*8.314;

%Feed based on inlet as 0.2667 kmol/s

F(1)=0.25*0.2667;

F(2)=0.75*0.2667;

F(3)=0*0.2667;

F(4)=0*0.2667;

F(5)=0*0.2667;

F(6)=0*0.2667;

F(7)=0*0.2667;

%Partial pressure for each component

pCO2=(F(1)/FT)*Ptot;

pH2=(F(2)/FT)*Ptot;

pCO=(F(3)/FT)*Ptot;

pH2O=(F(4)/FT)*Ptot;

pDME=(F(5)/FT)*Ptot;

pCH3OH=(F(6)/FT)*Ptot

pN2=(F(7)/FT)*Ptot;

%Kinetic parameters

k1=35.45*exp(-1.7069e4/(RT));

k2=7.3976*exp(-2.0436e4/(RT));

k3=8.2894e4*exp(-5.2940e4/(RT));

%Adsorption constant

kH2=0.249*exp(3.4394e4/(RT));

kCO2=1.02*10^-7*exp(6.74*10^4/(RT));

kCO=7.99*10^-7*exp(5.81*10^4/(RT));

%Equilibrium constant

kp1=exp((4213/T)-(5.752*log(T))-(1.707e-3*T)+(2.682e-6*(T^2))-(7.232e-10*(T^3))+17.6);

kp2=exp((2167/T)-(0.5194*log(T))+(1.037e-3*T)-(2.331e-7*(T^2))-1.2777);

kp3=exp((4019/T)+(3.707*log(T))-(2.783e-3*T)+(3.81e-7*(T^2))-(6.561e4/(T^3))-26.64);

%Rate of Reaction

r1=k1*(pCO2*pH2*(1-((pCH3OH*pH2O)/(kp1*pCO2*pH2^3))))/(1+kCO2*pCO2+kCO*pCO+sqrt(kH2*pH2))^3;

r2=k2*((((pCO2*pH2)/(kp2*pCO))-pH2O)/((1+(kCO2*pCO2)+(kCO*pCO)+(sqrt(kH2*pH2)))^3));

r3=k3*((pCH3OH^2/pH2O)-(pDME/kp3)); %MeOH dehydration

%Trans-membrane molar flux

permN2=1.3333e-8;

permH2O=1e-6;

prN2=0*Ptot;

prH2O=0*Ptot;

ppN2=0.79320;

ppH2O=0.20680;

JN2=permN2*(prN2-pN2);

JH2O=permH2O*(prH2O-pH2O);

%Mass balance for each component

vf=0.33;

Dsi=0.0198;

Dmo=14e-3;

%Mass balance for each component

A(1)=(pc*(1-vf)*(r1-r2+r3)*(pi*Dmo^2)-(JH2O*pi*Dmo));

A(2)=(-pc*(1-vf)*(3*r1-r2)*((pi/4)*(Dsi^2-Dmo^2)));

A(3)=(pc*(1-vf)*(r2)*((pi/4)*(Dsi^2-Dmo^2)));

A(4)=(-pc*(1-vf)*(r1-r2)*((pi/4)*(Dsi^2-Dmo^2)));

A(5)=(pc*(1-vf)*(r1-2*r3)*((pi/4)*(Dsi^2-Dmo^2)));

A(6)=(pc*(1-vf)*(r3)*((pi/4)*(Dsi^2-Dmo^2)));

A=A';

[Torsten]

You didn't include the script, but only the function to be used by ode45.

Answer 1

[Sam Chak]

Hi @NURUL AINA SYAHIRAH

First things first, check the values computed by the differential equations. All differential states return either NaN or Inf because the equations contain r2, or r3, or both. In case you don't know, r2 returns Inf due to a division by zero problem (pCO = 0), and r3 returns NaN due to the indeterminate term in pCH3OH^2/pH2O. Fix these two issues and the code should run.

Zspan   = [0 1];
F0      = [1 1 1 1 1 1];
[Z, F]  = ode45(@test_A, Zspan, F0);
plot(Z, F)

function A = test_A(Z,F)

    %Fixed parameters
    Ptot=50;
    T=523;
    pc=1900;
    
    %Feed based on inlet as 0.2667 kmol/s
    F(1)=0.25*0.2667;
    F(2)=0.75*0.2667;
    F(3)=0*0.2667;
    F(4)=0*0.2667;
    F(5)=0*0.2667;
    F(6)=0*0.2667;
    F(7)=0*0.2667;
    
    %Defining constant
    FT=F(1)+F(2)+F(3)+F(4)+F(5)+F(6)+F(7);
    RT=523*8.314;
      
    %Partial pressure for each component
    pCO2=(F(1)/FT)*Ptot;    % 12.5
    pH2=(F(2)/FT)*Ptot;     % 37.5
    pCO=(F(3)/FT)*Ptot  + 0.1;     % 0 plus a small non-zero value
    pH2O=(F(4)/FT)*Ptot + 0.1;     % 0 plus a small non-zero value
    pDME=(F(5)/FT)*Ptot;    % 0
    pCH3OH=(F(6)/FT)*Ptot;  % 0
    pN2=(F(7)/FT)*Ptot;     % 0
    
    %Kinetic parameters
    k1=35.45*exp(-1.7069e4/(RT)); 
    k2=7.3976*exp(-2.0436e4/(RT)); 
    k3=8.2894e4*exp(-5.2940e4/(RT));
    
    %Adsorption constant
    kH2=0.249*exp(3.4394e4/(RT));
    kCO2=1.02*10^-7*exp(6.74*10^4/(RT));
    kCO=7.99*10^-7*exp(5.81*10^4/(RT));
    
    %Equilibrium constant
    kp1=exp((4213/T)-(5.752*log(T))-(1.707e-3*T)+(2.682e-6*(T^2))-(7.232e-10*(T^3))+17.6);
    kp2=exp((2167/T)-(0.5194*log(T))+(1.037e-3*T)-(2.331e-7*(T^2))-1.2777);
    kp3=exp((4019/T)+(3.707*log(T))-(2.783e-3*T)+(3.81e-7*(T^2))-(6.561e4/(T^3))-26.64);
    
    %Rate of Reaction
    r1=k1*(pCO2*pH2*(1-((pCH3OH*pH2O)/(kp1*pCO2*pH2^3))))/(1+kCO2*pCO2+kCO*pCO+sqrt(kH2*pH2))^3;
    r2=k2*( (((pCO2*pH2)/(kp2*pCO)) - pH2O) / ((1 + kCO2*pCO2 + kCO*pCO + sqrt(kH2*pH2))^3) );
    r3=k3*((pCH3OH^2/pH2O)-(pDME/kp3)); %MeOH dehydration
    
    %Trans-membrane molar flux
    permN2=1.3333e-8;
    permH2O=1e-6;
    prN2=0*Ptot;
    prH2O=0*Ptot;
    ppN2=0.79320;
    ppH2O=0.20680;
    JN2=permN2*(prN2-pN2);
    JH2O=permH2O*(prH2O-pH2O);
    
    %Mass balance for each component
    vf=0.33; 
    Dsi=0.0198;
    Dmo=14e-3;
    
    %Mass balance for each component
    A(1)=( pc*(1-vf)*(r1 - r2 + r3)*(pi*Dmo^2)-(JH2O*pi*Dmo));   % NaN
    A(2)=(-pc*(1-vf)*(3*r1 - r2)   *((pi/4)*(Dsi^2-Dmo^2)));     % Inf
    A(3)=( pc*(1-vf)*(r2)*((pi/4)  *(Dsi^2-Dmo^2)));             % Inf
    A(4)=(-pc*(1-vf)*(r1 - r2)     *((pi/4)*(Dsi^2-Dmo^2)));     % Inf
    A(5)=( pc*(1-vf)*(r1 - 2*r3)   *((pi/4)*(Dsi^2-Dmo^2)));     % NaN
    A(6)=( pc*(1-vf)*(r3)*((pi/4)  *(Dsi^2-Dmo^2)));             % NaN
    A=A';

end

[NURUL AINA SYAHIRAH]

Thank you!

[Sam Chak]

You are welcome, @NURUL AINA SYAHIRAH.