clear all;
tic %start clock time
% 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');
disp('Length of the catalyst bed in 0.003 m');
h=input('Enter the Discretization number:');
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
%%%%%%%%%%%%%%%% 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
% 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
%%%%%%%%%%%%%%%% 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
% 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
% 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
%%%%%%%%% 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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('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
% 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
dZdt=dZdt';
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
