plotting profiles from differential equations
4 ビュー (過去 30 日間)
古いコメントを表示
Can anyone help me correct my code to produce a graphical solution to my differrential equations? Below is my code which is saved as Reactor3.m. I have attached my differretial equations in a pdf file along with the intial conditions.
function [dFdV,dTdV] = Reactor3(V,F,T)
%Reactor3 function is used to obtain mass and energy balances across the
%reactor
%Constants
P=2200; %inlet Pressure (kPa)
dcat=850000; %catalyst density g/m3
U=0.17603*3600; %overall heat transfer co (kj/m2 h K)
a=80; %(m^-1)
Tc=433.15; %coolant temperature isothermal (K)
CpA=64.7081297; %kj/kmol k
CpB=33.43828781; %kj/kmol k
CpC=73.79376; %kj/kmol k
CpD=45.25448885; %kj/kmol k
CpE=35.69429546; %kj/kmol k
Cp_arg= 21.08524537; %kj/kmol k
Cp_meth=48.47777972; %kj/kmol k
Cp_eth=81.59024981; %kj/kmol k
Cp_N2=29.91277101; %kj/kmol k
F_arg=2545.557663; %kmol/h
F_meth=25855.7708; %kmol/h
F_eth=342.8308447; %kmol/h
F_N2=3011.063741; %kmol/h
F_DCE=0.009610; %kmol/h
H1=-117.152*1000; %kj/kmol
H2=-1322.144*1000; %kj/kmol
%Kinetic Constants
k1=6.867*exp(-9260/(1.987*T(1)));
k2=107.3*exp(-10413/(1.987*T(1)));
k3=10.62*exp(-8792/(1.987*T(1)));
k4=39.59*exp(-9473/(1.987*T(1)));
k5=0.000001346*exp(9610/(1.987*T(1)));
k6=0.00000003109*exp(15044/(1.987*T(1)));
k7=0.19;
k8=0.07;
Ft = F(1)+F(2)+F(3)+F(4)+F(5)+ F_arg+ F_meth+ F_N2+ F_DCE; %total molar flowrate
P_O2 = P*( F(2)/Ft); %partial pressure of O2 (kPa)
P_E = P*( F(1)/Ft); %partial pressure of ethylene (kPa)
P_DCE= P*(F_DCE/Ft); %partial pressure of DCE (kPa)
r1=(dcat*((k1*P_O2*P_E)-(k2*P_O2*P_E*P_DCE^k7)))/1000*(1+k5*P_O2+k6*P_E); %kmol/ h m3
r2=(dcat*((k3*P_O2*P_E)-(k4*P_O2*P_E*P_DCE^k8)))/1000*(1+k5*P_O2+k6*P_E); %kmol/ h m3
%flowrate * heat capacity of each component
C_A= F(1)*CpA; %kj/h K
C_B= F(2)*CpB; %kj/h K
C_C= F(3)*CpC; %kj/h K
C_D= F(4)*CpD; %kj/h K
C_E= F(5)*CpE; %kj/h K
C_meth=F_meth*Cp_meth; %kj/h K
C_eth=F_eth*Cp_eth; %kj/h K
C_N2=F_N2*Cp_N2; %kj/h K
C_arg=F_arg*Cp_arg; %kj/h K
%sum of heat capacity flowrates
FC= C_A+C_B+C_C+C_D+C_E+C_meth+C_eth+C_N2+C_arg; %kj/h K
%Mass balance
dFdV(1) = r1+r2;
dFdV(2) = 0.5*r1+3*r2;
dFdV(3) = -r1;
dFdV(4) = -2*r2;
dFdV(5) = -2*r2;
%Energy balance
dTdV(1) = (U*a*(T-Tc)+dcat*(H1*r1+H2*r2))./FC;
end
[V,T] = ode45('Reactor3', [0 100], [523.15]);
plot(V,T)
xlabel('V(m3)');
ylabel('T (K))');
2 件のコメント
Alex Mcaulley
2019 年 3 月 14 日
What values of V, F and T are you using? What is the error are you getting?
回答 (1 件)
Torsten
2019 年 3 月 14 日
編集済み: Torsten
2019 年 3 月 14 日
function dFTdV = Reactor3(V,FT)
%Reactor3 function is used to obtain mass and energy balances across the
%reactor
F = FT(1:5);
T(1) = FT(6);
%Constants
P=2200; %inlet Pressure (kPa)
dcat=850000; %catalyst density g/m3
U=0.17603*3600; %overall heat transfer co (kj/m2 h K)
a=80; %(m^-1)
Tc=433.15; %coolant temperature isothermal (K)
CpA=64.7081297; %kj/kmol k
CpB=33.43828781; %kj/kmol k
CpC=73.79376; %kj/kmol k
CpD=45.25448885; %kj/kmol k
CpE=35.69429546; %kj/kmol k
Cp_arg= 21.08524537; %kj/kmol k
Cp_meth=48.47777972; %kj/kmol k
Cp_eth=81.59024981; %kj/kmol k
Cp_N2=29.91277101; %kj/kmol k
F_arg=2545.557663; %kmol/h
F_meth=25855.7708; %kmol/h
F_eth=342.8308447; %kmol/h
F_N2=3011.063741; %kmol/h
F_DCE=0.009610; %kmol/h
H1=-117.152*1000; %kj/kmol
H2=-1322.144*1000; %kj/kmol
%Kinetic Constants
k1=6.867*exp(-9260/(1.987*T(1)));
k2=107.3*exp(-10413/(1.987*T(1)));
k3=10.62*exp(-8792/(1.987*T(1)));
k4=39.59*exp(-9473/(1.987*T(1)));
k5=0.000001346*exp(9610/(1.987*T(1)));
k6=0.00000003109*exp(15044/(1.987*T(1)));
k7=0.19;
k8=0.07;
Ft = F(1)+F(2)+F(3)+F(4)+F(5)+ F_arg+ F_meth+ F_N2+ F_DCE; %total molar flowrate
P_O2 = P*( F(2)/Ft); %partial pressure of O2 (kPa)
P_E = P*( F(1)/Ft); %partial pressure of ethylene (kPa)
P_DCE= P*(F_DCE/Ft); %partial pressure of DCE (kPa)
r1=(dcat*((k1*P_O2*P_E)-(k2*P_O2*P_E*P_DCE^k7)))/1000*(1+k5*P_O2+k6*P_E); %kmol/ h m3
r2=(dcat*((k3*P_O2*P_E)-(k4*P_O2*P_E*P_DCE^k8)))/1000*(1+k5*P_O2+k6*P_E); %kmol/ h m3
%flowrate * heat capacity of each component
C_A= F(1)*CpA; %kj/h K
C_B= F(2)*CpB; %kj/h K
C_C= F(3)*CpC; %kj/h K
C_D= F(4)*CpD; %kj/h K
C_E= F(5)*CpE; %kj/h K
C_meth=F_meth*Cp_meth; %kj/h K
C_eth=F_eth*Cp_eth; %kj/h K
C_N2=F_N2*Cp_N2; %kj/h K
C_arg=F_arg*Cp_arg; %kj/h K
%sum of heat capacity flowrates
FC= C_A+C_B+C_C+C_D+C_E+C_meth+C_eth+C_N2+C_arg; %kj/h K
%Mass balance
dFdV(1) = r1+r2;
dFdV(2) = 0.5*r1+3*r2;
dFdV(3) = -r1;
dFdV(4) = -2*r2;
dFdV(5) = -2*r2;
%Energy balance
dTdV(1) = (U*a*(T-Tc)+dcat*(H1*r1+H2*r2))./FC;
dFTdV = [dFdV.';dTdV];
end
[V,FT] = ode15s(@Reactor3, [0 100], [13096.019;3156.3;0;32.38;0;523.15]);
plot(V,FT(:,6))
xlabel('V(m3)');
ylabel('T (K))');
2 件のコメント
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!