MATLAB Answers

0

plotting profiles from differential equations

Asked by Liam Gavin on 14 Mar 2019
Latest activity Edited by Torsten
on 14 Mar 2019
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 Comments

What values of V, F and T are you using? What is the error are you getting?
My values of V range from 0-100 m3, and im trying to deduce how the flowrate (F) and teperature (T) change across my volume.
Errors i get when i run the reactor3 file are:
Not enough input arguments.
Error in Reactor3 (line 29)
k1=6.867*exp(-9260/(1.987*T(1)));
Whereas the errors i get when i try to plot the graph include:
Index exceeds array bounds.
Error in Reactor3 (line 29)
k1=6.867*exp(-9260/(1.987*T(1)));
Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in energy_balance_solver (line 1)
[V,T] = ode45('Reactor3', [0 100], [523.15]);

Sign in to comment.

1 Answer

Answer by Torsten
on 14 Mar 2019
Edited by Torsten
on 14 Mar 2019

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 Comments

Hi, I really appreciate the help!
However the variables F and T are differrent. I know that all differential equations are dependant on T and most of them are dependant on F.
Is there anyway that i could plot how just T changes with V? Not FT?
T depends on F ; thus you have to solve all 6 equations simultaneoulsy.
The plot shows how T changes with V. I don't understand what you mean by your last question.

Sign in to comment.