Jump to content



Featured Articles

Check out the latest featured articles.

File Library

Check out the latest downloads available in the File Library.

New Article

Product Viscosity vs. Shear

Featured File

Vertical Tank Selection

New Blog Entry

Low Flow in Pipes- posted in Ankur's blog

Plug Flow Reactor Design With Matlab


This topic has been archived. This means that you cannot reply to this topic.
12 replies to this topic
Share this topic:
| More

#1 owx1985

owx1985

    Brand New Member

  • Members
  • 7 posts

Posted 27 February 2013 - 01:48 AM

I am designing a PFR reactor using MatLab, need some advise.
 
Below are my reaction
 
A--> B + C
2B --> D (undesired Products)
 
Kinetics are as below
1) rA = k1CA  k1 = 1e13*exp(-40000/RT)
2) rB = k2CB k2 = 1e26exp(-80000/RT)
 
I'm going to use performance eqn as below

FYI: my PFR will be integrated with the furnace, hence, my simulation will divide into 3 parts, the stack zone where the PFR will be preheated in the stack zone of the furnace follow by Convection zone and lastly the radiation zone.

 
dVdCA = - FaoCao/-ra
 
Fao = 30 kmol/hr; and Cao = 3379 mol/m3
This variable i extracted from Hysys during simulation
 
Now, i want to simulate the performance equation using matlab ode45 to acheive my CA and volume of the reactor. But i have a problem, how do i interlink both equation into my performance eqn?
 
 

Edited by owx1985, 27 February 2013 - 01:54 AM.


#2 thorium90

thorium90

    Gold Member

  • Members
  • 1,073 posts

Posted 27 February 2013 - 05:00 AM

Based on your need for a furnace, I assume both are endothermic reactions. Do they require the high temperatures in a furnace for the reaction to work? I assume the reactions are gas phase and catalytic. You might find it easier to assume a constant heat flux along the entire tube length in the radiant section. You will need your mass, energy and momentum balance equations for the matlab simulation not the performance equation for a PFR. Ergun equation is typically used for momentum balance, dP/dz (P for pressure, z for length). For mass balance, perhaps you can try based on conversion of each reaction, dX/dz (X for conversion of reaction j). Energy balance to solve by dT/dz. Take note that in the picture below, it does not have the heat flux term in it, but you can just add it in, just rmb to balance the units. You will also need Hess's law and the heats of formation. Components are i. Effectiveness factor is eta.

Attached Files


Edited by thorium90, 27 February 2013 - 05:02 AM.


#3 owx1985

owx1985

    Brand New Member

  • Members
  • 7 posts

Posted 27 February 2013 - 10:54 AM

My reaction will be in liquid phase, hence my reactor will be operating in high pressure condition, 23 bar. Maximum reaction temperature will be 340 C. Thermal cracking will start at around 250 C. If operate beyond the temperature, high amount of undesired products. My reaction kinetics is just first order reaction for first and 2nd reaction, hence it easy and can be simulate in the hysis. I manage to get the results. But for my case I want to optimize the whole reactor, in term best operating temp, diameter and length which mean volume with minimum energy input. My reaction will only starts at radiant zone, which means I assume my radiant zone is isothermal. My convection and shield zone will be a pre-heater for my reactor feed. I understand that material and energy balAnce are require for my matlab simulation, how to interlink this 2 eqn? What about the pressure eqn? Why do I need this eqn?

#4 owx1985

owx1985

    Brand New Member

  • Members
  • 7 posts

Posted 27 February 2013 - 11:00 AM

Can I have clearer picture on the eqns? Meaning the notation and how do i interlink the 2 reactions together or I just simply use below formulate

-rb = -ra+ rb?

#5 thorium90

thorium90

    Gold Member

  • Members
  • 1,073 posts

Posted 27 February 2013 - 11:06 AM

So the max rxn temperature is 340C but cracking starts at 250C and its liquid phase in the reactor? Seems odd... Do clarify what the inlet temperature is and it stays liquid in the entire reactor? Furnaces are normally used when high heat fluxes are required and if cracking already starts at 250C, you might be better off with an adiabatic reactor whose feed is heated in a shell and tube exchanger before entering. Since the reactions are endothermic, there is no "optimum" temperature, the higher the temperature the faster the kinetics and the smaller the reactor. The temperature will only be limited by things such as cracking, mechanical etc. You can get a bigger picture of the picture in post 2 by clicking on it...

 

The equations are calculated in every step, dz, the results of which are passed to the next step and so on in ode45. ode45 is like many IVP problems stringed together. You give the initial conditions, it calculates the results, this results goes to the initial guess for the next dz and so on until the reactor exit. You need the pressure equation as I assumed it is a catalytic reaction and there will be pressure drop through the reactor length

 

U is superficial velocity, phi is sphericity, dp is particle diameter, Cp is heat capacity, T_R is reference temperature, A is cross sectional area

 

what do you mean by,

 



-rb = -ra+ rb?


Edited by thorium90, 27 February 2013 - 11:21 AM.


#6 owx1985

owx1985

    Brand New Member

  • Members
  • 7 posts

Posted 27 February 2013 - 10:51 PM

Ok for first reaction:
ACO --> DO + AA (thermal cracking)
2nd reaction
2DO --> Gum
Production of DO = -r1- (-r2)
ACO = castor oil, DO = drying oil, Acetic acid.
Both reaction kinetics that I mentioned in post 1 applicable to liquid, hence my whole reaction need to operate at high pressure in order to make sure it is at liquid phase. Whereas for the pressure drop. My reaction no catalyst is required. Hence only high temperature reaction. Just like steam cracking. So meaning I just need to use material and energy balance for my reaction?

#7 thorium90

thorium90

    Gold Member

  • Members
  • 1,073 posts

Posted 27 February 2013 - 10:58 PM

A quick search off google reveals a number of relevant articles.

http://www.che.cemr....ng_oil/do-c.pdf

 

So you are designing both the furnace and the reactor? Based on the articles, if you design just the reactor, it would be a pretty straightforward task since you already have some kinetics on hand. The reactor is adiabatic and has inert packing to help radial mixing, which creates the pressure drop. This is a pressure driven process. Whenever a fluid flows from one place to another, there will be pressure drop.

Attached Files


Edited by thorium90, 27 February 2013 - 11:08 PM.


#8 owx1985

owx1985

    Brand New Member

  • Members
  • 7 posts

Posted 28 February 2013 - 10:20 AM

Hmm.. Thanks but to simulate the adiabatic reactor, is a bit too un realistic... Because the thermal cracking require high amount of heat. When the reaction starts... Temperature drop abruptly. Hence, prefer reactor will be isothermal. Hence, I simulate my reaction inside radiant zone of a furnace.

#9 thorium90

thorium90

    Gold Member

  • Members
  • 1,073 posts

Posted 28 February 2013 - 10:29 AM

I didnt calculate your heat of reaction, but you have done it in hysys already right? So you know if done in an adiabatic reactor what the temperature drop will be? Is it significant enough? The links suggest 40% per pass conversion, indicating large amounts of recycle are needed, its possible the temperature will not drop too much. You could try out both situations with some base conditions for comparison in hysys. The links also suggest not to go above 380C. You can do the simulation together in the furnace, you will need to modify the heat equation in post 2 slightly. In reality, it would appear an adiabatic reactor doesnt seem to be of much use if it can all be done in a furnace. I just felt the calculations will be easier to implement if done separately.

 

If you are ready to start, you can write out the full form of the equations in post 2 and post here, we can check if it is correct. (You need to expand it for components i, reactions j, references R, individual heat capacities and heats of formations).


Edited by thorium90, 28 February 2013 - 10:42 AM.


#10 owx1985

owx1985

    Brand New Member

  • Members
  • 7 posts

Posted 02 March 2013 - 11:25 PM

I tried to use the material balance in Matlab using

dX/dz = A*rj/Fo i convert into changes of in concentration. meaning i will get dCa/dz = pi*D^2*rA*Cao/Fo.

My initial conc for ACO which mean Cao = Fa/Vo = 3.47 kmol/m3. Using ODE45 in matlab my Conc decrease with length. The result doesnt seem to be the same as what i had in hysis.

As for the energy balance, should i use hess law or the 2nd eqn? if so what is deltaH? n? Cpi? Fi?

How to interlink the material balance and energybalnce?



#11 thorium90

thorium90

    Gold Member

  • Members
  • 1,073 posts

Posted 02 March 2013 - 11:57 PM

You have misunderstood the first equation. It is the rate of change of the conversion of reaction j at every step of z. Not the concentration. The concentration is calculated from the amount converted. Eg: if initially there is 1 mol of A and 100% A is converted in reaction 1, you will get 1mol B and 1 mol C.

 

Fi0 is the initial molar flow of i.

delta H_R(T) is the heat of reaction at temperature T. delta H_R(T_R) is the heat of reaction at reference conditions. delta H_f(T_R) is the heat of formation at reference conditions and these values are the ones you get from books like Perry's etc...

 

The n and m in the energy equation are the counters. Eg: from j=1 to j=n

 

At every step, you have to calculate the amount that is formed, ie: the example stated in the first paragraph of this reply. In addition, you have to calculate the new molar flows to the next dz, the new rate of reaction, the new concentration etc.

 

If you have done the matlab code already, perhaps you can paste it to notepad and attach here to review.


Edited by thorium90, 04 March 2013 - 07:22 AM.


#12 owx1985

owx1985

    Brand New Member

  • Members
  • 7 posts

Posted 26 March 2013 - 08:57 AM

For Heat of Reaction as below
Palmitic Acid (A)-- > 1-tetracedene [B] + Acetic Acid©
2 1-tetracedene --> Gum (D)
 
rA = -k1(CA); rB = k1(CA)-k2(CB^2); rC=k1(CA); rD=k2(CB^2)
 
Eg of one component
Heat of formation (B] = dH(B]  - dH(A); @ 25 oC
Heat of Reaction using Hess = dH1 + dH +dH2;
where dH1 = Int[Cp of A] from 613K to 298K.........................Cp = a + bT + cT^2 .....
dH = Heat of formation @ 298K or 25 oC.
dH2 = int[Cp of B] from 298K to 613K ...................................
 
I had 5 set of ODEs,
dCAdz= (A/Vo)*(rA);
dCBdz= (A/Vo)*(rB);
dCCdz=(A/Vo)*(rC);
dCDdz=(A/Vo)*(rD);
dTdz = -(A/Vo*Density*CpMass)*(-dHrxn1*rA + dHrxn2*rB +dHrxn3*rC+dHrxn4*rD+ [ 4U/D(Tsurr -T)]);
A = Area of the reactor pipe
Vo = Volumetric flow rate
CpMass = Specific Heat capacity of Mixture
U = Overall Heat Coefficient
D = diameter of reactor
T = Surrounding temperature
 
Command for ODE45
[X,Y]=ode45(@Rad2, [0:0.1:100], [3258;65.7;65.7;0;50]);
 
The matlab coding already send to your school mail

Attached Files


Edited by owx1985, 26 March 2013 - 08:59 AM.


#13 thorium90

thorium90

    Gold Member

  • Members
  • 1,073 posts

Posted 26 March 2013 - 01:55 PM

I think you should put the original matlab files here instead. Zip them up. Your code makes no sense. You put integration for t. t is not your temperature. dYdX(5) is your temperature. Y(1), Y(2) is not defined correctly either. It is obvious your code in its current state could not even run. The ode45 would be integrating a bunch of undefined constants... your code needs massive rewriting. Even passing of variables is not done correctly
I suggest you use conversion instead of concentration to reduce the number of ode to solve. Furthermore, your reactions involve phase change, assuming constant volume is incorrect. No attempt was made at calculating pressure changes too. No equation of state was used. No calculation for change in density and heat capacity as well as change in heat capacity was done either. The integration is not analytical, you should be using numerical integration like quad. 0.08m is not a standard diameter for tubes, although this is not too important, but 0.08m seems too large for such a small flow. Uniform radial profiles cannot possibly be a valid assumption in such a scenario. I think you have also messed up your mass and molar heat capacities too, you should check your units.
Theres quite a number of other problems too, but i suggest you follow through the above more critical problems first otherwise you will get very stuck...

Edited by thorium90, 26 March 2013 - 06:11 PM.





Similar Topics