Why do i get a empty matrix as a result?

조회 수: 3 (최근 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

카테고리

Help CenterFile Exchange에서 Programming에 대해 자세히 알아보기

제품


릴리스

R2019a

Community Treasure Hunt

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

Start Hunting!

Translated by