Why am I getting ode45 error?

Question

Rishita Srivastava 2023-10-27

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2039351-why-am-i-getting-ode45-error

评论： Walter Roberson 2023-10-27

在 MATLAB Online 中打开

My code is shown below. There are 2 files, one is a function file and the other is the execution file.

The execution file on running keeps returning an ode45 error.

function dW=PCDDtrial(~,f)
% A=PCDD1; B=O2; C=H2O; D=CO2; E=HCl
fA=f(1);fB=f(2);fC=f(3);fD=f(4);fE=f(5);y=f(6);T=f(7);Ta=f(8);
%% Constant Parameter Values
R=8.314; %Universal Gas Constant (J/mol.K)
To=500; %Inlet Temperature (K) 
Tao= 900; % heating temp (K)
Po=1.5; %Inlet Pressure (bar)  
N=1500; %Number of Tubes      %% to be changed
rho_gas=1.03E-02; %Inlet Density of Gas (kg/m3)   
rho_cat=1200; %Catalyst Apparent Density (kg/m3)  
rho_bulk=600; %Catalyst Bulk Density (kg/m3)     
phi=0.5; %Porosity   
Dt=0.05; %Tube Diameter (m)  
Dp=0.005; %Catalyst Diameter (m)      0.15-0.18mm
Ac=Dt^2*pi/4; %Tube Cross-sectional Area (m2)
Ua=4500; %Heat Transfer Coefficient for Area Per Unit Volume (W/(m3.K))      %% to be researched
%% Feed conditions 
%Molar Flowrate Per Tube
fAo=0.8772/N; %Inlet Molar Flowrate of PCDD1 Per Tube (mol/s)
fBo=12752.6/N; %Inlet Molar Flowrate of O2 Per Tube (mol/s)
fCo=0/N; %Inlet Molar Flowrate of H2O Per Tube (mol/s)
fDo= 0/N; %Inlet Molar Flowrate of CO2 Per Tube (mol/s)
fEo= 0/N; %Inlet Molar Flowrate of HCl Per Tube (mol/s)
fTo=fAo+fBo+fCo+fDo+fEo; %Inlet Total Molar Flowrate Per Tube (mol/s)
fT=fA+fB+fC+fD+fE; %Outlet Total Molar Flowrate Per Tube (mol/s)
%Inlet Mass Flowrate Per Tube
fAom=(fAo*356.4)/1000; %Inlet Mass Flowrate of PCDD1 Per Tube (kg/s)
fBom=(fBo*32)/1000; %Inlet Mass Flowrate of O2 Per Tube (kg/s)
fCom=(fDo*18)/1000; %Inlet Mass Flowrate of H2O Per Tube (kg/s)
fDom=(fDo*44)/1000; %Inlet Mass Flowrate of CO2 Per Tube (kg/s)
fEom=(fEo*36)/1000; %Inlet Mass Flowrate of HCl Per Tube (kg/s)
fTom=fAom+fBom+fCom+fDom+fEom;%Inlet Total Mass Flowrate Per Tube (kg/s)
%Viscosity Calculation 
viscA= 6.391E-14*To^4 -2.147E-10*To^3 +2.676E-07*To^2 -0.0001471*To + 0.03032; %Viscosity of PCDD1 (Pa.s)
viscB=(1.101E-06*To^0.5634)/(1+(96.3/To)); %Viscosity of O2 (Pa.s) ok
viscC=(1.7096E-06*To^1.1146); %Viscosity of H2O (Pa.s) ok
viscD=(2.148E-06*To^0.46)/(1+(290/To)); %Viscosity of CO2 (Pa.s) ok
viscE=(4.924E-07*To^0.6702)/(1+(157.7/To)); %Viscosity of HCl (Pa.s) ok
visc_avg=(viscA*(fAo/fTo)+viscB*(fBo/fTo)+viscC*(fCo/fTo)+viscD*(fDo/fTo)+viscE*(fEo/fTo)); %Average Viscosity (Pa.s)
%Alpha Calculation
G=fTom/(Ac); %Gas Mass Velocity (kg/m2.s)
Beta=((G*(1-phi))/(rho_gas*1*Dp*phi^3))*((150*(1-phi)*visc_avg)/Dp+1.75*G); %Pressure Per Length (Pa/m)
Alpha=(2*Beta)/(Ac*rho_cat*(1-phi)*Po*10^5);%Pressure Term (1/kg)
%Concentration Per Tube
CTo=Po/(R*To); %Inlet Concentration (mol/m3) CTo=Po/(R*To)
vo=fTo/CTo; %Volumetric Flowrate Per Tube (m3/s) fTo/CTo
%Outlet Concentration Calculation
CA=CTo*(To/T)*y*(fA/fT); %Outlet Concentration of PCDD1 (mol/m3)
CB=CTo*(To/T)*y*(fB/fT); %Outlet Concentration of O2 (mol/m3)
CC=CTo*(To/T)*y*(fC/fT); %Outlet Concentration of H2O (mol/m3)
CD=CTo*(To/T)*y*(fD/fT); %Outlet Concentration of CO2 (mol/m3)
CE=CTo*(To/T)*y*(fE/fT); %Outlet Concentration of HCl (mol/m3)
%% Arrhenius equations for rate constants
k1=7.8E02*exp(-15000/(R*T)); %k for PCDD1 gas (ml/g.s)
%% Reaction kinetics
R1= (k1*CA)/1000; %(mol/kg.s)   
%% Net reaction rate
RA=-R1;
RB=-12.5*R1;
RC=3*R1;
RD= 12*R1;
RE= R1;
%% Heat calculations  
HR1o=-10280.92; %Heat of Reaction 1 at 298.15 K (kJ/mol)
%Heat Capacity  %% to be changed
Tx=298.15; %Reference Temperature (K)
CpA=(-0.2962*To^2+691.4*To+2.724E04)/1000; %Heat Capacity of PCDD1 (J/mol.K)
CpB=(175430-6152.3*To+113.92*To^2-0.92382*To^3+0.0027963*To^4)/1000; %Heat Capacity of O2 (J/mol.K)
CpC=(276370-2090.1*To+8.125*To^2-0.014116*To^3+9.3701E-06*To^4)/1000; %Heat Capacity of H2O (J/mol.K)
CpD=(-8304300+104370*To-433.33*To^2+0.60052*To^3)/1000; %Heat Capacity of CO2 (J/mol.K)
CpE=(47300+90*To)/1000; %Heat Capacity of HCl (J/mol.K)
HR1= (HR1o+((6*CpC+24*CpD+2*CpE-25*CpB-2*CpA)*(T-Tx)/1000))*1000; %(J/mol)
%% ODE Equation
dW(1)=RA;dW(2)=RB;dW(3)=RC;dW(4)=RD; dW(5)=RE;
dW(6)=-(Alpha*fT*T)/(2*y*fTo*To);
dW(7)=((HR1*R1))/((fA*CpA+fB*CpB+fC*CpC+fD*CpD+fE*CpE)); %K/kg 
dW=[dW(1);dW(2);dW(3);dW(4);dW(5);dW(6);dW(7);];
%%Execution file 
Wi=0:0.1:25; %Range of W
%Initial Condition
To=500;N=1500;R=8.314;rho_cat=1200;
fAo=0.8772/N; % (mol/s)
fBo=12752.6/N; %(mol/s)
fCo=0/N; %(mol/s)
fDo= 0/N; %(mol/s)
fEo= 0/N; % (mol/s)
Fi=[fAo;fBo;fCo;fDo;fEo;1;To]; %Initial Conditions Tao
[W,f]=ode45('PCDDtrial',Wi,Fi);
fa1=f(:,1)*N; %mol/s
fb1=f(:,2)*N;
fc1=f(:,3)*N;
fd1=f(:,4)*N;
fe1=f(:,5)*N;
y=f(:,6);
T1=f(:,7);
%Ta1=f(:,8);
z=[fa1 fb1 fc1 fd1 fe1 y T1]; %Results Ta1
XA=((fAo*N-fa1)/(fAo*N))*100; %Conversion of PCDD1
%Plot Graphs
subplot(3,3,1)
plot(W,fa1)
xlabel('Weight of Catalyst Per Tube , W (kg)')
ylabel('Molar Flowrate of PCDD1 (mol/s)')
ytickformat('%.5f');
ax.YAxis.Exponent = 0;
subplot(3,3,2)
plot(W,fb1)
xlabel('Weight of Catalyst Per Tube , W (kg)')
ylabel('Molar Flowrate of O2 (mol/s)')
ytickformat('%.5f');
ax.YAxis.Exponent = 0;
subplot(3,3,3)
plot(W,fc1)
xlabel('Weight of Catalyst Per Tube , W (kg)')
ylabel('Molar Flowrate of H2O (mol/s)')
ytickformat('%.5f');
ax.YAxis.Exponent = 0;
subplot(3,3,4)
plot(W,fd1)
xlabel('Weight of Catalyst Per Tube , W (kg)')
ylabel('Molar Flowrate of CO2 (mol/s)')
subplot(3,3,5)
plot(W,fe1)
xlabel('Weight of Catalyst Per Tube , W (kg)')
ylabel('Molar Flowrate of HCl (mol/s)')
ytickformat('%.5f');
ax.YAxis.Exponent = 0;
subplot(3,3,6)
plot(W,y)
xlabel('Weight of Catalyst Per Tube , W (kg)')
ylabel('Pressure Ratio, y')
subplot(3,3,7)
plot(W,T1) %W,Ta1
xlabel('Weight of Catalyst Per Tube , W (kg)')
ylabel('Temperature (K)')
legend('Reactor Temperature','Cooling Water Temperature')
subplot(3,3,8)
plot(W,XA)
xlabel('Weight of Catalyst Per Tube , W (kg)')
ylabel('Conversion of C2H4, X')
for i1=1:8
subplot(3,3,i1)
grid on
grid minor
end

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Walter Roberson 2023-10-27

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2039351-why-am-i-getting-ode45-error#answer_1341666

编辑：Walter Roberson 2023-10-27

在 MATLAB Online 中打开

I am still testing the below; some minor changes need to be made to some of your plotting code.

%%Execution file 
Wi=0:0.1:25; %Range of W
%Initial Condition
To=500;N=1500;R=8.314;rho_cat=1200;
fAo=0.8772/N; % (mol/s)
fBo=12752.6/N; %(mol/s)
fCo=0/N; %(mol/s)
fDo= 0/N; %(mol/s)
fEo= 0/N; % (mol/s)
Fi=[fAo;fBo;fCo;fDo;fEo;1;To]; %Initial Conditions Tao
[W,f]=ode45(@PCDDtrial,Wi,Fi);
fa1=f(:,1)*N; %mol/s
fb1=f(:,2)*N;
fc1=f(:,3)*N;
fd1=f(:,4)*N;
fe1=f(:,5)*N;
y=f(:,6);
T1=f(:,7);
%Ta1=f(:,8);
z=[fa1 fb1 fc1 fd1 fe1 y T1]; %Results Ta1
XA=((fAo*N-fa1)/(fAo*N))*100; %Conversion of PCDD1
%Plot Graphs
subplot(3,3,1)
plot(W,fa1)
xlabel('Weight of Catalyst Per Tube , W (kg)')
ylabel('Molar Flowrate of PCDD1 (mol/s)')
ytickformat('%.5f');
ax.YAxis.Exponent = 0;
subplot(3,3,2)
plot(W,fb1)
xlabel('Weight of Catalyst Per Tube , W (kg)')
ylabel('Molar Flowrate of O2 (mol/s)')
ytickformat('%.5f');
ax.YAxis.Exponent = 0;
subplot(3,3,3)
plot(W,fc1)
xlabel('Weight of Catalyst Per Tube , W (kg)')
ylabel('Molar Flowrate of H2O (mol/s)')
ytickformat('%.5f');
ax.YAxis.Exponent = 0;
subplot(3,3,4)
plot(W,fd1)
xlabel('Weight of Catalyst Per Tube , W (kg)')
ylabel('Molar Flowrate of CO2 (mol/s)')
subplot(3,3,5)
plot(W,fe1)
xlabel('Weight of Catalyst Per Tube , W (kg)')
ylabel('Molar Flowrate of HCl (mol/s)')
ytickformat('%.5f');
ax.YAxis.Exponent = 0;
subplot(3,3,6)
plot(W,y)
xlabel('Weight of Catalyst Per Tube , W (kg)')
ylabel('Pressure Ratio, y')
subplot(3,3,7)
plot(W,T1) %W,Ta1
xlabel('Weight of Catalyst Per Tube , W (kg)')
ylabel('Temperature (K)')
legend('Reactor Temperature','Cooling Water Temperature')
subplot(3,3,8)
plot(W,XA)
xlabel('Weight of Catalyst Per Tube , W (kg)')
ylabel('Conversion of C2H4, X')
for i1=1:8
subplot(3,3,i1)
grid on
grid minor
end
function dW=PCDDtrial(~,f)
% A=PCDD1; B=O2; C=H2O; D=CO2; E=HCl
fA=f(1);fB=f(2);fC=f(3);fD=f(4);fE=f(5);y=f(6);T=f(7);
%Ta=f(8);
%% Constant Parameter Values
R=8.314; %Universal Gas Constant (J/mol.K)
To=500; %Inlet Temperature (K) 
Tao= 900; % heating temp (K)
Po=1.5; %Inlet Pressure (bar)  
N=1500; %Number of Tubes      %% to be changed
rho_gas=1.03E-02; %Inlet Density of Gas (kg/m3)   
rho_cat=1200; %Catalyst Apparent Density (kg/m3)  
rho_bulk=600; %Catalyst Bulk Density (kg/m3)     
phi=0.5; %Porosity   
Dt=0.05; %Tube Diameter (m)  
Dp=0.005; %Catalyst Diameter (m)      0.15-0.18mm
Ac=Dt^2*pi/4; %Tube Cross-sectional Area (m2)
Ua=4500; %Heat Transfer Coefficient for Area Per Unit Volume (W/(m3.K))      %% to be researched
%% Feed conditions 
%Molar Flowrate Per Tube
fAo=0.8772/N; %Inlet Molar Flowrate of PCDD1 Per Tube (mol/s)
fBo=12752.6/N; %Inlet Molar Flowrate of O2 Per Tube (mol/s)
fCo=0/N; %Inlet Molar Flowrate of H2O Per Tube (mol/s)
fDo= 0/N; %Inlet Molar Flowrate of CO2 Per Tube (mol/s)
fEo= 0/N; %Inlet Molar Flowrate of HCl Per Tube (mol/s)
fTo=fAo+fBo+fCo+fDo+fEo; %Inlet Total Molar Flowrate Per Tube (mol/s)
fT=fA+fB+fC+fD+fE; %Outlet Total Molar Flowrate Per Tube (mol/s)
%Inlet Mass Flowrate Per Tube
fAom=(fAo*356.4)/1000; %Inlet Mass Flowrate of PCDD1 Per Tube (kg/s)
fBom=(fBo*32)/1000; %Inlet Mass Flowrate of O2 Per Tube (kg/s)
fCom=(fDo*18)/1000; %Inlet Mass Flowrate of H2O Per Tube (kg/s)
fDom=(fDo*44)/1000; %Inlet Mass Flowrate of CO2 Per Tube (kg/s)
fEom=(fEo*36)/1000; %Inlet Mass Flowrate of HCl Per Tube (kg/s)
fTom=fAom+fBom+fCom+fDom+fEom;%Inlet Total Mass Flowrate Per Tube (kg/s)
%Viscosity Calculation 
viscA= 6.391E-14*To^4 -2.147E-10*To^3 +2.676E-07*To^2 -0.0001471*To + 0.03032; %Viscosity of PCDD1 (Pa.s)
viscB=(1.101E-06*To^0.5634)/(1+(96.3/To)); %Viscosity of O2 (Pa.s) ok
viscC=(1.7096E-06*To^1.1146); %Viscosity of H2O (Pa.s) ok
viscD=(2.148E-06*To^0.46)/(1+(290/To)); %Viscosity of CO2 (Pa.s) ok
viscE=(4.924E-07*To^0.6702)/(1+(157.7/To)); %Viscosity of HCl (Pa.s) ok
visc_avg=(viscA*(fAo/fTo)+viscB*(fBo/fTo)+viscC*(fCo/fTo)+viscD*(fDo/fTo)+viscE*(fEo/fTo)); %Average Viscosity (Pa.s)
%Alpha Calculation
G=fTom/(Ac); %Gas Mass Velocity (kg/m2.s)
Beta=((G*(1-phi))/(rho_gas*1*Dp*phi^3))*((150*(1-phi)*visc_avg)/Dp+1.75*G); %Pressure Per Length (Pa/m)
Alpha=(2*Beta)/(Ac*rho_cat*(1-phi)*Po*10^5);%Pressure Term (1/kg)
%Concentration Per Tube
CTo=Po/(R*To); %Inlet Concentration (mol/m3) CTo=Po/(R*To)
vo=fTo/CTo; %Volumetric Flowrate Per Tube (m3/s) fTo/CTo
%Outlet Concentration Calculation
CA=CTo*(To/T)*y*(fA/fT); %Outlet Concentration of PCDD1 (mol/m3)
CB=CTo*(To/T)*y*(fB/fT); %Outlet Concentration of O2 (mol/m3)
CC=CTo*(To/T)*y*(fC/fT); %Outlet Concentration of H2O (mol/m3)
CD=CTo*(To/T)*y*(fD/fT); %Outlet Concentration of CO2 (mol/m3)
CE=CTo*(To/T)*y*(fE/fT); %Outlet Concentration of HCl (mol/m3)
%% Arrhenius equations for rate constants
k1=7.8E02*exp(-15000/(R*T)); %k for PCDD1 gas (ml/g.s)
%% Reaction kinetics
R1= (k1*CA)/1000; %(mol/kg.s)   
%% Net reaction rate
RA=-R1;
RB=-12.5*R1;
RC=3*R1;
RD= 12*R1;
RE= R1;
%% Heat calculations  
HR1o=-10280.92; %Heat of Reaction 1 at 298.15 K (kJ/mol)
%Heat Capacity  %% to be changed
Tx=298.15; %Reference Temperature (K)
CpA=(-0.2962*To^2+691.4*To+2.724E04)/1000; %Heat Capacity of PCDD1 (J/mol.K)
CpB=(175430-6152.3*To+113.92*To^2-0.92382*To^3+0.0027963*To^4)/1000; %Heat Capacity of O2 (J/mol.K)
CpC=(276370-2090.1*To+8.125*To^2-0.014116*To^3+9.3701E-06*To^4)/1000; %Heat Capacity of H2O (J/mol.K)
CpD=(-8304300+104370*To-433.33*To^2+0.60052*To^3)/1000; %Heat Capacity of CO2 (J/mol.K)
CpE=(47300+90*To)/1000; %Heat Capacity of HCl (J/mol.K)
HR1= (HR1o+((6*CpC+24*CpD+2*CpE-25*CpB-2*CpA)*(T-Tx)/1000))*1000; %(J/mol)
%% ODE Equation
dW(1)=RA;dW(2)=RB;dW(3)=RC;dW(4)=RD; dW(5)=RE;
dW(6)=-(Alpha*fT*T)/(2*y*fTo*To);
dW(7)=((HR1*R1))/((fA*CpA+fB*CpB+fC*CpC+fD*CpD+fE*CpE)); %K/kg 
dW=[dW(1);dW(2);dW(3);dW(4);dW(5);dW(6);dW(7);];
end

2 个评论
显示无隐藏无

Walter Roberson 2023-10-27

在 MATLAB Online 中打开

If you try to run this, it will take a long long long time.

Essentially you have a stiff system.

However.. if you switch to ode15s or ode23s then at a time earlier than 2e-5 it will fail integration.

I am currently running it with ode45 and time limit 2e-5 and it has taken an hour or so already.

As well I had to make some minor code changes; the current version is

%%Execution file 
Wi=0:0.1:25; %Range of W
%Initial Condition
To=500;N=1500;R=8.314;rho_cat=1200;
fAo=0.8772/N; % (mol/s)
fBo=12752.6/N; %(mol/s)
fCo=0/N; %(mol/s)
fDo= 0/N; %(mol/s)
fEo= 0/N; % (mol/s)
Fi=[fAo;fBo;fCo;fDo;fEo;1;To]; %Initial Conditions Tao
Wi = [0 1e-4];
[W,f]=ode45(@PCDDtrial,Wi,Fi);
fa1=f(:,1)*N; %mol/s
fb1=f(:,2)*N;
fc1=f(:,3)*N;
fd1=f(:,4)*N;
fe1=f(:,5)*N;
y=f(:,6);
T1=f(:,7);
%Ta1=f(:,8);
z=[fa1 fb1 fc1 fd1 fe1 y T1]; %Results Ta1
XA=((fAo*N-fa1)/(fAo*N))*100; %Conversion of PCDD1
%Plot Graphs
subplot(3,3,1)
plot(W,fa1)
xlabel('Weight of Catalyst Per Tube , W (kg)')
ylabel('Molar Flowrate of PCDD1 (mol/s)')
ytickformat('%.5f');
ax = gca;
ax.YAxis.Exponent = 0;
subplot(3,3,2)
plot(W,fb1)
xlabel('Weight of Catalyst Per Tube , W (kg)')
ylabel('Molar Flowrate of O2 (mol/s)')
ytickformat('%.5f');
ax = gca;
ax.YAxis.Exponent = 0;
subplot(3,3,3)
plot(W,fc1)
xlabel('Weight of Catalyst Per Tube , W (kg)')
ylabel('Molar Flowrate of H2O (mol/s)')
ytickformat('%.5f');
ax = gca;
ax.YAxis.Exponent = 0;
subplot(3,3,4)
plot(W,fd1)
xlabel('Weight of Catalyst Per Tube , W (kg)')
ylabel('Molar Flowrate of CO2 (mol/s)')
subplot(3,3,5)
plot(W,fe1)
xlabel('Weight of Catalyst Per Tube , W (kg)')
ylabel('Molar Flowrate of HCl (mol/s)')
ytickformat('%.5f');
ax = gca;
ax.YAxis.Exponent = 0;
subplot(3,3,6)
plot(W,y)
xlabel('Weight of Catalyst Per Tube , W (kg)')
ylabel('Pressure Ratio, y')
subplot(3,3,7)
plot(W,T1) %W,Ta1
xlabel('Weight of Catalyst Per Tube , W (kg)')
ylabel('Temperature (K)')
%legend('Reactor Temperature','Cooling Water Temperature')
subplot(3,3,8)
plot(W,XA)
xlabel('Weight of Catalyst Per Tube , W (kg)')
ylabel('Conversion of C2H4, X')
for i1=1:8
subplot(3,3,i1)
grid on
grid minor
end
function dW=PCDDtrial(~,f)
% A=PCDD1; B=O2; C=H2O; D=CO2; E=HCl
fA=f(1);fB=f(2);fC=f(3);fD=f(4);fE=f(5);y=f(6);T=f(7);
%Ta=f(8);
%% Constant Parameter Values
R=8.314; %Universal Gas Constant (J/mol.K)
To=500; %Inlet Temperature (K) 
Tao= 900; % heating temp (K)
Po=1.5; %Inlet Pressure (bar)  
N=1500; %Number of Tubes      %% to be changed
rho_gas=1.03E-02; %Inlet Density of Gas (kg/m3)   
rho_cat=1200; %Catalyst Apparent Density (kg/m3)  
rho_bulk=600; %Catalyst Bulk Density (kg/m3)     
phi=0.5; %Porosity   
Dt=0.05; %Tube Diameter (m)  
Dp=0.005; %Catalyst Diameter (m)      0.15-0.18mm
Ac=Dt^2*pi/4; %Tube Cross-sectional Area (m2)
Ua=4500; %Heat Transfer Coefficient for Area Per Unit Volume (W/(m3.K))      %% to be researched
%% Feed conditions 
%Molar Flowrate Per Tube
fAo=0.8772/N; %Inlet Molar Flowrate of PCDD1 Per Tube (mol/s)
fBo=12752.6/N; %Inlet Molar Flowrate of O2 Per Tube (mol/s)
fCo=0/N; %Inlet Molar Flowrate of H2O Per Tube (mol/s)
fDo= 0/N; %Inlet Molar Flowrate of CO2 Per Tube (mol/s)
fEo= 0/N; %Inlet Molar Flowrate of HCl Per Tube (mol/s)
fTo=fAo+fBo+fCo+fDo+fEo; %Inlet Total Molar Flowrate Per Tube (mol/s)
fT=fA+fB+fC+fD+fE; %Outlet Total Molar Flowrate Per Tube (mol/s)
%Inlet Mass Flowrate Per Tube
fAom=(fAo*356.4)/1000; %Inlet Mass Flowrate of PCDD1 Per Tube (kg/s)
fBom=(fBo*32)/1000; %Inlet Mass Flowrate of O2 Per Tube (kg/s)
fCom=(fDo*18)/1000; %Inlet Mass Flowrate of H2O Per Tube (kg/s)
fDom=(fDo*44)/1000; %Inlet Mass Flowrate of CO2 Per Tube (kg/s)
fEom=(fEo*36)/1000; %Inlet Mass Flowrate of HCl Per Tube (kg/s)
fTom=fAom+fBom+fCom+fDom+fEom;%Inlet Total Mass Flowrate Per Tube (kg/s)
%Viscosity Calculation 
viscA= 6.391E-14*To^4 -2.147E-10*To^3 +2.676E-07*To^2 -0.0001471*To + 0.03032; %Viscosity of PCDD1 (Pa.s)
viscB=(1.101E-06*To^0.5634)/(1+(96.3/To)); %Viscosity of O2 (Pa.s) ok
viscC=(1.7096E-06*To^1.1146); %Viscosity of H2O (Pa.s) ok
viscD=(2.148E-06*To^0.46)/(1+(290/To)); %Viscosity of CO2 (Pa.s) ok
viscE=(4.924E-07*To^0.6702)/(1+(157.7/To)); %Viscosity of HCl (Pa.s) ok
visc_avg=(viscA*(fAo/fTo)+viscB*(fBo/fTo)+viscC*(fCo/fTo)+viscD*(fDo/fTo)+viscE*(fEo/fTo)); %Average Viscosity (Pa.s)
%Alpha Calculation
G=fTom/(Ac); %Gas Mass Velocity (kg/m2.s)
Beta=((G*(1-phi))/(rho_gas*1*Dp*phi^3))*((150*(1-phi)*visc_avg)/Dp+1.75*G); %Pressure Per Length (Pa/m)
Alpha=(2*Beta)/(Ac*rho_cat*(1-phi)*Po*10^5);%Pressure Term (1/kg)
%Concentration Per Tube
CTo=Po/(R*To); %Inlet Concentration (mol/m3) CTo=Po/(R*To)
vo=fTo/CTo; %Volumetric Flowrate Per Tube (m3/s) fTo/CTo
%Outlet Concentration Calculation
CA=CTo*(To/T)*y*(fA/fT); %Outlet Concentration of PCDD1 (mol/m3)
CB=CTo*(To/T)*y*(fB/fT); %Outlet Concentration of O2 (mol/m3)
CC=CTo*(To/T)*y*(fC/fT); %Outlet Concentration of H2O (mol/m3)
CD=CTo*(To/T)*y*(fD/fT); %Outlet Concentration of CO2 (mol/m3)
CE=CTo*(To/T)*y*(fE/fT); %Outlet Concentration of HCl (mol/m3)
%% Arrhenius equations for rate constants
k1=7.8E02*exp(-15000/(R*T)); %k for PCDD1 gas (ml/g.s)
%% Reaction kinetics
R1= (k1*CA)/1000; %(mol/kg.s)   
%% Net reaction rate
RA=-R1;
RB=-12.5*R1;
RC=3*R1;
RD= 12*R1;
RE= R1;
%% Heat calculations  
HR1o=-10280.92; %Heat of Reaction 1 at 298.15 K (kJ/mol)
%Heat Capacity  %% to be changed
Tx=298.15; %Reference Temperature (K)
CpA=(-0.2962*To^2+691.4*To+2.724E04)/1000; %Heat Capacity of PCDD1 (J/mol.K)
CpB=(175430-6152.3*To+113.92*To^2-0.92382*To^3+0.0027963*To^4)/1000; %Heat Capacity of O2 (J/mol.K)
CpC=(276370-2090.1*To+8.125*To^2-0.014116*To^3+9.3701E-06*To^4)/1000; %Heat Capacity of H2O (J/mol.K)
CpD=(-8304300+104370*To-433.33*To^2+0.60052*To^3)/1000; %Heat Capacity of CO2 (J/mol.K)
CpE=(47300+90*To)/1000; %Heat Capacity of HCl (J/mol.K)
HR1= (HR1o+((6*CpC+24*CpD+2*CpE-25*CpB-2*CpA)*(T-Tx)/1000))*1000; %(J/mol)
%% ODE Equation
dW(1)=RA;dW(2)=RB;dW(3)=RC;dW(4)=RD; dW(5)=RE;
dW(6)=-(Alpha*fT*T)/(2*y*fTo*To);
dW(7)=((HR1*R1))/((fA*CpA+fB*CpB+fC*CpC+fD*CpD+fE*CpE)); %K/kg 
dW=[dW(1);dW(2);dW(3);dW(4);dW(5);dW(6);dW(7);];
end