Why do i get a empty matrix as a result?

조회 수: 2(최근 30일)
Mathias Becerra Rissi
Mathias Becerra Rissi 2021년 11월 10일
답변: Sulaymon Eshkabilov 2021년 11월 10일
First of all, sorry for the typos/vocabulary, english is not my first.
The thing is, i'm modeling a reactor, with differential equations, that i defined as v(1), v(2) and v(3) for conversion, pressure and temperature. The code run smoothly in general but i need to graph another variable, thats depends on two of these, v(1) and v(2), but the result is a empty matrix.
This is the code
function dv=trioxidoo(t,v)
k1=exp(12.16-5473/v(3));
k2=exp(-9.953+8619/v(3));
k3=exp(-71.745+52596/v(3));
kp=10^(5022/v(3)-4.765);
pso2=cso20*((thetaso2-v(1))/(1+epsilon*v(1)))*(v(2)/P0)*R*T0;
po2=cso20*((thetao2-0.5*v(1))/(1+epsilon*v(1)))*(v(2)/P0)*R*T0;
pso3=cso20*((thetaso3+v(1))/(1+epsilon*v(1)))*(v(2)/P0)*R*T0;
rso2=-k1*po2*pso2*(1-(pso3/(pso2*(po2^0.5)*kp)))/(22.414*(1+k2*pso2+k3*pso3));
dv1=-rso2/Fso20;
dv2=(-B0/((1-FI)*Ac*Rhoc))*(v(2)/P0)*(v(3)/T0)*(1/(1-epsilon*v(1)));
dv3=((-rso2)*(-Dhrx))/((Fso20*thetaporcp)+(Fso20*v(1)*Dcp));
Ke1=(v(1)/(1-v(1)));
Ke2=(1-(0.5*yso20*v(1)));
Ke3=(yo20-(0.5*yso20*v(1)));
ln=log(Ke1*((Ke2/Ke3)^0.5)*(v(2))^0.5);
Te=(-B/(A+Re*ln));
This are the variables
h=25; %m
D=12; %m
T0=673;% [K]
P0=101.3; % 303.975[kpa]=3atm 1013.25kpa=10atm 1bar=101.3Kpa
R=8.314; %[kpa*m3/kmol*K]
yso20=0.11;
yo20=0.14;
yn20=0.75;
yso30=0;
Mso2=64; %kg/kmol
Mo2=32; %kg/kmol
Mn2=28; %kg/kmol
Mso3=80; %kg/kmol
PMm=Mso2*yso20+Mo2*yo20+Mn2*yn20+Mso3*yso30; %kg/kmol
rho0=P0*PMm/(R*T0); %kg/m3
v0=350000; %m3/h
m0=v0*rho0; %kg/s
F0=m0/PMm; %kmol/s
Fso20=F0*yso20;
cso20=Fso20/v0; %[kmol/m3]
thetaso2=yso20/yso20;
thetao2=yo20/yso20;
thetan2=yn20/yso20;
thetaso3=yso30/yso20;
delta=(1-0.5-1);
epsilon=delta*yso20;
B0=0.5; %[kpa/m]
FI=0.5;
Ac=(pi/4)*D^2;
Rhoc=3.357;
Dhrx=-99000; %[kJ/kmol]
Cpso2=51.62; %[kJ/kmol*K]
Cpso3=76.87; %[kJ/kmol*K]
Cpo2=33.49; %[kJ/kmol*K]
Cpn2=31.15; %[kJ/kmol*K]
Dcp=Cpso3-0.5*Cpo2-Cpo2; %[kj/kmol*K]
thetaporcp=thetaso2*Cpso2+thetaso3*Cpso3+thetao2*Cpo2+thetan2*Cpn2;%[kj/kmol*K]
wspan=[0:1:50];
[t,v]=ode45('trioxidoo',wspan,[0;P0;T0])
A=0.0938;
B=-98.591;
Re=R;
The variable or result that i need to graph is Te, but as i said earlier, it gives me a empty matrix and i dont know why.

답변(1개)

Sulaymon Eshkabilov
Sulaymon Eshkabilov 2021년 11월 10일
There are a few errs (where the parameters are defined and how the output variable dv() is defined) in defining the function file trioxidoo.m. Here is the corrected two codes:
% Code to recall and execute trioxidoo.m and plot the computed results:
clearvars; close all; clc
T0=673;% [K]
P0=101.3; % 303.975[kpa]=3atm 1013.25kpa=10atm 1bar=101.3Kpa
wspan=0:1:50;
[t,v]=ode45('trioxidoo',wspan,[0;P0;T0]);
plot(t,v(:,1), 'r-', t,v(:,2), 'g-', t,v(:,3), 'b-','linewidth', 2)
grid on; xlabel('t'); ylabel('v_1, v_2, v_3')
legend('v_1', 'v_2', 'v_3', 'location', 'best')
%% You can also consider to plot the solutions in this way.
% Since the magnitude of v_1 w.r.t v_2 and v_3 is so small.
figure
yyaxis left
plot(t,v(:,1), 'b-', 'linewidth', 2)
ylabel('v_1')
grid on
yyaxis right
plot(t,v(:,2), 'r--', t,v(:,3), 'g-.','linewidth', 2)
ylabel('v_2, v_3'), ylim([50, 1000])
xlabel('t')
legend('v_1', 'v_2', 'v_3', 'location', 'best')
Here is the corrected trioxidoo.m function file:
function dv=trioxidoo(t,v)
h=25; %m
D=12; %m
T0=673;% [K]
P0=101.3; % 303.975[kpa]=3atm 1013.25kpa=10atm 1bar=101.3Kpa
R=8.314; %[kpa*m3/kmol*K]
yso20=0.11;
yo20=0.14;
yn20=0.75;
yso30=0;
Mso2=64; %kg/kmol
Mo2=32; %kg/kmol
Mn2=28; %kg/kmol
Mso3=80; %kg/kmol
PMm=Mso2*yso20+Mo2*yo20+Mn2*yn20+Mso3*yso30; %kg/kmol
rho0=P0*PMm/(R*T0); %kg/m3
v0=350000; %m3/h
m0=v0*rho0; %kg/s
F0=m0/PMm; %kmol/s
Fso20=F0*yso20;
cso20=Fso20/v0; %[kmol/m3]
thetaso2=yso20/yso20;
thetao2=yo20/yso20;
thetan2=yn20/yso20;
thetaso3=yso30/yso20;
delta=(1-0.5-1);
epsilon=delta*yso20;
B0=0.5; %[kpa/m]
FI=0.5;
Ac=(pi/4)*D^2;
Rhoc=3.357;
Dhrx=-99000; %[kJ/kmol]
Cpso2=51.62; %[kJ/kmol*K]
Cpso3=76.87; %[kJ/kmol*K]
Cpo2=33.49; %[kJ/kmol*K]
Cpn2=31.15; %[kJ/kmol*K]
Dcp=Cpso3-0.5*Cpo2-Cpo2; %[kj/kmol*K]
thetaporcp=thetaso2*Cpso2+thetaso3*Cpso3+thetao2*Cpo2+thetan2*Cpn2;%[kj/kmol*K]
A=0.0938;
B=-98.591;
Re=R;
k1=exp(12.16-5473/v(3));
k2=exp(-9.953+8619/v(3));
k3=exp(-71.745+52596/v(3));
kp=10^(5022/v(3)-4.765);
pso2=cso20*((thetaso2-v(1))/(1+epsilon*v(1)))*(v(2)/P0)*R*T0;
po2=cso20*((thetao2-0.5*v(1))/(1+epsilon*v(1)))*(v(2)/P0)*R*T0;
pso3=cso20*((thetaso3+v(1))/(1+epsilon*v(1)))*(v(2)/P0)*R*T0;
rso2=-k1*po2*pso2*(1-(pso3/(pso2*(po2^0.5)*kp)))/(22.414*(1+k2*pso2+k3*pso3));
dv(1,:)=-rso2/Fso20;
dv(2,:)=(-B0/((1-FI)*Ac*Rhoc))*(v(2)/P0)*(v(3)/T0)*(1/(1-epsilon*v(1)));
dv(3,:)=((-rso2)*(-Dhrx))/((Fso20*thetaporcp)+(Fso20*v(1)*Dcp));
Ke1=(v(1)/(1-v(1)));
Ke2=(1-(0.5*yso20*v(1)));
Ke3=(yo20-(0.5*yso20*v(1)));
ln=log(Ke1*((Ke2/Ke3)^0.5)*(v(2))^0.5);
Te=(-B/(A+Re*ln));
end

제품


릴리스

R2019a

Community Treasure Hunt

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

Start Hunting!

Translated by