Error using odearguments ; Error in ode15s (line 153); Error in anand_shrivastan (line 625) ;below is the code

4 views (last 30 days)
% Defining Global Variables
global E
global DiffMix
global Deff
global Tia
global i
global Ru
global Pres
global RR
global EffDia
global EffMol
global EffKBEp
global u
global L
global dx
global Concc
global h
%%%%%%%%%%%%%%%%%%%%%% MODELING PARAMETERS %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Universal Gas Constant (J/K/mol)
% Ref: From Fundamentals of Heat and Mass Transfer by Frank Incorporation
Ru=8.314;
% Defining values of coefficients from Chemkin thermodynamics data which
% needs to be incorporated for the calculation of CP and Enthalpy
% The values given below can be used when the temperatures are from 300 to
% 1000K for the respective species
a11=3.376541; a12=0.0012530634; a13=-0.000003302750;
a14=0.000000005217810; a15=-0.000000000002446262; a16=9817.961;
% Coefficients of NO
a21=3.262451; a22=0.0015119409; a23=-0.000003881755;
a24=0.000000005581944; a25=-0.000000000002474951; a26=-14310.539;
% Coefficients of CO
a31=3.298677; a32=0.0014082404; a33=-0.000003963222;
a34=0.000000005641515; a35=-0.000000000002444854; a36=-1020.8999;
% Coefficients of N2
a41=2.275724; a42=0.009922072; a43=-0.000010409113;
a44=0.000000006866686; a45=-0.000000000002117280; a46=-48373.14;
% Coefficients of CO2
a51=2.500000; a52=0.000000000; a53=0.0000000000000;
a54=0.0000000000000; a55=0.0000000000000; a56=-745.3750;
% Coefficients of He
% Ref: Granger et al.2001_TPDstudiesofNO-CORn_Topicsofcatalysis_V16-17_394.
% Granger et al.1998_KineticsofNO-CORnoverRh/Al2O3_J.Catal.V175_194-203
% Density of particle density (kg/m^3)
DenP = 1000;
% Thermal Conductivity of the pellets based on Prater Number (W/m/k)
Kp=0.22;
% Thermal Conductivity of the gas (NO+CO+He) at initial conditions (W/m/K)
Kg=0.14884;
% Constant Specific Heat for pellets at Constant Pressure (J/KgK)
Cpp = 813.25;
% Length of Pellet bed Reactor (m) -> (6 cm)
L=0.06;
% Length of Catalyst bed used for modeling (m) -> (3 mm)
Lcb = 0.003;
% Diameter of Packed bed reactor (m) -> (1.2 cm)
D=0.012;
% Space velocity (per hour)
SV = 93000;
% Diameter of the pellets (m)
db = 8*10^-5;
% Pore Diameter of the pellets (m)
dp = 1.7*10^-9;
% Pressure of the packed bed reactor (N/m^2) (or) 1 atm
Pres=101325;
% Geometric Surface area per unit volume (m^2/m^3)
Gca= 1;
% Inlet Velocity for packed bed calculated from space velocity (m/s)
u = (SV/3600)*L;
% Porosity of the packed bed reactor from Muller's Expression
E = 0.379 + (0.078/((D/db)-1.80));
% Since, we don't have data for the tortuous path. We can use Ho and
% Strider expression for finding out tortuosity factor
TorF = 1- (0.5*log(E));
% Since, we don't have data for the tortuous path. We can use Ho and
% Strieder expression for finding out Knusden tortuosity factor
KnTorF = (9/8)-(0.5*log(E))+(((13/8)-(9/8))*((E)^0.4));
% Tortuosity of the flow in packed bed reactor
Tor = sqrt (TorF);
% dynamic viscosity of the fluid in (Kg/m.s)
DyVis= 2.1974*10^-5;
% Values for calculating molecular diffusivity
% Ref: Diffusion Mass Transfer in Fluid System by Cussler
% Collision diameter in angstroms
Dia(1)=3.492; % Diameter of NO
Dia(2)=3.690; % Diameter of CO
Dia(3)=3.798; % Diameter of N2
Dia(4)=3.941; % Diameter of CO2
Dia(5)=2.551; % Diameter of He
% Ratio of Epsilon by kb in degree Kelvin
EpKB(1)=116.7; % Epsilon by Kb for NO
EpKB(2)=91.7; % Epsilon by Kb for CO
EpKB(3)=71.4; % Epsilon by Kb for N2
EpKB(4)=195.2; % Epsilon by Kb for CO2
EpKB(5)=10.22; % Epsilon by Kb for He
% Molecular weight of the species
Mol(1)=30.0061; % Molecular Weight of NO (Kg/Kmol)
Mol(2)=28.01; % Molecular Weight of CO (Kg/Kmol)
Mol(3)=28.01348; % Molecular Weight of N2 (Kg/kmol)
Mol(4)=44.0095; % Molecular Weight of CO2 (Kg/Kmol)
Mol(5)=4.003; % Molecular Weight of He (Kg/Kmol)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% DISCRETIZATION OF THE REACTOR %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Welcome to 1-D Packed Bed Modeling for NO reduction reaction by CO over rhodium/alumina catalysts');
Welcome to 1-D Packed Bed Modeling for NO reduction reaction by CO over rhodium/alumina catalysts
disp('Length of the catalyst bed in 0.003 m');
Length of the catalyst bed in 0.003 m
%h=input('Enter the Discretization number:');
h = 10;
while h > 10
disp('For this modeling the length of the catalyst bed is 0.003 m and so,please discritize accordingly to reduce computational time');
h=input('Enter the Discretization number:');
if isempty(h)
h=input('Please enter the Discretization number less than 10:');
end
end
% Discretization Length of the reactor (m)
dx = Lcb/(h-1);
%%%%%%%% INITIALIZE MOLE FRACTIONS AND TEMPERATURE OF THE REACTOR %%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initial Mole fractions are equal to the initial partial pressures of
% NO and CO provided as the reactor operates at 1 atm pressure condition
% Partial Pressure of the reactants (atm)
P(1)=0.005; % Inlet NO partial pressure
P(2)=0.005; % Inlet CO partial pressure
P(3)=0; % Inlet N2 partial pressure
P(4)=0; % Inlet CO2 partial pressure
P(5)=0.99; % Inlet He partial pressure
% Initial Mole fractions of the reactants to check Biot Number limitations
X(1,1)=0.005; % Inlet NO Mole fraction
X(1,2)=0.005; % Inlet CO Mole fraction
X(1,3)=0; % Inlet N2 Mole fraction
X(1,4)=0; % Inlet CO2 Mole fraction
X(1,5)=0.99; % Inlet He Mole fraction
% Initial Temperature across the reactor bed (K)
Tin=303;
%%%%%%%%%%%%%%%%% CALCULATION OF MOLECULAR DIFFUSIVITY %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialization Molecular mass mixture (Kg/Kmol)
MolMassMix=0;
MolMass = zeros(5);
for j=1:5
MolMass(j)=0;
for k=1:5
% Calculation of Effective Collision Diameter
EffDia(j,k)=0.5*(Dia(j)+Dia(k));
% Calculation of Effective Molecular Weight
EffMol(j,k)=((Mol(j)+Mol(k))/(Mol(j)*Mol(k)))^0.5;
if j~=k
% Effective Energy Calculation
EffKBEp(j,k)=(1/((EpKB(j)*EpKB(k))^0.5));
end
end
% Molecular Mass of individual Species in 1 kmol of mixture (Kg/Kmol)
MolMass(j)= X(1,j)*Mol(j);
% Calculation of Molecular Mass of Mixture (Kg/Kmol)
MolMassMix=MolMassMix+MolMass(j);
end
% Gas Constant of Mixture (J/Kg/K)
Rmass=(Ru*1000)/MolMassMix;
% Density of the fluid (Kg/m^3)
DenF=Pres/(Rmass*Tin);
% Initializing inlet concentrations and mass fractions for all species
% at the initial condition
Cin = zeros(1,5);
Y = zeros(1,5);
KBTEp = zeros(5,5);
Ohm = zeros(5,5);
Diff = zeros(5,5);
DiffMix = zeros(5);
Dm = zeros(5);
DKn = zeros(5);
for j=1:5
% Inlet Concentration of Individual Species at first node (mol/m^3)
Cin(1,j)=(DenF/MolMassMix)*X(1,j)*1000;
% Mass Fraction of Individual Species at first node
Y(1,j)=(X(1,j)*Mol(j))/MolMassMix;
end
% Calculating Binary Molecular diffusivity of the species
for j=1:5
for k=1:5
% Inverse Energies
KBTEp(j,k)=(EffKBEp(j,k)*Tin);
% Calculation of Omega
if KBTEp(j,k)<5
Ohm(j,k)=1.4803*(KBTEp(j,k)^-0.397);
else
Ohm(j,k)=1.0765*(KBTEp(j,k)^-0.16);
end
% Calculation of Indiviual Diffusion (m^2/sec)
Diff(j,k)=(0.000000186*Tin^1.5*EffMol(j,k))/(Pres*(9.87*10^-6)*EffDia(j,k)^2*Ohm(j,k));
% Assign binary diffusions equal to zero when it goes to infinity
if Diff(j,k)==inf
Diff(j,k)=0;
end
Naa=isnan(Diff(j,k));
if Naa==1
Diff (j,k)=0;
end
end
end
% Calculating Molecular Diffusivity of the species in mixture
for j=1:5
for k=1:5
if j~=k
% Diffusion in Mixture
DiffMix(j)= DiffMix(j)+(X(1,k)/Diff(j,k));
end
end
DiffMix(j)=((1-Y(1,j))/DiffMix(j));
if DiffMix(j)==inf
DiffMix(j)=0;
end
Na=isnan(DiffMix(j));
if Na==1
DiffMix(j)=0;
end
end
% Molecular diffusivity of each species (m^2/sec)
for j=1:5
Dm(j) = DiffMix(j);
end
%%%%%%%%%%%%%%%%% CALCULATION OF KNUDSEN DIFFUSIVITY %%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculation of Knudsen Diffusivity of each species (m^2/sec)
% Ref: Hayes and Kolaczkowski_ Introduction to catalytic combustion
for j=1:5
DKn(j) = (dp/3)*sqrt((8* Ru * 10^3 * Tin)/(3.14* MolMass(j)));
if DKn(j)==inf
DKn(j)=0;
end
Na=isnan(DKn(j));
if Na==1
DKn(j)=0;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% To check the ratio of the diffusivities of each species less than 0.1 for
% having reasonably very less deviation on the radial direction
% Ref: Desmet et al.2003_Chem.Eng.Sci_V58_3203
e(j)=0;
for j=1:5
e(j) = DKn(j) /Dm(j);
end
if e(j) < 0.1
disp ('The ratio of diffusivity is less than 0.1 and lumped model of the species equation can be used');
end
The ratio of diffusivity is less than 0.1 and lumped model of the species equation can be used
%%%%%%%%%%%%%%%% CALCULATION OF EFFECTIVE DIFFUSIVITY %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Effective diffusivities of the species in m^2/sec (series Pore Model)
Deff(j)=0;
for j=1:5
Deff(j) = ((DKn(j)*Dm(j))*E)/((DKn(j)*TorF)+(KnTorF*Dm(j)));
end
% Schmidt number for each species
Sc(j)=0;
for j=1:5
Sc(j) = (DyVis)/(DenF*Deff(j));
if Sc(j)==inf
Sc(j)=0;
end
Na=isnan(Sc(j));
if Na==1
Sc(j)=0;
end
end
%%%%%%%%%%%%%%%% CALCULATION OF REYNOLDS NUMBER %%%%%%%%%%%%%%%%%%%%%%%%
% Reynolds number for packed bed reactor
% When Reynolds number of the reactor is greater than 50, we can neglect
% axial dispersion of heat in the energy equation. Ref: Borkink,J.G.H_1991
Re = (D*u*DenF)/(DyVis);
if Re > 20
disp('Neglect Axial dispersion effect in the energy equation');
end
Neglect Axial dispersion effect in the energy equation
% Reynolds number for packed bed particles
Rep = (db*u*DenF)/(DyVis);
%%%%%%%%%%%%%%%% CALCULATION OF MASS BIOT NUMBER %%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculation of j-factor for finding species mass transfer
jD = (0.667)/(E*((Rep)^0.481));
% Calculation of Species mass transfer for each species (m/s)
km = zeros(5);
for j=1:5
km(j) = (jD*E*u)/((Sc(j))^0.66667);
if km(j)==inf
km(j)=0;
end
Na=isnan(km(j));
if Na==1
km(j)=0;
end
end
% Characteristic length (m)
Lm = (Tor*(db/2));
% To ensure the model validity - Checking Mass Transfer Biot number for
% each species
Bim = zeros(5);
for j=1:5
Bim(j) = ((km(j) * Lm)/Deff(j));
Na=isnan(Bim(j));
if Na==1
Bim(j)=0;
end
end
% Mass Biot number for each species should be less than 10 for neglecting
% radial diffusion effect-Ref: Venderbosch et al.Chem.Eng.Sci._V53_19_3355.
if Bim(j) < 10
disp ('Mass Transfer Biot number is less than 10');
disp ('Neglect diffusion effect of the species in the radial direction');
else
disp ('1-D model is invalid')
disp ('Please use 2-D model with radial effect')
pause(30)
end
Mass Transfer Biot number is less than 10
Neglect diffusion effect of the species in the radial direction
%%%%%%%%%%%%%%%% CALCULATION OF HEAT BIOT NUMBER %%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Characteristic length for finding thermal biot number in (m)
Lt = D/4;
% Nusselt number for finding heat transfer coefficient using Li and
% Finlayson Correlation (1977)
Nu = 0.17*(Rep)^0.79;
% Heat transfer coefficient on outside of the reactor (W/m2/K)
hcoff = (Nu*Kg)/Lt;
% Parallel flow model to find effective conductivity (W/m/K)
% Ref: Nield and Bejan_2006_Convection in Porous Media
Keff = Kg*(1+(3*(1-E)*((Kp/Kg)-1))/(3+E*((Kp/Kg)-1)));
% Keff = (E*Kg)+((1-E)*Kp);
% To neglect the external heat transfer term check thermal Biot number
% Ref: Finlayson1971_Packedbedanalysis_Chem.Eng.Sci_V26_1081, Ref: Ferguson
% andFinlayson1974_NO-COmodeling_AIChE_V20_3_539
Bih = (hcoff*D)/(2*Keff);
if Bih < 1
disp ('Wall Heat transfer Biot Number is less than 1 and so, neglect the radial effect of temperature change in the reactor');
else
disp ('1-D model is invalid')
disp ('Please use 2-D model with radial effect')
pause(30)
end
Wall Heat transfer Biot Number is less than 1 and so, neglect the radial effect of temperature change in the reactor
% To neglect intraparticle gradients (or) lumping pellets and bulk into
% same temperature. The Particle Biot number condition will check this
% statement. Ref: Ruud.J.Wijngaarden and
% K.RoelWesterterp1993_Chem.Eng.Technol.V16_363
Bihp = (hcoff*db)/Kp;
if Bihp < 0.1
disp('Temperature in bulk and pellet surface can be lumped in the energy equation')
else
disp('Two phase temperatures for pellet and bulk should be used for the energy equation')
pause(30)
end
Temperature in bulk and pellet surface can be lumped in the energy equation
% If the Biot number ratios of mass and heat transfer is approximately
% equal to ten, then the concentration and temperature profiles over the
% catalyst radius have very less gradients. ref: Venderbosch et al.
% Chem.Eng.Sci._V53_19_3355.
Biratio(j)=0;
for j=1:5
Biratio(j) = Bim(j)/Bih;
Biratio(j)= round(Biratio(j));
end
if ((Biratio(1))&&(Biratio(2))) < 10
disp ('The lumped model derived for the packed bed model can be used');
end
The lumped model derived for the packed bed model can be used
%%%%%%%%% Calculating Constant Specific Heat for Gas at Node 1: %%%%%%%%%%%
Cp(1,1)=Ru*(a11+a12*Tin+a13*Tin*Tin+a14*Tin*Tin*Tin+a15*Tin*Tin*Tin*Tin);
Cp(1,2)=Ru*(a21+a22*Tin+a23*Tin*Tin+a24*Tin*Tin*Tin+a25*Tin*Tin*Tin*Tin);
Cp(1,3)=Ru*(a31+a32*Tin+a33*Tin*Tin+a34*Tin*Tin*Tin+a35*Tin*Tin*Tin*Tin);
Cp(1,4)=Ru*(a41+a42*Tin+a43*Tin*Tin+a44*Tin*Tin*Tin+a45*Tin*Tin*Tin*Tin);
Cp(1,5)=Ru*(a51+a52*Tin+a53*Tin*Tin+a54*Tin*Tin*Tin+a55*Tin*Tin*Tin*Tin);
Cpmix=0;
% Calculation of Constant Pressure Specific Heat Mol Basis (KJ/Kmol.K)
for j=1:5
Cpmix=Cpmix+Cp(1,j)*X(1,j);
end
% Calculation of Constant Pressure Specific Heat kg basis (J/Kg.K)
Cpmixmass=Cpmix*1000/MolMassMix;
% Volumetric ratio of heat capacities by Ferguson
% andFinlayson1974_NO-COmodeling_AIChE_V20_3_539
Cr = (E*DenF*Cpmixmass)/((1-E)*DenP*Cpp);
if Cr < 0.002
disp ('Fluid phase heat capacity term in the energy equation can be neglected');
else
disp ('Fluid phase heat capacity term in the energy equation is considered in the energy equation and not neglected');
end
Fluid phase heat capacity term in the energy equation can be neglected
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Simulation paused for 40 secs to view the limitation conditions used for the model');
Simulation paused for 40 secs to view the limitation conditions used for the model
%pause(40)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialization of various parameters of the reactor bed before simulation
for i=1:h
% Initial Temperature Values Across the reactor in K
Tia(i)= 303;
% Initialization of Reaction Rate at Various Nodes
RR(i)=0;
% Partial Pressure of the reactants (atm)
P(1)=0; % Inlet NO partial pressure
P(2)=0; % Inlet CO partial pressure
P(3)=0; % Inlet N2 partial pressure
P(4)=0; % Inlet CO2 partial pressure
P(5)=1; % Inlet NO partial pressure
% Initial Mole fractions of the reactants
X(i,1)=0; % Inlet NO Mole fraction
X(i,2)=0; % Inlet CO Mole fraction
X(i,3)=0; % Inlet N2 Mole fraction
X(i,4)=0; % Inlet CO2 Mole fraction
X(i,5)=1; % Inlet He Mole fraction
% Initializing inlet concentrations and mass fractions for all
% species at the initial condition
MolMassMix=0;
for j=1:5
% Molecular Mass of Individual Species in 1 Kmol of Mixture(Kg/Kmol)
MolMasss(i,j)= X(i,j)*Mol(j);
% Calculation of Molecular Mass of Mixture (Kg/Kmol)
MolMassMix=MolMassMix+MolMasss(i,j);
end
% Gas Constant of Mixture (J/Kg/K)
Rmass=(Ru*1000)/MolMassMix;
% Density of the fluid (Kg/m^3)
DenF=Pres/(Rmass*Tin);
% Inlet Concentration of Individual Species at first node (mol/m^3)
for j=1:5
Concc(i,j)=(DenF/MolMassMix)*X(i,j)*1000;
Na=isnan(Cin(1,j));
if Na==1
Concc(i,j)=0;
end
% Mass Fraction of Individual Species at first node
Y(i,j)=(X(i,j)*Mol(j))/MolMassMix;
Na=isnan(Y(i,j));
if Na==1
Y(i,j)=0;
end
end
end
% Option for ODE Solver
options=odeset('RelTol',1e-6,'Stats','on');
%%%%%%% Starting calculation for Initial time step at first node %%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initial Mole fractions of the reactants.
X(1,1)=0.005; % Inlet NO Mole fraction
X(1,2)=0.005; % Inlet CO Mole fraction
X(1,3)=0; % Inlet N2 Mole fraction
X(1,4)=0; % Inlet CO2 Mole fraction
X(1,5)=0.99; % Inlet He Mole fraction
% Initializing inlet concentrations and mass fractions for all
% species at the initial condition
MolMassMix=0;
Cinp = zeros(10,5);
Cinpp = zeros(10,5);
for j=1:5
% Molecular Mass of Individual Species in 1 Kmol of mixture (Kg/Kmol)
MolMasss(1,j)= X(1,j)*Mol(j);
% Calculation of Molecular Mass of Mixture (Kg/Kmol)
MolMassMix=MolMassMix+MolMasss(1,j);
end
% Temperature at node 1 (K)
Tia(1)=423;
% Gas Constant of Mixture (J/Kg/K)
Rmass=(Ru*1000)/MolMassMix;
% Density of the fluid (Kg/m^3)
DenF=Pres/(Rmass*Tia(1));
% Inlet Concentration of Individual Species at first node (mol/m^3)
for j=1:5
Concc(1,j)=(DenF/MolMassMix)*X(1,j)*1000;
Na=isnan(Cin(1,j));
if Na==1
Concc(1,j)=0;
end
end
for i=2:h
for j=1:5
Cinp(i,j)=Concc(i,j);
Cinpp(i,j)=Concc(i,j);
end
end
%%%%%%%% Calculating Constant Specific Heat for Gas at Node 1: %%%%%%%%%%%
Cp(1,1)=Ru*(a11+a12*Tia(1)+a13*Tia(1)*Tia(1)+a14*Tia(1)*Tia(1)*Tia(1)+a15*Tia(1)*Tia(1)*Tia(1)*Tia(1));
Cp(1,2)=Ru*(a21+a22*Tia(1)+a23*Tia(1)*Tia(1)+a24*Tia(1)*Tia(1)*Tia(1)+a25*Tia(1)*Tia(1)*Tia(1)*Tia(1));
Cp(1,3)=Ru*(a31+a32*Tia(1)+a33*Tia(1)*Tia(1)+a34*Tia(1)*Tia(1)*Tia(1)+a35*Tia(1)*Tia(1)*Tia(1)*Tia(1));
Cp(1,4)=Ru*(a41+a42*Tia(1)+a43*Tia(1)*Tia(1)+a44*Tia(1)*Tia(1)*Tia(1)+a45*Tia(1)*Tia(1)*Tia(1)*Tia(1));
Cp(1,5)=Ru*(a51+a52*Tia(1)+a53*Tia(1)*Tia(1)+a54*Tia(1)*Tia(1)*Tia(1)+a55*Tia(1)*Tia(1)*Tia(1)*Tia(1));
Cpmix=0;
% Calculation of Constant Pressure Specific Heat Mol Basis (KJ/Kmol.K)
for j=1:5
Cpmix=Cpmix+Cp(1,j)*X(1,j);
end
% Calculation of Constant Pressure Specific Heat kg basis (J/Kg.K)
Cpmixmass=Cpmix*1000/MolMassMix;
% 'AConstant' is used in the energy equation
Aconstant = ((1- E)* DenP * Cpp)+(E*DenF*Cpmixmass);
% 'Bconstant' for finding time step in energy equation
Bconstant = (E*DenF*Cpmixmass)/(((1- E)* DenP * Cpp)+(E*DenF*Cpmixmass));
%%%%%%%%%% Defining the Time Difference for Temperature Modeling %%%%%%%%%%
delt=(2*Keff)/(Aconstant*(Bconstant^2)*(u^2));
delt1 = ((dx^2)*(Aconstant))/(Keff);
% The above equation is used for finding time step which provides the
% minimum value for the time step and thus, used for modeling purposes
% Defining the Time Step as Per Input File
if delt1 > 0.25
TimeStep=0.25;
else
TimeStep=0.1;
TempStep = (2/60)* TimeStep;
disp ('Change the Temperature step value to:');
disp(TempStep);
pause(30)
end
if (TimeStep < delt)&&(TimeStep < delt1)
dt=TimeStep;
else
disp ('Check the time step for stability conditions');
pause (100)
end
% Defining the Initial Value for Time
a=1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% TEMPERATURE AND SPECIES MODELING FROM NODES '2' to 'h' FOR %%%%%
%%%%% EACH TIME STEP - INCREMENT BASED ON STABILITY CONDITIONS %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Pre-allocate variables for the maximum possible time and space steps
Hp = zeros(10,5);
THR = zeros(10,1);
Hcond = zeros(10,1);
Hconv = zeros(10,1);
Tf = zeros(10,1);
Cinn = zeros(10,5);
B = zeros(5);
C = zeros(10,5);
MolMasss = zeros(10,5);
Dm = zeros(10,5);
DKn = zeros(10,5);
for T=423:0.0083333:573 % 2 degree K/min
% If the time step(dt) is 0.1, then use temperature rise step as 0.0033
% Initial Mole fractions is equal to the initial partial pressures of
% NO and CO provided and the reactor operates at 1 atm pressure
% Calculation of Pellet Temperature at Next Time Step.
for i=2:h
% Calculation of Enthalpy for each species (J/kg)
% Calculation of Enthalpy for NO
Hp(i,1)=Ru*Tia(i)*(a11+((a12*Tia(i))/2)+(a13*Tia(i)*Tia(i)/3)+(a14*Tia(i)*Tia(i)*Tia(i)/4)+(a15*Tia(i)*Tia(i)*Tia(i)*Tia(i)/5)+(a16/Tia(i)));
% Calculation of Enthalpy for CO
Hp(i,2)=Ru*Tia(i)*(a21+((a22*Tia(i))/2)+(a23*Tia(i)*Tia(i)/3)+(a24*Tia(i)*Tia(i)*Tia(i)/4)+(a25*Tia(i)*Tia(i)*Tia(i)*Tia(i)/5)+(a26/Tia(i)));
% Calculation of Enthalpy for N2
Hp(i,3)=Ru*Tia(i)*(a31+((a32*Tia(i))/2)+(a33*Tia(i)*Tia(i)/3)+(a34*Tia(i)*Tia(i)*Tia(i)/4)+(a35*Tia(i)*Tia(i)*Tia(i)*Tia(i)/5)+(a36/Tia(i)));
% Calculation of Enthalpy for CO2
Hp(i,4)=Ru*Tia(i)*(a41+((a42*Tia(i))/2)+(a43*Tia(i)*Tia(i)/3)+(a44*Tia(i)*Tia(i)*Tia(i)/4)+(a45*Tia(i)*Tia(i)*Tia(i)*Tia(i)/5)+(a46/Tia(i)));
% Calculation of Enthalpy for He
Hp(i,5)=Ru*Tia(i)*(a51+((a52*Tia(i))/2)+(a53*Tia(i)*Tia(i)/3)+(a54*Tia(i)*Tia(i)*Tia(i)/4)+(a55*Tia(i)*Tia(i)*Tia(i)*Tia(i)/5)+(a56/Tia(i)));
% Total Heat Released from NO - CO reaction
THR(i)=(-Hp(i,1)*RR(i))-(Hp(i,2)*RR(i))+(Hp(i,3)*RR(i)*0.5)+(Hp(i,4)*RR(i));
if i<h
% Heat From Conduction
Hcond(i)=(((Keff*dt)/(Aconstant*dx*dx))*(Tia(i+1)-(2*Tia(i))+(Tia(i-1))));
% Heat From Convection
Hconv(i)=(((u*Bconstant*dt)/(2*dx))*(Tia(i+1)-Tia(i-1)));
% Calculations of pellet reactor at new node
Tf(i)=Tia(i)-Hconv(i)+ Hcond(i)-((dt/Aconstant)*THR(i));
else
% Heat From Conduction
Hcond(i)=(((2*Keff*dt)/(Aconstant*dx*dx))*(Tia(i-1)-Tia(i)));
% Heat From Convection
Hconv(i)=0;
% End Boundary Condition for pellet reactor at the end
Tf(i)=Tia(i)-Hconv(i)+ Hcond(i)-((dt/Aconstant)*THR(i));
end
end
Tia(1)=T;
for i=2:h
Tia(i)=Tf(i);
end
%%%%%%%%%%%%%%%%%%%%%%% SPECIES MODELING %%%%%%%%%%%%%%%%%%%%%%%%%
% While loop to make the species concentrations from nodes '2' to 'h'
% steady state for the corresponding temperature
delc = 1;
while delc > 0.0000000001
for i=2:h
for j=1:5
% Molecular Mass of Individual Species in 1 Kmol of Mixture (Kg/Kmol)
MolMasss(i,j)= X(i,j)*Mol(j);
% Calculation of Molecular Mass of Mixture (Kg/Kmol)
MolMassMix=MolMassMix+MolMasss(i,j);
end
% Calling Molecular diffusivity calculator Function
[DiffMix]=MolecularDiffusion(Tia,X,Y);
for j=1:5
% Molecular diffusivity of each species at each node (m^2/sec)
Dm(i,j)=DiffMix(j);
% Knudsen Diffussivity of each species at each node (m^2/sec)
DKn(i,j) = (dp/3)*sqrt((8* Ru * 10^3 * Tia(i))/(3.14* MolMasss(i,j)));
if DKn(i,j)==inf
DKn(i,j)=0;
end
Na=isnan(DKn(i,j));
if Na==1
DKn(i,j)=0;
end
% Effective diffusivity of each species at each node (m^2/sec)
Deff(i,j) = ((DKn(i,j)*Dm(i,j))*E)/(DKn(i,j)*TorF)+(KnTorF*Dm(i,j));
end
% Solving Concentration at each nodes for the corresponding temperature
delb = 1;
% While loop used to obtain steady state solution
while delb > 0.0000000001
% Calling ODE Solver for the corresponding node
[t,Z]=ode15s(@SpeciesODE, [0 dt],Cinp(i,:), options);
[m,n]=size(Z);
% Putting the Values back in designated species
Cinn(i,1)=Z(m,1);
Cinn(i,2)=Z(m,2);
Cinn(i,3)=Z(m,3);
Cinn(i,4)=Z(m,4);
Cinn(i,5)=Z(m,5);
for j=1:5
% Calculate error to make concentration steady state at that node
B(j)=abs(Cinn(i,j)-Cinp(i,j));
Concc(i,j)=Cinn(i,j);
Cinp(i,j)=Cinn(i,j);
end
% The maximum value in the array is compared with 'delb'
delb=max(B);
end
end
for i=2:h
for j=1:5
% Calculate error to make concentration steady state at all nodes
C(i,j)=abs(Concc(i,j)-Cinpp(i,j));
Cinpp(i,j)=Concc(i,j);
end
end
% The maximum value in the array is compared with 'delc'
delc=max(C);
end
% Store Temperature each cycle
Tempp(a,1)=T;
% Store: NO and CO conversions of the packed bed reactor (%)
NOConv(a,1)=((Concc(1,1)-Cinpp(h,1))/Concc(1,1))*100;
COConv(a,2)=((Concc(1,2)-Cinpp(h,2))/Concc(1,2))*100;
% Increment the counter step by 1
a = a + 1;
end
Error using odearguments
SPECIESODE must return a column vector.

Error in ode15s (line 153)
odearguments(odeIsFuncHandle, odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
% Plot for Temperature Vs Concentration of NO
figure;
plot(Tempp(:,1),NOConv(:,1))
xlabel('Inlet Temperature (K)')
ylabel('Conversion of NO in (%)')
legend('NO Conversion Curve')
axis([423 573 0 100])
title('Conversion Curve for NO reduction by CO')
% Plot for Temperature Vs Concenration of CO
figure;
plot(Tempp(:,1),COConv(:,2))
xlabel('Inlet Temperature (K)')
ylabel('Conversion of CO in (%)')
legend('CO Conversion Curve')
axis([423 573 0 100])
title('Conversion Curve for NO reduction by CO')
% Provides total process time of the model
SimulationTime = toc;
%End Program
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%% ODE SOLVER FOR SPECIES EQUATION %%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ODE Function to solve concentration of the species for the corresponding
% node and temperature
function dZdt = SpeciesODE(t,Z)
% Defining Global Variables
global E
global Deff
global Tia
global i
global u
global Concc
global Ru
global dx
global RR
global h
global Pres
% Assigning Numbers to Species
% 1:NO
% 2:CO
% 3:N2
% 4:CO2
% 5:He
% Getting the Initial Values from array
Cs(1)=Z(1);
Cs(2)=Z(2);
Cs(3)=Z(3);
Cs(4)=Z(4);
Cs(5)=Z(5);
% Values from Granger et. al.
% Pre-Exponential factor (mol/m^2.sec) => Conversion from (mol/gm/hr)
% ((7.4453657*10^21)*0.023)/(3600*3.14*0.012*0.003)
PreExp=0.33*10^21;
% Activation Energy in J/mol which is 47kcal/mol for NO dissociation step
Ea =196648;
% Per Atmosphere
ANO=2.4*10^-2;
% The Enthalpy is calculated using the equation given by Granger et al.
% 2001_TPDstudies_Topicsofcatalysis_V16-17_394.
% Enthalpy for NO and CO are calculated based on their partial pressures
% and correlation factors given in Granger et al.2001. This correlation
% takes care of the influence of partial pressures and temperature on the
% reactor
% Enthalpy of NO in J/mol
HNO=-48576;
% Equilibrium Constant of NO
KNO=ANO*exp(((-1)*(HNO))/(Ru*Tia(i)));
% Per Atmosphere
ACO=1.4*10^-2;
% Enthalpy of CO in J/mol
HCO=-32554;
% Equilibrium Constant of CO
KCO=ACO*exp(((-1)*(HCO))/(Ru*Tia(i)));
%Partial pressures for calculating reaction rates
Ctot = 0;
for j=1:5
Ctot = Ctot + Cs(j);
end
PP(1)=(9.86923*(10^-6)*Pres)*(Cs(1)/Ctot);
PP(2)=(9.86923*(10^-6)*Pres)*(Cs(2)/Ctot);
if PP(1)<0
PP(1)=0;
end
if PP(2)<0
PP(2)=0;
end
% Reaction Rate Expression for NO dissociation step
% Ref: AnandSrinivasanandDr.Depcik2010_Cat.Rev.Sci.Eng._V52_1-32
RR(i)=((PreExp * exp(-1*((Ea)/(Ru*Tia(i))))*(KNO*PP(1)))/((1+(KNO*PP(1))+(KCO*PP(2)))^2));
%Calculation of Concentrations for the corresponding time step
if i<h
dZdt(1)=-((((Cs(1))-(Concc((i-1),1)))*E*u)/(dx))+(Deff(i,1)*E*((Concc((i+1),1))-(2*Cs(1))+(Concc((i-1),1)))/(dx^2))-(RR(i));
dZdt(2)=-((((Cs(2))-(Concc((i-1),2)))*E*u)/(dx))+(Deff(i,2)*E*((Concc((i+1),2))-(2*Cs(2))+(Concc((i-1),2)))/(dx^2))-(RR(i));
dZdt(3)=-((((Cs(3))-(Concc((i-1),3)))*E*u)/(dx))+(Deff(i,3)*E*((Concc((i+1),3))-(2*Cs(3))+(Concc((i-1),3)))/(dx^2))+(0.5*RR(i));
dZdt(4)=-((((Cs(4))-(Concc((i-1),4)))*E*u)/(dx))+(Deff(i,4)*E*((Concc((i+1),4))-(2*Cs(4))+(Concc((i-1),4)))/(dx^2))+(RR(i));
dZdt(5)=-((((Cs(5))-(Concc((i-1),5)))*E*u)/(dx))+(Deff(i,5)*E*((Concc((i+1),5))-(2*Cs(5))+(Concc((i-1),5)))/(dx^2));
% Concentration calculation for the node i=h (or) last node of
% the reactor using Neumann boundary condition
else
dZdt(1)=-((((Cs(1))-(Concc((i-1),1)))*E*u)/(dx))+(Deff(i,1)*E*((2*Concc((i-1),1))-(2*Cs(1)))/(dx^2))-(RR(i));
dZdt(2)=-((((Cs(2))-(Concc((i-1),2)))*E*u)/(dx))+(Deff(i,2)*E*((2*Concc((i-1),2))-(2*Cs(2)))/(dx^2))-(RR(i));
dZdt(3)=-((((Cs(3))-(Concc((i-1),3)))*E*u)/(dx))+(Deff(i,3)*E*((2*Concc((i-1),3))-(2*Cs(3)))/(dx^2))+(0.5*RR(i));
dZdt(4)=-((((Cs(4))-(Concc((i-1),4)))*E*u)/(dx))+(Deff(i,4)*E*((2*Concc((i-1),4))-(2*Cs(4)))/(dx^2))+(RR(i));
dZdt(5)=-((((Cs(5))-(Concc((i-1),5)))*E*u)/(dx))+(Deff(i,5)*E*((2*Concc((i-1),5))-(2*Cs(5)))/(dx^2));
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%% MOLECULAR DIFFUSIVITY %%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function for calculating Molecular Diffusion for NO-CO species
function [DiffMix]=MolecularDiffusion(Tia,X,Y)
global i
global EffDia
global EffMol
global EffKBEp
global Pres
KBTEp=zeros(5,5);
Ohm=zeros(5,5);
Diff=zeros(5,5);
DiffMix=zeros(5);
% Binary Diffusion of Species (m^2/sec)
for n=1:5
for k=1:5
% Inverse Energies
KBTEp(n,k)=(EffKBEp(n,k)*Tia(i));
% Calculation of Omega
if KBTEp(n,k)<5
% Equation from Curve Fitting in Excel
Ohm(n,k)=1.4803*(KBTEp(n,k)^-0.397);
else
% Equation from Curve Fitting in Excel
Ohm(n,k)=1.0765*(KBTEp(n,k)^-0.16);
end
% Calculation of Binary Diffusion (m^2/sec)
Diff(n,k)=(0.000000186*Tia(i)^1.5*EffMol(n,k))/(Pres*(9.87*10^-6)*EffDia(n,k)^2*Ohm(n,k));
% 'Not-a-number' values are checked in the array and assigned to be
% zero
if Diff(n,k)==inf
Diff(n,k)=0;
end
Naa=isnan(Diff(n,k));
if Naa==1
Diff (n,k)=0;
end
end
end
% Diffusion of Species in mixture
for n=1:5
% Initialization
DiffMix(n)=0;
for k=1:5
if n~=k
% Diffusion in Mixture
DiffMix(n)= DiffMix(n)+(X(i,k)/Diff(n,k));
end
end
DiffMix(n)=((1-Y(i,n))/DiffMix(n));
% 'Not-a-number' values are checked in the array and assigned to be zero
if DiffMix(n)==inf
DiffMix(n)=0;
end
Na=isnan(DiffMix(n));
if Na==1
DiffMix(n)=0;
end
end
end

Accepted Answer

Torsten
Torsten on 28 Mar 2023
Add
dZdt = dZdt.';
as last line in function "SpeciesODE".
  2 Comments
Anand
Anand on 28 Mar 2023
thank you !!
however i have encountered another problem,
Warning: Failure at t=0.000000e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (7.905050e-323) at time t.
> In ode15s (line 661)
In anand_shrivastan (line 625)
Torsten
Torsten on 29 Mar 2023
Edited: Torsten on 29 Mar 2023
This means that technically, your code is correct, but the integrator has problems integrating your equations because of many possible reasons. The first thing I would do is to check whether the differential equations you want to solve are set up correctly in your code.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!