12 replies to this topic

[Huba]

Hi,

 

I need to simulate a catalytic packed bed reactor. Its for the oxidation of SO2, so according to litterature, it's supposed to be adiabatic. I'm using the ODE solver on Matlab for this problem, but how do I input the energy balance , if it's not even an ODE?

 

According to the book I've used in my reaction course (Fogler, S. Elements of Chemical Reaction Engineering), the energy balance is

 

 

T = T0 + (Hrx*X)/(Σ θi* Cpi)

 

I apologize if the question seems more like a lack of practice on Matlab, but I've never had to solve an adiabatic reactor this way before.

[thorium90]

You put your balances wrt to the reactor length, z.

Attached Files

•  equations.jpg   48.4KB   23 downloads

Edited by thorium90, 24 November 2013 - 12:45 AM.

[Huba]

Ok, and eta is effectiveness factor?

[thorium90]

Yes.

[Huba]

I have tried with the above equations, but frankly the results are strange, the conversion is around 25% for an increase in temperature of less than 10 degrees. Theoretically, conversions go up to 70% with changes of temperature from 1260 to 1500 degrees R.

 

Can someone review my matlab code?

 

The kinetic constants are taken from an article where they also do a simulation on Matlab, but the reactor is multitubular (non adiabatic). If I manually test these equations with certain temperature, the resulting rate is pretty low

 

Also, I'm a bit confused about the Ergun equation. if epsilon is the catalyst porosity, what is phi?

 

I apologize if most of the comments are in french, I'd really appreciate any kind of help ! I've been working on this for a while now...

Attached Files

•  Matlab_SO2.zip   2.54KB   94 downloads

Edited by QDanh, 27 November 2013 - 10:42 PM.

[breizh]

Info about Ergun equation .

 

Breizh

[curious_cat]

 

Also, I'm a bit confused about the Ergun equation. if epsilon is the catalyst porosity, what is phi?

 

Spheracity, I think. 

[thorium90]

The cat is correct. I thought the form of the equations was pretty obvious for a student of chemical engineering, nevertheless;

 

phi is sphericity.

epsilon is bed porosity, not catalyst porosity.

eta is effectiveness factor

http://www.google.co...9,d.bmk&cad=rja

 

I'm finding it hard to understand your code, firstly because its not in english although the units give some clue as to what it is and so I just have to assume the numbers and units are correct, secondly, its not mentioned what the numbers are.

 

For gas phase reactions, superficial velocity, U, changes wrt density which changes as a function of pressure and temperature and the amount of the different components passing that dz.

 

You should check the formulation of the ODEs, it is incorrect. Also, it should have the units of per unit length, Temperature/Length and Pressure/Length.

Edited by thorium90, 28 November 2013 - 01:42 AM.

[curious_cat]

Adding to @thoriums last comment, it would help if you elaborated what your solution strategy is. It's somewhat hard to try and interpret your code. 

 

Write down your equations and post. How many ODE's do you have? How many variables? what are initial conditions etc. 

 

Is your material balance tallying? These are failrly standard checks to test the goodness of your code. Bugs are a bigger problem than accuracy of kinetic data. 

[Huba]

Basically, I have 3 ODEs (dX/dz, dT/dz, dP/dz), z for bed length (or height). So I have 3 variables: X (conversion), T(temperature), P(pressure).

 

The initial conditions are : X=0 % of oxidized SO2, feed temperature of around 430-500 C(converted in R), initial pressure of 1.2 atm.

 

The feed flow is 56113.46 m^3/h, which i converted to m^3/s.

The initial composition is (vol. %) : 12% SO2, 9% O2, 79% N2

 

With this, I can calculate molar flows.

 

The bed volume is 11.8 m^3, and bulk density is 33.8 lb/ft^3.

So it's 11.8 * 35.3146667 to convert in ft^3

 

 

 

I think the strategy here is mainly, for each step:

 

1) Set initial conditions and parameters

 

2) Calculate temperature dependant variables (k, Kp, Cp, dHRx)

 

 

3) Calculate reaction rate:

 

- for below 5 % of oxidized SO2, the rate equation is : -r = k*(0.848 - 0.012/Kp^2)            (1)

 

 

- for above 5%:    r = -k * (P_SO2/P_SO3)^0.5 * [P_O2 - (P_SO3/(Kp*P_SO2))^2]             (2)

 

where P_i = P_SO20 * (theta_i*+mu*X) / (1+e*X) * (P/P0) * (T0/T)

 

Here i assumed that e = 0 (since dPdz does not depend on conversion). So by substituting this last equation in equation (2), I should get

 

 

r = -k * (1-X/X)^0.5 * [ (P/P0) * P_SO20 * (thetaO2-X/2) * (T0/T) - (X/(1-X))^2 * (1/Kp^2) ] 

 

(although I took out the (T0/T) for testing, it didn't make much of a difference).

** r is expressed in lbmol SO2/(lb. cat*s).

 

 

 

4) Solve balance equations

 

Mass : dX/dz = (-r)/F_SO20 * rho_k * (surface), which is in (-)/ft,       

 

Since r is per unit mass, I need to multiply by density then by surface to have unit length. (1-epsilon) is to substract the void fraction.

 

 

Energy: dT/dz = [ ((-r)*dHRx) * rho_k * (surface) ] / sum (F_i*Cp_i),   in degrees R/ft

 

 

Momentum : basically Ergun equation.

[curious_cat]

Your strategy seems ok. 

 

So at <5% stage rate is independent of all concentrations? Your Eq. 1 doesn't seem to have a conc. term. Just checking. 

 

At 5% do (1) & (2) give you identical rates? 

[curious_cat]

You can try simulating a isothermal case ignoring Pressure drop & simulate that first. If it works sensibly at least you've validated the kinetic equations. 

 

Then add the complexity of the adiabatic case & finally the ergun equation. 

 

Always easier to debug in smaller bites. 

[thorium90]

Im still guessing what is (0.848 - 0.012) in your rate equation.
For the second part of your rate equation, it uses partial pressures. You need to calculate partial pressures by calculating from the conversion for that dz, not directly put X. X is conversion not partial pressure.
 
Im also puzzling on
PA0=(FA0/FT0)*1.2;               % en atm
 
I take it your reaction is S + O2 => SO2  ?
 
Reactants have 2 moles and products have one mole so clearly you cant use a fixed molar flow rate. You have to calculate the molar flow for each dz.
 
For the same reason as above, density and superficial velocity will be different which results in pressure drop and therefore your ergun equation calculates pressure that way.

As i described, it would therefore actually be more complicated to simulate an isothermal and or isobaric reactor, since the reaction is exothermic and the pressure changes, one would therefore have to add additional terms to remove heat for isothermal and maintain pressure (somehow) with additional terms too.
 
Ta=1265.67;                                                      % en R (298 K)
dHRx0 = -99*0.9478/453.59;                                       % en BTU/lbmol
 
Is 298K = 1265.67R ?
How is dHRx0 = -99*0.9478/453.59; equals heat of formation?
Either take it from a table or calculate it using the formula for heat of formation.
http://chemwiki.ucda...py_Of_Formation
 
I assume this is supposed to be Hess Law?
dHRx = dHRx0 + (Cp(4)-Cp(1)-Cp(2))*(T-Ta);                       % en BTU/lbmol (tient pas compte de N2**)
 
That is incorrect. It should have been;
Heat of reaction(T) = (Std heat of formation) + (Enthalpy change(T) - Enthalpy change(298K))
http://en.wikipedia....wiki/Hess's_law
 
The ODEs, you have put epsilon all over the place in the wrong places.
 
Also, the things mentioned in my post 8.
 
You can repost here when you have made the necessary corrections to your code.

Edited by thorium90, 29 November 2013 - 05:10 AM.