Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN. with Ode15s.
6 views (last 30 days)
Show older comments
I'm writing a code about kinetic reaction in reactor solve 9 differiential equation, but now i'm stuck in Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN.
Here this code that I write
P0 = [0.5] ;%bar (Intial pressure)
Ta0 = [823] ;%K (Intial air heat temperature)
repeat = [1] ;%Repeat condition Ta constant and Ta non constant
sheet = 0 ;%Sheet in excel
condition = 0 ;%Number of condition
for a = 1:length(repeat)
for i = length(P0)
for j = length(Ta0)
T0 = Ta0(j) ;%K (Intial temperature before entering bed)
yj0 = [0.1 0.1 0 0 0 0 0.8]'
pj0 = yj0.*P0(i) ;%bar (Partial pressure)
Fj0 = (pj0.*v0_feed)./(R_feed.*T0) ;%mol/s (Molar flowrate)
Lspan = [0 0.03] ;%m (Bed length range)
dep_var0 = [Fj0' T0 Ta0(j)] ;%[mol/s,K,K] (Intial dependent variable)
condition = condition + 1 ;
var = [P0(i) Ta0(j) condition T0] ;%[bar,K,-,K] (Another variable that nessery)
[L ,dep_var] = ode23s(@func_Ni,Lspan,dep_var0,[],var) ;
FCH4 = dep_var(:,1) ;%mol/s
FCO2 = dep_var(:,2) ;%mol/s
FCO = dep_var(:,3) ;%mol/s
FH2 = dep_var(:,4) ;%mol/s
FH2O = dep_var(:,5) ;%mol/s
FC = dep_var(:,6) ;%mol/s
FN2 = dep_var(:,7) ;%mol/s
T = dep_var(:,8) ;%K
Ta = dep_var(:,9) ;%K
For the functions
function dydx = func_Ni(y,x,op)
load Data_import_Ni
FCH4 = x(1,1) ;
FCO2 = x(2,1) ;
FCO = x(3,1) ;
FH2 = x(4,1) ;
FH2O = x(5,1) ;
FC = x(6,1) ;
FN2 = x(7,1) ;
T = x(8,1) ;
Ta = x(9,1) ;
P0 = op(1) ;%bar (Intial pressure)
Tac = op(2) ;%K (Intial air heat temperature constant)
condition = op(3) ;
T0 = op(4) ;%K (Intial temperature)
FT = FCH4+FCO2+FCO+FH2+FH2O+FC+FN2 ;%mol/s (Total molar flowrate)
%% Partial pressure (bar)
pCH4 = P0*(FCH4/FT)*(T0/T) ;
pCO2 = P0*(FCO2/FT)*(T0/T) ;
pCO = P0*(FCO/FT)*(T0/T) ;
pH2 = P0*(FH2/FT)*(T0/T) ;
pH2O = P0*(FH2O/FT)*(T0/T) ;
pC = P0*(FC/FT)*(T0/T) ;
pN2 = P0*(FN2/FT)*(T0/T) ;
%% Heat capacity (J/(mol*K))
cpj = heatcp(:,1) + heatcp(:,2)*T + heatcp(:,3)*(T^2) + heatcp(:,4)*(T^-2) ;
%% Del_Heat of reaction
delH1 = hrxi_Tr(1)+((2.*cpj(3)+2.*cpj(4)-cpj(1)-cpj(2)).*(T-Tr)) ;
delH2 = hrxi_Tr(2)+((cpj(3)+cpj(5)-cpj(2)-cpj(4)).*(T-Tr)) ;
delH3 = hrxi_Tr(3)+((cpj(6)+2.*cpj(4)-cpj(1)).*(T-Tr)) ;
delH4 = hrxi_Tr(4)+((cpj(3)+cpj(4)-cpj(5)-cpj(6)).*(T-Tr)) ;
delH5 = hrxi_Tr(5)+((2*cpj(3)-cpj(6)-cpj(2)).*(T-Tr)) ;
%% Reaction rate constant (mol/(kg*s))
ki = k0i.*exp(-Eai./T) ;
%% Equilibrium constant
Kpi = Kp0i.*exp(-dHrxi_Tr./T) ;%[bar^2,-,bar,bar,bar]
%% Adsorption equilibrium constant
KCH4_1 = K0j1(1)*exp(hdx1(1)/T) ;%mol/(kg*s*bar)
KCO2_1 = K0j1(2)*exp(hdx1(2)/T) ;%mol/(kg*s*bar)
KCO2_2 = K0j2(2)*exp(hdx2(2)/T) ;%bar^-1
KH2_2 = K0j2(4)*exp(hdx2(4)/T) ;%bar^-1
KCH4_3 = K0j3(1)*exp(hdx3(1)/T) ;%bar^-1
KH2_3 = K0j3(4)*exp(hdx3(4)/T) ;%bar^1.5
KCH4_4 = K0j4(1)*exp(hdx4(1)/T) ;%bar^-1
KH2_4 = K0j4(4)*exp(hdx4(4)/T) ;%bar^1.5
KH2O_4 = K0j4(5)*exp(hdx4(5)/T) ;%-
KCO2_5 = K0j5(2)*exp(hdx5(2)/T) ;%bar
KCO_5 = K0j5(3)*exp(hdx5(3)/T) ;%bar^-1
%% Rate reaction (mol/(kg*s)
r1 = ((ki(1)*KCH4_1*KCO2_1*pCH4*pCO2)/((1+KCO2_1*pCO2+KCH4_1*pCH4)^2))*(1-(((pCO*pH2)^2)/(Kpi(1)*pCH4*pCO2))) ;
r2 = ((ki(2)*KCO2_2*KH2_2*pCO2*pH2)/((1+KCO2_2*pCO2+KH2_2*pH2)^2))*(1-((pCO*pH2O)/(Kpi(2)*pCO2*pH2))) ;
r3 = ((ki(3)*KCH4_3)*(pCH4-((pH2^2)/Kpi(3))))/((1+KCH4_3*pCH4+(pH2^1.5/KH2_3))^2) ;
r4 = ((ki(4)/KH2O_4)*((pH2O/pH2)-(pCO/Kpi(4))))/((1+KCH4_4*pCH4+(pH2O/(KH2O_4*pH2))+(pH2^1.5/KH2_4))^2) ;
r5 = ((ki(5)/(KCO_5*KCO2_5))*((pCO2/pCO)-(pCO/Kpi(5))))/((1+KCO_5*pCO+(pCO2/(KCO_5*KCO2_5*pCO)))^2) ;
%% rnet of species (mol/(kg*s))
rn_CH4 = -r1-r3 ;
rn_CO2 = -r1-r2-r5 ;
rn_CO = 2*r1+r2+r4+2*r5 ;
rn_H2 = 2*r1-r2+2*r3+r4 ;
rn_H2O = r2-r4 ;
rn_C = r3-r4-r5 ;
rn_N2 = 0 ;
%%
sumFjcpj = FCH4*cpj(1)+ FCO2*cpj(2)+ FCO*cpj(3)+ FH2*cpj(4)+ FH2O*cpj(5)+ FC*cpj(6)+ FN2*cpj(7) ;%Watt/K
sumrihi = r1*delH1+ r2*delH2+ r3*delH3+ r4*delH4+ r5*delH5 ;%Watt/kg
%% ODE equations
dFCH4dL = Rho_bed*Ac*rn_CH4 ;%mol/(m*s)
dFCO2dL = Rho_bed*Ac*rn_CO2 ;%mol/(m*s)
dFCOdL = Rho_bed*Ac*rn_CO ;%mol/(m*s)
dFH2dL = Rho_bed*Ac*rn_H2 ;%mol/(m*s)
dFH2OdL = Rho_bed*Ac*rn_H2O ;%mol/(m*s)
dFCdL = Rho_bed*Ac*rn_C ;%mol/(m*s)
dFN2dL = Rho_bed*Ac*rn_N2 ;%mol/(m*s)
if condition <= length(T0)*length(P0)
dTdL = ((Rho_bed*Ac*sumrihi)-(U*ain*Ac*(T-Tac)))/sumFjcpj ;%K/m
dTadL = 0 ;%K/m
else
dTdL = ((Rho_bed*Ac*sumrihi)-(U*ain*Ac*(T-Ta)))/sumFjcpj ;%K/m
dTadL = (U*aout*Ac_air*(T-Ta))/(mc_air*cpj(8)) ;%K/m
end
%% Calculations
dydx(1,1) = dFCH4dL ;
dydx(2,1) = dFCO2dL ;
dydx(3,1) = dFCOdL ;
dydx(4,1) = dFH2dL ;
dydx(5,1) = dFH2OdL ;
dydx(6,1) = dFCdL ;
dydx(7,1) = dFN2dL ;
dydx(8,1) = dTdL ;
dydx(9,1) = dTadL ;
end
7 Comments
Torsten
on 8 Mar 2022
The calling part of your program (the part above the function part) is still incomplete.
Accepted Answer
Torsten
on 9 Mar 2022
P0 = [0.5] %bar (Intial pressure)
Ta0 = [823] ;%K (Intial air heat temperature)
repeat = [1] ;%Repeat condition Ta constant and Ta non constant
sheet = 0 ;%Sheet in excel
condition = 0 ;%Number of condition
for a = 1:length(repeat)
for i = length(P0)
for j = length(Ta0)
T0 = Ta0(j) ;%K (Intial temperature before entering bed)
yj0 = [0.1 0.1 1e-8 1e-8 1e-8 1e-8 0.8]'
pj0 = yj0.*P0(i) ;%bar (Partial pressure)
Fj0 = (pj0.*v0_feed)./(R_feed.*T0) ;%mol/s (Molar flowrate)
Lspan = [0 0.03] ;%m (Bed length range)
dep_var0 = [Fj0' T0 Ta0(j)] ;%[mol/s,K,K] (Intial dependent variable)
condition = condition + 1 ;
var = [P0(i) Ta0(j) condition T0] ;%[bar,K,-,K] (Another variable that nessery)
[L ,dep_var] = ode15s(@(y,x)func_Ni(y,x,var),Lspan,dep_var0) ;
FCH4 = dep_var(:,1) ;%mol/s
FCO2 = dep_var(:,2) ;%mol/s
FCO = dep_var(:,3) ;%mol/s
FH2 = dep_var(:,4) ;%mol/s
FH2O = dep_var(:,5) ;%mol/s
FC = dep_var(:,6) ;%mol/s
FN2 = dep_var(:,7) ;%mol/s
T = dep_var(:,8) ;%K
Ta = dep_var(:,9) ;%K
end
end
end
By changing the initial conditions for the species and by modifying the call to the ODE solver, NaN values for the derivatives are no longer computed.
But you'll have to study the output from your function more carefully to make the solver really advance in time.
2 Comments
Torsten
on 9 Mar 2022
I live in Germany. I'm happy if I could help - no need to send a gift.
Best
Torsten.
More Answers (0)
See Also
Categories
Find more on Linear Algebra in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!