Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN. with Ode15s.
7 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
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 Commenti
Torsten
il 8 Mar 2022
The calling part of your program (the part above the function part) is still incomplete.
Risposta accettata
Torsten
il 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 Commenti
Torsten
il 9 Mar 2022
I live in Germany. I'm happy if I could help - no need to send a gift.
Best
Torsten.
Più risposte (0)
Vedere anche
Categorie
Scopri di più su Linear Algebra in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!