Delta should go from -1 to 1, however the membrane reactor is showing better ..... Pi *Q. R*T. As a function of conversion, the partial pressures are defined at ...
WORCESTER POLYTECHNIC INSTITUTE
Membrane Reactor Modeling for Hydrogen Production through Methane Steam Reforming A Dissertation Submitted to the Faculty of the WORCESTER POLYTECHNIC INSTITUTE in partial fulfillment of the requirements for the Degree of Master of Science In Chemical Engineering By Jean-‐François ROUX 5/1/2011
Approved:
Professor N. Kazantzis, Advisor
Professor D. DiBiasio, Departement head
Acknowledgement Life is a funny thing. During these past two years, I had to be different person: a student, a teacher, a talker, a listener, a searcher, a finder… Every new experience brought me something unique, taught me new things, and made me who I am today. However, I could not have done and learn all of this alone, and it is with a deep pleasure I would like to thank all of them. My first thanks go to my advisor, professor Kazantzis, for his kindness and guidance. His support and belief in my abilities were invaluable and made me finish what I started when I arrived here. I would like to thank also my co-advisor, Professor Ma, for his councils and the financial support I received from him as a Research Assistant. And I like to thank the entire Chemical Engineering Department for its academicals and financial support, in particular all the teachers I worked for as a Teaching Assistant: professor DiBiasio and his French music taste, professor Deskins and his boxing fight, professor Clark and his quote(s), professor Dixon who let me grade in the “European way”, and professor Datta. It has been my pleasure to be teaching. I’m not forgetting Engin, whose work inspired my own research and who supported me for more than a year. Also, a real thank you to Ivan, Reyyan, Natalie, Alex, Mike and Susan with whom I shared group meetings and who gave me multiple help on my research. My thanks to Felicia and Tiffany for their smiles, and to Jack for his help in lab session. But work was not all. Indeed, I wish to thank everyone who support me during these two years, and survive all my speeches! Therefore, thank you Sophia and Will, Florent and Sarah, Olga and Oleski, Neal, Dan, Yuxian, Ahsan, Ivan, I. Oljora, Matt, Angel, Deepack, Angela, Sylvain, Federico, Evan, Chris, Ivan M., Marco, Kat, Jeff L., Karan, Vi, Alessandro and Giuseppe. My thanks also to Gaëlle, Mylène, Romain, Joanna, Aurore, Tipi, Ella, Vanessa, Évy, Alice and Mark for their thoughts and messages sent during these two years. And a very special thank you to Luca, the best Italian ever, Kathleen, my food partner, Thomas, my awesome roommate, and Jeff, my drinks partner. You are unique and my friends. I would finally like also to thank my mother for always believing in me whatever the situation, my brother for everything he did for me, my future-sister-in-law, and all my family for their continuous support. I would have not gone that far without all of you, I love you all. Thank you Dad, you will always be with me, in my heart. 1
Abstract A mathematical modeling framework for the methane steam reforming reaction operating in steady state has been developed. Performances are compared between the classic catalytic packed bed reactor and a Pd-based catalytic membrane reactor. Isothermal simulations on MATLAB © has first been conducted and showed a higher performance of the membrane reactor over the packed bed reactor. Methane conversion of 1 can be reached for lower temperatures than used with industrial PBR, and better performances are reached with an increase of the operating pressure. Optimum conditions were defined for Temperature (500-600 Celsius), reaction side pressure (16-40 bars), membrane thickness (1-7 micrometers), steam/methane ratio (3-4), reactor length (5-10 meters) and permeate sweep ratio (20 or more). This model was validated through multiple recognized sources (experimental and models). Adiabatic simulations were conducted in order to develop a mathematical model base for nonisothermal simulations. The results are showing how the energy consumption by reaction is happening at the beginning of the reactor. The maximum conversion of methane for the PBR is quickly reached while the conversion evolution observed for the membrane reactor is smoother. The convection through the membrane is high and the reaction and permeate temperatures are similar. Finally, the heat transfer between the exterior and the reactor has been considered in order to conduct non-isothermal simulations the overall heat transfer coefficient for the exterior/reactor convection energy transfer has been considered constant (UOR=817kJ/(m2.h.K)). Results showed high methane conversion (X>0.95) reached for high pressure (P>16bars) and medium temperature (Twall=Tinlet=600Celsisus), with the parameters range already indicated.
2
Table of Contents Acknowledgement.................................................................................................................. 1 Abstract.................................................................................................................................. 2 Figures list .............................................................................................................................. 4 Tables list ............................................................................................................................... 6 Introduction ........................................................................................................................... 7 Main problem......................................................................................................................... 8 1D Isothermal Steady-‐State model ....................................................................................... 12 PBR analysis ................................................................................................................................... 12 Fixed parameters................................................................................................................................12 Kinetic .................................................................................................................................................14 Mass balances.....................................................................................................................................14 MR analysis.................................................................................................................................... 15 Hydrogen permeability.......................................................................................................................15 Differential system .............................................................................................................................16 MATLAB model .............................................................................................................................. 17 Main program.....................................................................................................................................17 Functions and results representations ...............................................................................................18 Convergence issues ............................................................................................................................19 Results ........................................................................................................................................... 19 Effect of reactor length.......................................................................................................................20 Effect of total pressure .......................................................................................................................23 Effect of pressure drop .......................................................................................................................26 Effect of membrane thickness ............................................................................................................27 Effect of H2O ratio..............................................................................................................................30 Effect of Eta ........................................................................................................................................32 Effect of alpha.....................................................................................................................................34 Effect of sweep ...................................................................................................................................37 Comparison with sources ...................................................................................................................40 Optimum parameters .........................................................................................................................46 1D Adiabatic Steady-‐State model.......................................................................................... 49 PBR analysis ................................................................................................................................... 49 Heat transfers .....................................................................................................................................49 Heat balance.......................................................................................................................................51 3
MR analysis.................................................................................................................................... 53 Heat transfers and shell/tube problem ..............................................................................................53 Heat balances .....................................................................................................................................55 MATLAB model .............................................................................................................................. 59 Changes from isothermal model ........................................................................................................59 Functions and changes of volumetric flow rate .................................................................................60 Results ........................................................................................................................................... 61 PBR adiabatic ......................................................................................................................................61 Comparison PBR isothermal vs Adiabatic ...........................................................................................66 MR Adiabatic ......................................................................................................................................67 Comparison adiabatic MR vs adiabatic PBR .......................................................................................70
Non-‐isothermal .................................................................................................................... 78 Differential system......................................................................................................................... 78 U overall ........................................................................................................................................ 79 U overall constant .......................................................................................................................... 81 Temperature along the reactor, Tw=T0=500Celsius ..........................................................................82 Temperature along the reactor, Tw=T0=600Celsius ..........................................................................85 Conclusion ............................................................................................................................ 88 References............................................................................................................................ 89 Appendices........................................................................................................................... 90 MATLAB Code – Isothermal Main program..................................................................................... 90 MATLAB Code – Isothermal Functions ............................................................................................ 94 MATLAB Code – Adiabatic Main Code ............................................................................................ 95 MATLAB Code – Adiabatic Functions .............................................................................................. 97
Figures list Figure 1 : Packed Bed Reactor drawing ......................................................................................................10 Figure 2 : Membrane Reactor drawing.......................................................................................................10 Figure 3 : Effect of length, 350 C ................................................................................................................21 Figure 4 : Effect of length, 500C .................................................................................................................21 Figure 5 : Effect of length, 750C .................................................................................................................22 Figure 6 : Effect of length on Delta.............................................................................................................23 Figure 7 : Effect of pressure, MR vs PBR.....................................................................................................24 Figure 8 : Effect of pressure on Delta .........................................................................................................24 Figure 9 : 3D plot of Delta=f(Pt,T) side a.....................................................................................................25 Figure 10: 3D plot of Delta=f(Pt,T) side b ...................................................................................................26
4
Figure 11 : Effect of thickness, 350C...........................................................................................................28 Figure 12 : Effect of thickness, 527C...........................................................................................................29 Figure 13 : Effect of thickness, 700C...........................................................................................................29 Figure 14 : Effect of m at Pt=30bars, MR vs PBR ........................................................................................31 Figure 15 : D=f(T,m) at Pt=30 bars..............................................................................................................31 Figure 16 : D=f(T,Pt) in our model ..............................................................................................................33 Figure 17 : D=f(T,Pt) for Eta=1 ....................................................................................................................33 Figure 18: MR vs PBR, effect of alpha, Pt=30 bars......................................................................................35 Figure 19: MR vs PBR, effect of alpha, Pt=2 bars........................................................................................35 Figure 20: Delta=f(T,alpha) .........................................................................................................................36 Figure 21 : effect of sweep MR vs PBR, 350C .............................................................................................37 Figure 22: D=f(Pt,s), 350C ...........................................................................................................................38 Figure 23: effect of sweep MR vs PBR, 500C ..............................................................................................38 Figure 24: D=f(Pt,s), 500C ...........................................................................................................................39 Figure 25: effect of sweep MR vs PBR, 600C ..............................................................................................39 Figure 26: D=f(Pt,s), 600C ...........................................................................................................................40 Figure 27: Comparison MR vs PBR, Engin vs our model .............................................................................41 Figure 28: D=f(T,Pt), Engin vs our model ....................................................................................................41 Figure 29: PFR CH4 conv. vs T, Barbieri vs Shu vs our model .....................................................................42 Figure 30: MR CH4 conv. vs T, Barbieri vs Shu vs our model......................................................................43 Figure 31: PFR CH4 conv. vs T, Gallucci vs Shu vs our model .....................................................................44 Figure 32: PFR vs MR, CH4 conv. vs T, Gallucci vs Shu vs our model..........................................................45 Figure 33: Hydrogen flux for membrane reactor, effect of Pt and alpha ...................................................47 Figure 34:PBR Adiabatic, CH4 conv. vs inlet Temperature .........................................................................61 Figure 35:PBR adiabatic, Exit Temperature vs Inlet temperature ..............................................................62 Figure 36: Effect of pressure, Temperature vs ksi ......................................................................................63 Figure 37: PBR adiabatic, CH4 conv and Temp vs ksi, Pt=15 bars, T = 700C...............................................64 Figure 38: PBR adiabatic, CH4 conv and Temp vs ksi, Pt=30 bars, T = 700C...............................................64 Figure 39: Temperature vs ksi, at different inlet temperature...................................................................65 Figure 40: iso vs adia, effect of pressure in a PBR ......................................................................................66 Figure 41: CH4 conv. vs inlet temperature, effect of pressure for MR.......................................................67 Figure 42: Outlet Temperature vs inlet temperature, effect of pressure for MR.......................................68 Figure 43: MR adiabatic, CH4 conv and Temp vs ksi, Pt=30 bars, T = 700C................................................69 Figure 44: MR adiabatic, permeate and reactor temperatures vs ksi, T=700C and Pt=30bars ..................70 Figure 45: PBR vs MR, T out vs T in, effect of pressure ..............................................................................71 Figure 46: MR vs PBR, CH4 conv. vs T, effect of pressure ..........................................................................72 Figure 47:PBR vs MR, high pressures, conversion vs temperature ............................................................73 Figure 48:PBR vs MR, high pressures, out temp vs temp ...........................................................................73 Figure 49: Delta vs T in, effect of pressure (2 and 30 bars) ........................................................................74 Figure 50:Delta vs T, effect of pressure (15 and 30 bars) ...........................................................................74 Figure 51:Effect of catalyst deactivation on Delta, Pt=30 bars...................................................................75 5
Figure 52:Effect of catalyst deactivation on methane conversion, Pt=30 bars ..........................................76 Figure 53:Effect of catalyst deactivation on Temperature, Pt=30 bars ......................................................76 Figure 54:Non-‐isothermal MR and PBR reactors, Temperature vs ksi, Pt=30bars .....................................82 Figure 55:Non-‐isothermal MR and PBR reactors, CH4 conversion vs ksi, Pt=30bars .................................83 Figure 56:Non-‐isothermal MR and PBR, Temperature vs ksi, Pt=16bars ...................................................84 Figure 57:Non-‐isothermal MR and PBR, CH4 conversion vs ksi, Pt=16bars ...............................................84 Figure 58:Non-‐Isothermal MR and PBR, CH4 conv vs ksi ...........................................................................85 Figure 59:Non-‐Isothermal MR and PBR, Temperature vs ksi .....................................................................86 Figure 60:Non-‐Isothermal MR, Reaction and Permeate Temperature vs ksi, Tw=T0=600C, Pt=30bars ....87
Tables list Table 1: Independent parameters ..............................................................................................................12 Table 2: Kinetics constants .........................................................................................................................13 Table 3: CP coefficients ...............................................................................................................................50
6
Introduction
With the development of new technologies such fuel cells, the demand in hydrogen keeps on increasing. Methane Steam Reforming has been for the past decades and still is one of the major ways of producing hydrogen as fuel or for multiple chemical engineering processes. This tendency should not change in the future since the main reactants for the system Methane Steam Reforming (MSR) / Water Gas-Shift (WGS) are steam and natural gases, available easily in nature (compare to other energy sources). Therefore, research on this subject in order to optimize it is necessary. Numerous researches have been done on the MSR/WGS system and among them, membrane reactors are a promising technology in order to enhance hydrogen production. The use of Pd- and Pd/alloy- based catalytic membranes are showing remarkable results in terms of conversion and reactor performances, due to the hydrogen permeability through the Palladium. Limitations are however encountered with the membrane production for the thickness required, its behavior to impurities such as sulfur and the palladium price. In light of these considerations and with today’s computational processors, comprehensive process systems modeling frameworks are recognized as significant elements in this research field. The results obtained by these mathematical modeling methods help identifying the optimum system parameters such as reactor design, inlet composition and general state. They are also characterizing the inherently complex behavior of catalytic membrane reactors with a set of fes simplifications, giving advices for a future industrial realization of this process.
This work has been in order to answer the problems created by the process, such as: what simplifications have to be made? What is the influence of Temperature, pressure and inlet composition? Do the results are confirmed by experiments? Is the membrane reactor model presenting a realistic vision of the process? The mathematical modeling framework created to answer these questions has been develop on the software MATLAB © and is the heart of this thesis. By a set of simplifications, the process model has been accurately design, optimum parameters for the reaction have been determined and results have been compared to other models and experimental results. Our main source of comparison has been Engin Ayturk work and his isothermal 1D steady state model.
7
Main problem The system Methane Steam Reforming / Water Gas Shift appears as a simple model. Few species are involved as reactants – Steam and Methane – and the products are only Carbon Monoxide (which can react with steam), Carbon Dioxide and Hydrogen. However, multiple side reactions are in fact involved in the whole process and the kinetic associated to the system becomes difficult to use. The reactions impose a complexity of the system representation such as adsorption, carbon formation, and side reactions. From these different reactions, multiple species should be added in the process model. Xu and Froment (1989) worked on this problem in the 80ies and they tried to find the major reactions involved in this process, in order to neglect the others. In papers they published, they gave elements on which other researches have been done. From their results, three main reactions have to be considered in order to have the whole aspect of the process. Add to that, they define a specific kinetic for these reactions based on a Langmuir-Hinshelwood mechanism, making the task of modeling the process easier. These three reactions represent the methane steam reforming, the water-gas shift and the methanation:
We can see how the third reaction is the direct sum of the two previous one. This would be pushing us to eliminate it from our model. However, this set of reaction has been found from an intensive research from Xu and Froment (1989), and we will keep the reactions as they have been exposed. We will use in this thesis the following index:
Therefore, the associated reaction rates can be written: 8
PH3 2 .PCO ⎞ k1 ⎛ ⎟ DEN 2 r1 = 2.5 .⎜ PCH 4 .PH 2O − PH 2 ⎝ K1 ⎠ PH .PCO2 ⎞ k1 ⎛ ⎟ DEN 2 .⎜ PCO .PH 2O − 2 PH 2 ⎝ K 2 ⎠ PH4 .PCO2 ⎞ k1 ⎛ ⎟ DEN 2 r3 = 3.5 .⎜ PCH 4 .PH2 2O − 2 PH 2 ⎝ K 3 ⎠ r2 =
DEN = 1+ KCO .PCO + K H 2 .PH 2 + KCH 4 .PCH 4 + K H 2O .
PH 2O PH 2
€
The partial pressures are defined by the initial inlet, the conversion of methane (consumption), the conversion of carbon dioxide (creation), and the hydrogen flux through the membrane. These parameters are expressed by:
9
It is possible to draw both systems in order to have a visual representation of the situation:
Figure 1 : Packed Bed Reactor drawing
Figure 2 : Membrane Reactor drawing
The only but important difference between the two systems is the permeation zone, present in the Membrane Reactor at the other side of the membrane. Physically speaking, the permeate side is a tube of radius r0 inside the reactor of radius Ri. The membrane is covering the tube over a determined length L, which is the same length where the catalyst is added in the reaction side. Before and after this zone, only inert quartz is
10
filled in the reactor so the reaction is starting only in the defined zone. The catalyst used is a Nickel-based Alumina catalyst so that the equilibrium is quickly reached. In the tube side of the membrane reactor, a sweep gas is used in order to control the pressure and push out the hydrogen crossing over the membrane. The sweep ratio is an important parameter of our system. In both reactors, a Plug Flow analysis is used. We are considering an element of volume dV as our system. The area of the slice at x is a constant and we will take for the PBR the same area than the MR in order to get the same order in the results, allowing us to compare both reactors. This area is defined by:
It is interesting to notice that we can keep the same values for most of the parameters between the Packed Bed reactor and the Membrane reactor. The area of contact, the void fraction or the overall heat coefficient won’t be affected by the system we are working on. The initial conditions depend on the inlet gas injected, defined by its state: composition, pressure and temperature. The reactor parameters, such as length, diameter, void fraction or catalyst density will be the same between the two reactors. Therefore, in order the mass balances involved in the mathematical functions of our PBR model will be identical to the Membrane Reactor model’s ones, by considering in the Packed Bed:
If this is true for the Mass Balances, the situation becomes more complex for the heat balances. Therefore, the PBR situation will be separated from the MR situation in our modeling process. However, the isothermal mass balances will be used as a base to develop the adiabatic model.
11
1D Isothermal Steady-State model
The simplest model of our problem is one-dimensional, isothermal and steady state. By considering an isothermal reactor, we avoid the heat balances (which are difficult to represent accurately as we’ll see later), and by limiting the geometry in one dimension, the diffusion is considered as negligible. One will not forget that by taking these simplifications, the results obtained will be inaccurate and may be far from reality. However, the fact is that these models help us getting a range for experimental results, and depending on the simplifications made, this range is increased or limited. By suppressing the simplifications one after another, we will reach an accurate visualization of reality. This problem has to be always considered when using numerical models.
PBR analysis
The packed bed reactor is simpler to represent than the membrane reactor, therefore it will be our first model. Fixed parameters
Among all the parameters to consider, a lot of them can be considered as constant. Some are the reactor characteristics (table 1) while others are kinetic constants (table2). One important thing to notice is the effectiveness factor; element representing hwo difficult adsorption is for some molecules due to their size and the pores size. Values of Eta can change from 0.01 to 1 between papers, for all three reactions. In our case, we will follow the values given by Xu and Froment (1989) and Shu et al. (1994):
Table 1: Independent parameters symbol r0 Ri ε ρ cata Pp
12
value 2.54 6.35 0.5 2100 1
unit [cm] [cm] [-] [kg.m-3] [bar]
description Outer tube radius (without membrane) Inner shell radius Bed porosity Catalyst density Permeate side pressure
Table 2: Kinetics constants element/reaction CH4
H2O
CO
CO2 H2 reaction 1 reaction 2 reaction 3
constant Kabs0 ΔH G0 A G0 B G0 C Kabs0 ΔH G0 A G0 B G0 C Kabs0 ΔH G0 A G0 B G0 C G0 A G0 B G0 C Kabs0 ΔH k0 E η k0 E η k0 E η
value 1.8 e-1 -38.3 -75.3 7.6 e-2 1.9 e-5 0.4 88.7 -241.7 4.2 e-2 7.4 e-6 40.9 -70.7 -109.9 -9.2 e-2 1.5 e-6 -393.4 3.8 e-3 1.3 e-6 2.9 e-2 -82.9 1.8 e-4 240.1 0.01 7.6 67.1 0.3 2.2 e-5 243.9 0.01
unit [1/bar] at 823K [kJ/mol] [kJ/mol] [kJ/(mol.K)] [kJ/(mol.K2)] [-] at 823K [kJ/mol] [kJ/mol] [kJ/(mol.K)] [kJ/(mol.K2)] [1/bar] at 648K [kJ/mol] [kJ/mol] [kJ/(mol.K)] [kJ/(mol.K2)] [kJ/mol] [kJ/(mol.K)] [kJ/(mol.K2)] [1/bar] at 648K [kJ/mol] [kmol.bar0.5/(kgcat.h)] [kJ/mol] [-] [kmol/(bar.kgcat.h)] [kJ/mol] [-] [kmol.bar0.5/(kgcat.h)] [kJ/mol] [-]
13
Kinetic
Kinetic parameters are expressed by the Arrhenius law:
⎛ −Δ rG ⎞ K e,i = exp⎜ ⎟ ⎝ R.T ⎠ ⎧ Δ G = ∑ν .G i, j j ⎪ r i, j with : ⎨ € ⎪G = A + B * T + C * T 2 ⎩ j j j j
All of these parameters depend on the temperature (in Kelvin). In the isothermal case, the temperature is considered to be constant all over the length of the reactor and so equal to the inlet € temperature. Therefore, these parameters will be calculated for T0 and considered constant during the simulation.
Mass balances
From what we said in the introduction, only three reactions have to be considered in order to build our model. In a general situation, that would mean using three sets of differential equation. However, the reactions are not independent since the third one is the sum of the two others. Therefore, the system can be mathematically defined by two differential equations: one on the Methane conversion, the other on the Carbone Dioxide conversion. The mass balances corresponding to these two species are written and from them, the differential equations are found. As said in chapter 1, balances are considered on a volume dV, between x and x+dx and for a constant surface Ac:
14
In our MATLAB simulation, the conversion of methane will be the first variable and the conversion of carbon dioxide will be the second.
Therefore:
In our MATLAB simulation, the conversion of methane will be the first variable and the conversion of carbon dioxide will be the second.
MR analysis
Now that the PBR has been defined, we can build a model for the Membrane Reactor. The only difference for this reactor is the hydrogen exiting through the membrane. By decreasing the partial pressure of hydrogen in the retentate side, the equilibrium is pushed to the right and to an increase in the hydrogen production.
Hydrogen permeability
Numerous people investigated a mathematical expression for the hydrogen permeability. This parameter indeed depends on the quality of the membrane, its composition and can change 15
significantly in the case of an alloy. A general agreement is that an Arrhenius expression can be used in order to describe this variable. In our model, we consider a pure palladium membrane. Ayturk (2007) reported a comparison between experimental results found by multiple scientists and defined a general Arrhenius expression for QPd , within a 16% error margin:
QPd
⎛ −E 0 ⎞ = Q0 * exp⎜ ⎟ ⎝ R * T ⎠
⎧Q0 = 6.3227e −3 m 3 µm.m −2 .h.bar 0.5 ⎨ −1 ⎩ E 0 = 15630 J.mol
This expression has been used in all of our membrane reactor models.
€ Differential system
The flux of hydrogen through the membrane is described by the Sievert law:
And the differential equation resulting from this flux is defined by:
Therefore, by considering the parameter YH2 as the ratio of hydrogen flow rate over the initial methane flow rate, the third differential equation can be written:
dYH 2 dξ
=
QPd .2π .( r0 + δ ).L 0 FCH .δ .22,4 4
(
. PH0.52 ret − PH0.52 perm
)
Note: the number “22,4” is used to get the correct units in our equation.
€ The final differential system is now:
16
MATLAB model
We can now write our MATLAB model. Since our goal is to compare the PBR with the MR, our model needs to simulate both reactors, stores the results and plots a visual comparison. However, by representing conversion over length, we can see the Main program
The main program is separated in parts. The first part defines the constant parameters. However, we need to separate these constants in three categories: the ones depending on Temperature, the working parameters and the others. Our first part will be explicit equations of these last ones – constant parameters independent of Temperature and non-working parameters. The second part defines the working parameters:
The Temperature is a specific working parameter, due to the kinetic aspect of the reactions. Therefore, its value will usually be defined between 350 and 850 degree Celcius – 500 17
different values. In order to express this, the initial temperature is given and a loop will increase the temperature by one degree for 500 points. On the other hand, the other parameters are associated with only a few different values since we just want to see their general influence on our system. No more than six values are therefore given to each parameter. In terms of coding, a six-steps loop is including inside the temperature loop, associated to 6 values for these parameters. Two explicit equations are then written for each parameter: a unique value (if we keep it constant) and a 6 elements vector (if it’s our working parameter). This allows us to modify the final results, if we want to see the influence of pressure and methane/steam ratio for a constant temperature for example. The third part is inside the loops. First are calculated the kinetic values: adsorption constants, reaction rate constants, equilibrium constants and Hydrogen permeability. These parameters depend on temperature; therefore we are calculating them only now since T is one of the main parameters of our system. Also calculated here two constants defined as:
These two constants are calculated in the main program and not in the functions, avoiding them to be calculated at every step of the numerical resolution and saving solving time. Then, the increment is defined and we solve our differential equation systems with MATLAB functions (ode15s in our case). The conversions of methane, carbon dioxide and hydrogen are stored in matrix. The last part is the graphic representation of the results depending on the parameters.
Functions and results representations
We use ode15s to solve our differential equations, which is a 4th order Runge Kutta algorithm. This will give us a good accuracy for the solvation. The first things to write in the functions are the expressions for the partial pressures since they are functions of the conversion. Then, the reaction rates are calculated with these partial pressures, which allow the differential system to be written numerically.
18
An interesting point is the fact that for the iteration n, the reaction rates are calculated from partial pressures defined by n-1 conversions. Mathematically speaking, this is a correct approach to solve a differential problem. Convergence issues
In order for our system to converge, we need to avoid any impossible numerical value. In our model, this problem can appear due to the expression of our rates of reaction. Indeed, a pressure of hydrogen equal to 0 will make the convergence impossible to achieve. Since our system is producing hydrogen, this problem has to be considered only initially. The solution we used to avoid this problem is to consider the initial hydrogen flow rate different from 0. This is a mathematical representation, not a physical one; therefore the value taken for this initial hydrogen flow has to be small enough to be negligible in front of the other flow, but large enough to avoid divergences issues in our mathematical model. Shu et al. (1994) showed that accurate results are obtained with this method as long as we use a small iteration step; they even produced a table linking the precision with this ratio. Multiple values of precision and initial hydrogen ratio have been tried and we chose the highest values possible in order to avoid long simulations: In our model, we use:
θ H2 =
FH0 2 0 FCH 4
= 0.01 ; ε (ξ) = 1.e −4
Isothermal simulations took usually 5 minutes (for 500 Temperature points-increments, 5 pressures and 4 different alpha). €
Results
As said above, multiple parameters are involved in the optimization of our models. Pressure, Temperature and composition of the inlet flow of gas are obviously capital but the dimension of the reactor and thus the membrane are also important. In order to determine the optimum parameters of our system, we’ll consider the variation of only one variable while the others remain constant. The temperature and the dimensionless length are specific parameters however. Indeed, the temperature dependence of the kinetic parameters (adsorbance, reaction, equilibrium, permeation) imposes the system to consider the 19
T-variation as a description of our results. The same can be said for
since the system is 1
dimension dependent. In order to compare the Packed Bed reactor with the Membrane reactor, the parameter Delta will be used. It is defined by:
This parameter is representing the efficiency of the membrane reactor over the packed bed reactor; therefore high values of Delta would indicate the state where our membrane reactor is the most efficient. Delta should go from -1 to 1, however the membrane reactor is showing better conversion for any working parameter changes, therefore Delta goes from 0 to 1. Effect of reactor length
In a first time, we’ll consider the variation of the reactor size i.e. its length since we are in 1D.
The characteristics of our model are:
Three different temperatures [350 ; 500 ; 750] and six pressures [2 ; 5 ; 10 ; 20 ; 30] have been considered (T in Celsius degrees, P in bar). The membrane reactor conversion is indicated in straight lines while the packed bed reactor conversion is in dotted line.
20
Figure 3 : Effect of length, 350 C
Figure 4 : Effect of length, 500C
21
Figure 5 : Effect of length, 750C
The first element we can notice is the difference of conversion between the PBR and the MR. After around two meters, the conversion of methane in the membrane reactor is higher than in the case of the PBR. Another point is that the maximum conversion of methane is obtained faster for the PBR than for the MR, due to the constant flow of hydrogen through the membrane. However, the conversion obtained is higher in the MR case and depends on the total pressure, but we will talk about this later.
The interesting thing we can see from these results is that while working at high temperature (750degC), the maximum conversion is obtained for a dimensionless length equal to 0.2, equivalent to a length of 2.5 meters. Even at 500degC, the methane is totally consumed for high pressure (in the case of the MR) after a length of 10 meters (ksi equal to 0.8). Thus the reactor length (and the membrane length) can be optimized depending on the operating temperature. Since the reactor price is for a large part due to its membrane, a decrease of its length will be an interesting point.
22
Figure 6 : Effect of length on Delta
Comparison between MR and PBR methane conversions shows how by increasing the length of the reactor large values of Delta are reached for low temperature. We can also see that a reactor length of 5 meters allows us to obtain a high delta for a temperature of 560degC.
As a result, we will consider a length of 7 meters for our next simulations.
Effect of total pressure
Let’s now consider the effect of pressure in our system. As before, the dotted lines are representing the PFR while the plain lines are representing the MR. The characteristics of our 0 ⎧ L = 7 m;δ = 5 µm;m = 3;FCH = 1kmol.h −1 ;s = 10 4 ⎪ representations are: ⎨ P ∈ [2 : 5 :10 : 20 : 30] ⎪ ⎩T ∈ [ 350 : 750]
€
23
Figure 7 : Effect of pressure, MR vs PBR
Figure 8 : Effect of pressure on Delta
24
The effect of pressure on the MR is the opposite of on the PBR. When the pressure is increased, the conversion in the Packed Bed reactor is decreasing while for the Membrane reactor, this conversion is increasing with the pressure. These results are confirmed by the literature and were expected when we consider the differential equation on the hydrogen flux (depending on the hydrogen pressure). This is one of the main reasons of using a membrane. By using a graphic representation of Delta, we can see the link between these two parameters. Indeed, the higher the pressure, the higher Delta. It also appears that these high values of Delta are obtained for medium temperature. In order to have a larger visualization of the effect of pressure and temperature, a 3D plot presents interesting things. In order to create this representation, a temperature variation of five degrees and a pressure variation of one bar are used:
Figure 9 : 3D plot of Delta=f(Pt,T) side a
25
Figure 10: 3D plot of Delta=f(Pt,T) side b
The optimum state pressure/temperature is easily observable on this representation (the reddest parts are the highest delta). The 3D results are showing a value of Delta superior to 0.8 for T=[500,600] and P=[16,40]. These data are important since they are giving us a range of values for the operating pressure and the temperature for optimizing the others parameters.
Effect of pressure drop
Considering the length of our system, the first thought was to evaluate the variation of pressure due to the pressure drop. However, the calculations are showing no influence in our system so the pressure drop evaluation has not been considered any longer in our simulations. As a reminder, we will write how the pressure drop has been calculated. The equation is a typical Ergun expression of the differential variation of pressure:
26
Multiple parameters are considered in this expression, and have to be calculated at each step since they depends on the characteristics of the system at the point considered: Gx is the mass specific gas flow rate, defined by
Rep is the modified Reynolds Number defined by
The mixture density is calculated as followed:
Epsilon is the void fraction of our bed The mixture viscosity is a complex parameter and its calculation referred to Poling, Parausmitw and O’Connell’s book “the properties of Gases and Liquids”
with
The different viscosities have been calculated with the following equation:
The different parameters are expressed for each species in Poling et al (2007) book. As we said earlier, these complex mathematical calculations have no incidence on our system since the pressure drop is negligible. They are therefore not considered in our system, which is reducing the processing time to solve our problem. The only interesting thing is that the mixture viscosity expression is used in heat transfer coefficient correlations. We will see this in the non-isothermal chapter.
Effect of membrane thickness
The next parameter to study is the membrane thickness. We are representing Delta for three different temperatures (350, 527 and 700 Celsius degrees), six different pressures (1, 2, 5, 27
10, 20 and 30 bars) and the membrane thickness value is starting from 1 micrometer and goes to 30 micrometers:
Figure 11 : Effect of thickness, 350C
28
Figure 12 : Effect of thickness, 527C
Figure 13 : Effect of thickness, 700C
29
Two different behaviors appeared from our results. At low temperature (350degC), an increase of the membrane thickness will provoke an immediate decrease of Delta whereas at high temperature (700degC), the ratio Delta barely changes with this parameter. This can mathematically be explained by the expression of the hydrogen transfer through the membrane, since at low temperature the expression of QPd is small and therefore the total expression is highly dependent on the membrane thickness. We saw previously that the maximum of Delta was obtained for 527Celsius degrees (for a pressure of 30 bars); this is why the membrane thickness effect has been plotted for this temperature. From the graphic representation, we can see the how a membrane thickness of five micrometers is giving a maximum for Delta as long as the pressure is superior to 20 bars. Therefore, this thickness will be conserved in our simulations. It is interesting to keep in mind that the reactor price is due for a large part to the membrane. But add to the fact our results are for a specific case (isothermal, 1D, steady state), there are physical limitations in order to build membranes. A perfect membrane would be extremely thin, but building it might be technologically difficult. This is another problem however, that cannot be included in our simulations.
Effect of H2O ratio
As said above, methane and water are the only species introduces in the inlet gas, hydrogen is present only to avoid mathematical divergences. Therefore, the ratio steam/methane (called m) is an important parameter of our system since it is defining this inlet flow. The flow of methane is a constant in our model; therefore an increase of this m means an increase in the flow of steam used in our reactor. Until now, the ratio has been considered as 3 for our simulations, value used in different sources. We will see now the influence of m on our system:
30
Figure 14 : Effect of m at Pt=30bars, MR vs PBR
Figure 15 : D=f(T,m) at Pt=30 bars
31
The results are showing a similar behavior for every cases, except for m=1. At a ratio steam/methane superior to 3, the membrane reactor isn’t presenting any change while the methane conversion is increased for packed bed reactor, due to the equilibrium definition (ratio of partial pressures). These graphic representations are confirming our choice of taking 3 for this ratio, since this value is giving the best available results. Taking m=2 is giving the same kind of conversion but for a slightly superior temperature, therefore it is understandable than some papers are using this ratio. In our case however, we’ll try to obtain the best parameters for a low temperature.
Effect of Eta
As we said above, the effectiveness factors Eta have to be taken in account in order to obtain an accurate representation of our system. Physically, these parameters are characterizing diffusion through the pores and micro-pores and adsorption constraints on the catalyst. These elements depend on temperature and so does each effectiveness factors. However, the Eta values are considered as independent of temperature in multiple sources, and instead are associated with typical values for each reaction, values given by Xu and Froment (1989):
. It is interesting to notice the effectiveness
factors for reaction 1 and 3 are equal. Considering the molecular equations, this should means that the methane is adsorbed on the catalyst and this adsorption is the limiting parameter for Eta (value smaller for the first and third reactions, where the methane is involved). Let’s see how this parameter affects our system:
32
Figure 16 : D=f(T,Pt) in our model
Figure 17 : D=f(T,Pt) for Eta=1
33
The first plot represents Delta when we take in account the values of Eta given previously. The second plot shows the results for Eta assimilated as one (perfect effectiveness). At high temperature (superior to 600degC), the results appear to be the same. At low temperature however, a non-negligible variation is observed. The results are showing a difference up to 50% on Delta at these temperatures. In the state we are interested in (highest delta); the maximum value is obtained for a twenty degrees difference. For 30 bars, this maximum (almost 0.9) is obtained at 485degC while considering Eta, whereas we can reach it at 460degC without considering the effectiveness factor. This is proving how values of Eta have important consequences on our model, and that this parameter has to be considered in order to obtain more accurate results. However, we have to keep in mind that Eta is in fact a function of temperature; therefore considering Eta as independent of it is a simplification that should be avoid in a complete model.
Effect of alpha Until now, we considered the catalyst to be stable in any situation. Eta is expressing the difficulties molecules have to reach the catalyst only. Physically speaking, the catalyst efficiency is decreasing over time due to different reasons. In our case, since the methane appears to be one of the species adsorbed, it can react in a side reaction and form Coke on the catalyst. These Carbon structure will bind to catalyst sites and decrease its efficiency. The sites can also be poisoned by impurities in our inlet flow, such as sulfur. Whatever the physical reason(s) for this, we can represent this decrease in efficiency with a parameter alpha. Since all reaction are happening on the catalyst surface, this parameter is the same for all three reactions and is integrated to our mathematical model as followed:
Since our model is at steady state, we will simulate it for different values of alpha (but with all other parameters constant). Let’s represent the methane conversion for both systems due to this parameter:
34
Figure 18: MR vs PBR, effect of alpha, Pt=30 bars
Figure 19: MR vs PBR, effect of alpha, Pt=2 bars
35
As we can see, the catalyst deactivation is having different effects on the system, depending on the temperature and pressure we are working on. For the packed bed reactor, variation of a few degrees is observed at low temperature, when the conversion of methane is low. The effect is more visible at low pressure than at high pressure, which is in accordance to the results obtained for the methane conversion of the PBR observed previously. For the membrane reactor, the variation of conversion is higher than for the packed bed reactor. This variation seems to be independent of the temperature; it just decreases at low temperature (when the conversion is very low) or at high temperature (when the conversion reaches 1). The parameter Delta is plotted for low and high pressure:
Figure 20: Delta=f(T,alpha)
The general behavior of Delta is unchanged. However, the maximum value obtained decreases with alpha, and it is reached for higher temperatures. Even at half efficiency, we can still reach a Delta of 0.86 at 30 bars and 555 Celsius, so in the range of temperature and pressure defined previously.
36
Effect of sweep
The sweep gas is used only in the membrane reactor. It is used to help the hydrogen being pushed out of the tube, allowing us to control the pressure in the permeate side. For the isothermal situation, the composition of this gas isn’t important.
Figure 21 : effect of sweep MR vs PBR, 350C
37
Figure 22: D=f(Pt,s), 350C
Figure 23: effect of sweep MR vs PBR, 500C
38
Figure 24: D=f(Pt,s), 500C
Figure 25: effect of sweep MR vs PBR, 600C
39
Figure 26: D=f(Pt,s), 600C
From these graphic representations, we can see the influence of the sweep ratio on our system. At low temperature (350C), the methane conversion is continuously increasing with this ratio. As a result, Delta is behaving in the same way since the sweep ratio has no incidence on the Packed Bed reactor case. We can also notice how the general behavior of the plots is identical for any pressure. At medium and high temperature (500 and 600C), two things can be observed. First of all, the high pressures are at their maximum conversion/Delta immediately. And secondly, when this ratio is above 20, all pressures are showing a conversion/Delta equal to their maximum (difference of less than 5%). Since we realized that the best performances of our membrane reactor are obtained for a temperature higher than 500 Celsius, a sweep factor of 20 will be chosen as optimum ratio.
Comparison with sources
Multiple works have been done on reactor modeling for the methane steam reforming. It is however impossible to compare all of them on the same graphic representation since every
40
results published have been found for a different set of conditions. Therefore, we will use these conditions in our model in order to make 1/1 comparison. The first work we will use as a comparison is Engin Ayturk (2009) work, since it was the base of our research. The parameters he used for his research are the same than us in order of reactor size, kinetics and membrane characteristics. The main difference results in the value of Eta, assimilated as 1 for his work:
Figure 27: Comparison MR vs PBR, Engin vs our model
Figure 28: D=f(T,Pt), Engin vs our model
41
As we can see, the results are basically similar, except our model is taking in account the effectiveness factor. The slight difference observed is in the same order of the difference we saw when we studied the effect of Eta.
Others results to consider are from Barbieri and Di Maio (1997) and Shu et al (1994) works. Barbieri and Di Maio (1997) were first comparing their model to the experiments made by Shu, Grandjean and Kaliaguine (1994) and they obtained accurate results. The conditions they used are, however, very different than ours:
Let’s compare their results to our model:
Figure 29: PFR CH4 conv. vs T, Barbieri vs Shu vs our model
42
It appears from this comparison that our model is not following the experimental data (and Barbieri’s model) when temperature is increased. However, multiple parameters are missing in Barbieri’s and Shu’s papers, like the characteristic of the bed or the effectiveness factor. Also, considering our simplifications (one dimension, steady state, etc), a non-perfect correspondence between experiments and model was expected.
Let’s compare the membrane reactor results now. Barbieri et al (1997) simulations were for 3 different flow in the tube side (sweep factor), and his model shows parallel flow configuration (solid lines) and counterflow configuration (dotted lines). The symbols are representing experimental results from Shu et al (1994), for a sweep flowrate of 40 SCCM:
Figure 30: MR CH4 conv. vs T, Barbieri vs Shu vs our model
These results highlight the accuracy of our model. Indeed, we can see that our model is following accurately the experimental results at the contrary of Barbieri’s results. In his model, 43
the sweep flow rate has an important impact on the methane conversion whereas for us, the conversion is only slightly changing with the sweep ratio. It is interesting to notice that considering the initial flow rate of methane (162 SCCM = 0.0476 kmol/h for 500C and 1.36 bars), the sweep ratio is really low compared:
⎧ 40 0.7 ⎫ ⎪ ⎪ ⎨80 SCCM ⇔ sweep ratio 1.45 ⎬ ⎪ ⎪ 3.6 ⎭ ⎩200 These values are lower than the ratio used in our optimum model (sweep=20).
€ Interesting works have also been done by Gallacci et al. (2004) in isothermal modeling. As did Berbieri et al. (1997), they compared their results with the experiments made by Shu et al. (1994) The inlet conditions are a little different than what we just saw but the reactor and the membrane characteristics are the same:
Figure 31: PFR CH4 conv. vs T, Gallucci vs Shu vs our model
As we can see, our PFR model is extremely close to Gallucci’s and the experimental results. Indeed, the maximum error is at 300 Celsius when the conversion obtained is 0.05 when
44
Shu et al. got 0.03. This is giving some validity to our isothermal PBR model (and is giving us some doubts on the conditions indicated in Barbieri’s paper).
Let’s now compare the results obtained by Gallucci et al in the membrane reactor case. Concerning the conditions, one of the major differences is the membrane thickness, equal to 50 microns in this case. The results has been put on the same graphic representation, with TR equivalent to Traditional Reactor (PBR):
Figure 32: PFR vs MR, CH4 conv. vs T, Gallucci vs Shu vs our model
The difference between the two models is more visible in this plot. Our model is showing a behavior equivalent to the experiments and Gallucci et al. (2004) model, however a maximum variation of 0.1 is observed at 400 Celsius. The variation between models can be explained by the conditions used and some parameters that have not been indicated. We saw how the effectiveness factors values were affecting our system at low temperatures; this can be a reason of this difference. The kinetics constants chosen can also been an issue if the sources used are different. Finally, Gallucci et al. (2004) expressed molar differential variations for every species and therefore, used a differential system of seven equations. Considering the step of integration used and the different parameters, this can explain the variation. 45
The difference between our model and the experimental results are probably the results of our simplifications, especially the one-dimension case (since the reactor is small, the isothermal simplification can be considered as accurate). Also, the results are made of only a few points; therefore the error may not have been reduced by repetition.
In conclusion, different sources are partly confirming the accuracy of our model and allow us to use it as a base to develop an adiabatic model.
Optimum parameters
The isothermal has been studied for multiple parameters. We can now define optimum parameters. We have to be careful when using the term “optimum”. Indeed, there are multiple ways of considering how good a reactor is. In our situation, it can be from the quality of the hydrogen produced (high quality avoid a separation unit after the reactor), the quantity of hydrogen produced, the energy required to work, material used, etc.
The first economical element to consider is the membrane. Indeed, Palladium is a very expensive material; reducing the size of this membrane will decrease the total initial price of our reactor. In our situation, the membrane thickness and length have been considered. We saw that a high conversion of methane can be obtained as long as the length is higher than 5 meters and for the thinnest thickness. Since obtaining a thin membrane is physically a challenge, a thickness of 5 microns is the average value of our optimum thickness. Concerning the length of the membrane, a smaller length means a higher operative temperature in order to produce the same results. Therefore, the economic comparison between initial cost and operative cost should be compared in order to obtain the optimum operating temperature and membrane length. We are considering only the range of these parameters.
Concerning the pressure, we saw that high conversion was obtained for high values but the operating cost of compressors is linked to this. Therefore, the optimum pressure of our system should be the lowest possible which can still give us high methane conversion. On this point, we can observe the effect of pressure (for an alpha value of 0.5 and 1) at different pressures: 46
Figure 33: Hydrogen flux for membrane reactor, effect of Pt and alpha
For a methane inlet flux of 1 kmol per hour, we can see a production of nearly 4 kmol per hour of hydrogen when the pressure is superior to 16 bars, even when the catalyst deactivation value is taken equal to 0.5 (T=575Celsius in this case). The optimum pressure is therefore confirmed to be in the range of pressure defined previously.
Concerning the steam/methane ratio, since the methane flux has been taken constant in our model, it represents only the steam flow. The system temperature is high enough to boil the water into steam. If we are using steam as a sweep, we need to condensate it after the reactor; therefore in heat exchanger will be used in order to save in energy. The only economical parameter to consider for the steam is the compressors and pumps; therefore a low ratio will decrease the operating cost. We saw that the total conversion stays the same for ratio superior to 3 therefore this value will be taken as our optimum parameter. Linked to this is the optimum flow of sweep. Considering that this ratio is giving a sweep flow value of multiple times the flow of methane, heat exchangers and compressors are needed for high value of sweep. In order to get a high conversion of methane and limit these operating costs, the optimum sweep ratio has to be the lowest value of the range found previously, i.e. 20.
47
The conclusion of the different simulations we did on all parameters is focused in the following optimum parameters for our reactor:
48
1D Adiabatic Steady-State model
The isothermal 1D steady-state model is a very simple representation of our system. One of the most important simplifications is to avoid considering the heat balance. However, our set of reactions shows an endothermic effect. Since multiple kinetics parameters are present, a change on the temperature should have large consequences on the results. This is why the next step in our model has to be the non-isothermal case. One major difficulty by considering the heat balance is the heat transfer coefficients values. Indeed, the energy parameters to consider for our system are: -
Heat of the flux Heat produced by the reaction Heat of the hydrogen flux going through the membrane Heat exchange between shell and tube Heat exchange between the shell and outside
Physically speaking, an adiabatic model is useless since it will never be build due to the system endothermic behavior. Mathematically speaking, an adiabatic model is a way to express the heat balance without taking in account the difficult expression of the heat transfer coefficient linked to the exchange between the shell and outside (furnace). It is therefore a capital step in our model.
PBR analysis As we did for the isothermal case, we will start with the PBR model. Heat transfers
Multiple parameters are involved in the heat transfers.
First of all, any change of temperature will involve a change in the kinetic parameters’ values. As a result, these kinetics parameters have to be re-calculated in each system dV considered. 49
Furthermore, each specie flow is bringing a heat to the system, heat expressed with its molar flow rate and its heat capacity. These Cp are functions of temperature and have to be calculated for each system, as the others kinetic parameters. Expressions used for Cp are depends on empirical constants and usually characterize an ideal gas. In our situation, the temperature and pressure conditions are low enough to be considered each gas as in its ideal state, allowing us to use this expression:
Cp j
Dj −3 2 −6 5 = A + B * T *10 + C * T *10 + j j j 2 *10 R T
With: Table 3: CP coefficients
€ A
B (103)
C (106)
D (10-5)
CH4
1.702
9.081
-2.164
0
H2O
3.470
1.450
0
0.121
CO
3.376
0.557
0
-0.031
CO2
5.457
1.045
0
-1.157
H2
3.249
0.422
0
0.083
As we just said, the species energy is defined by:
(
E j = Cp j * F j * T − Tref
)
In order to express each flow rate and considering the gas at their ideal state, we are using the law for ideal gas in the case of a PFR: € P *Q Pi * Q = Fi * R * T ⇔ Fi = i R*T As a function of conversion, the partial pressures are defined at each step, as for the temperature. The volumetric flow however has not been defined yet. In order to obtain its € expression, we have to remember that the pressure in our system is constant. Therefore:
This volumetric flow rate has to be calculated at each iteration since the molecular dilatation is affecting every flow rate, so the energy balances. The only available information is its initial value, defined by the inlet composition and state.
50
The last element to take in account in our heat balance is the energy produced by the reactions. It depends on the enthalpies for all three reactions and the rate of reaction for each case, since the rate of reaction are expressed as function of the conversion. The Enthalpy values at reference temperature (298K) are:
Δ r H 298K , reaction 1 = 206310 J.mol −1 Δ r H 298K , reaction 2 = −41200 kJ.mol −1 Δ r H 298K , reaction 3 = 165110 kJ.mol −1 The Enthalpy is a function of the temperature T: € dΔ r H i = ∑ν i, j * C p ij (T ).dT j T Δ r H i (T ) = Δ r H 298K + ∑ν i, j *
i Pj
∫ C (T ).dT
298K
j
T ⎛ ⎞ Dj Δ r H (T ) = Δ r H 298K + ∑ν i, j * ∫ ⎜ A j + B j * T *10 −3 + C j * T 2 *10 −6 + 2 *10 5 ⎟.dT T ⎠ 298K ⎝ j €
⎡ ⎛ 1 Bj Cj 1 ⎞⎤ Δ r H (T ) = Δ r H 298K + ∑ν i, j * ⎢ A j .(T − 298) + *10 −3 * (T 2 − 298 2 ) + *10 −6 * (T 3 − 298 3 ) − D j *10 5 * ⎜ − ⎟⎥ 2 3 ⎝ T T298 ⎠⎦ ⎣ j
€
€
We can now integrate in our heat balance the enthalpy for each reaction, and only as function of the temperature.
Heat balance
Now we defined all energy characteristics of our system, it is possible to write the heat balance on dV:
⎛ ⎞ Cp * F * T − T + −Δ H T * η * r * ρ ( ) ⎜ ∑ j j ∑ ( r i ) i i cata ⎟ * AC * L * dξ ref ⎝ i ⎠ j
(
)
(
) (
= ∑ Cp j * F j + dF j * T + dT − Tref j
(
)
) (
)
= ∑ Cp j * F j * T − Tref + ∑ Cp j * dF j * T − Tref + ∑ Cp j * F j * dT + o( dT.dξ) j
€
j
j
51
Mathematically speaking, the differential of length multiply by the temperature has the same differential order than the differential of temperature multiply by the flow rate. However, the product of both differential terms (temperature and dimensionless length) is negligible in front of the rest of the equation. Therefore, this differential equation can be simplified in:
dT = dξ
⎛ ⎞ dF j * T − Tref ⎜∑ ( −Δ r H i (T )) * ηi * ri * ρ cata ⎟ * AC * L − ∑ Cp j * dξ ⎝ i ⎠ j
(
∑ Cp
j
(
* F j * T − Tref
j
)
)
The reference Temperature is 298K. The differential flow term is described depending on the differential length, its expression € obtained through mass balances for every species:
⎧ F − (η .r + η .r ) * ρ * A * L * dξ = F + dF 1 1 3 3 bed C CH 4 CH 4 ⎪ CH 4 ⎪ FH O − (η1 .r1 + η2 .r2 + 2 * η3 .r3 ) * ρbed * AC * L * dξ = FH O + dFH O 2 2 ⎪ 2 ⎨ FCO + ( −η1 .r1 + η2 .r2 ) * ρbed * AC * L * dξ = FCO + dFCO ⎪ ⎪ FCO2 + (η2 .r2 + η3 .r3 ) * ρ bed * AC * L * dξ = FCO2 + dFCO2 ⎪ ⎩ FH 2 + ( 3* η1 .r1 + η2 .r2 + 4 * η3 .r3 ) * ρbed * AC * L * dξ = FH 2 + dFH 2
⎧ dFCH 4 = −(η1 .r1 + η3 .r3 ) * ρ bed * AC * L ⎪ ⎪ dξ ⎪ dFH O 2 = −(η1 .r1 + η2 .r2 + 2 * η3 .r3 ) * ρbed * AC * L ⎪ d ξ ⎪ ⎪ dF ⇔ ⎨ CO = ( −η1 .r1 + η2 .r2 ) * ρ bed * AC * L ⎪ dξ ⎪ dFCO 2 = (η2 .r2 + η3 .r3 ) * ρbed * AC * L ⎪ d ξ ⎪ ⎪ dFH 2 = ( 3* η1 .r1 + η2 .r2 + 4 * η3 .r3 ) * ρbed * AC * L ⎪ ⎩ dξ
Since the rates of reaction have already been expressed depending on constants, temperature and conversion, a numerical solution is possible. The system of differential equations becomes:
52
MR analysis
Now that the Packed Bed reactor adiabatic model has been made, we can develop the Membrane reactor model. In order to represent the system accurately, we will be talking about the reactor side as the shell or reaction side, and the other side of the membrane (inner side of the reactor) the tube or permeate side.
Heat transfers and shell/tube problem
There are multiple new data to consider in the membrane reactor heat balance that were not taken in account in the case of the packed bed reactor heat balance.
First of all, the flow of hydrogen going through the membrane has to be considered. Linked to this mass transfer, the intrinsic energy hold by this flow is exiting the shell side and enters the tube side. This will increase the energy of the permeate side, so its temperature. This last point is indicating that both temperatures have to be considered for each side of the membrane. As a result, two heat balances have to be written and two differential equations will be part of our differential system, giving us a set of five differential equations.
53
The hydrogen flux through the membrane is defined by the parameter YH2 and the coefficients for the heat capacity of hydrogen have already been written above. Since the flux is going from the shell side to the tube side, the temperature of this flux is the reaction temperature:
(
dQH 2 = Cp H 2 (Treac ) * df H 2 * Treac − Tref YH 2 =
f H2 0 FCH 4
)
0 ⇔ df H 2 = dYH 2 * FCH 4
The second heat transfer to consider is the convection due to temperature difference € between the shell and the tube results. Physically speaking, we should be considering multiple temperatures: temperature of the reaction side, temperature of the limit layer on the membrane at the reaction side, temperature of the membrane, temperature of the limit layer on the membrane at the permeate side and temperature of the permeate side. However, DeFalco et al. (2007) highlighted the fact that thermal conductivities in metallic membrane structures are high, which leads us to neglect thermal resistance for the membrane. Add to the fact that hydrogen is going from the shell to the tube side, the membrane temperature is estimated equal to the reaction side temperature, and the limit layers temperatures are neglected. Therefore, the convection term is written:
(
conv dQmemb = hmemb * dS * Treac − Tperm
)
The surface depends on the differential dimensionless length. The heat transfer coefficient is a function of multiple parameters. However, based on € Madia et al. work, it can be considered as constant:
(
)
dS = 2π * ( ro + δ ) * L * dξ hmemb = 8640 J.m −2 .h −1 .K −1 In the way this energy term has been written, it is defined as an out-flux i.e. the energy leaving the system. However, if the tube temperature is higher than the shell temperature, this € term will be negative and correspond to an accurate energy transfer (heat is going from hot to cold). Since the convection term is a function of both temperatures (shell side and tube side), both energy balances are bound and have to be solved together at each integration step.
The last heat energy term would be the convection between outside and the reaction side through the outer wall of the reactor. Since we are in the adiabatic case, this term is equal to zero.
54
Heat balances
As we said above, two heat balances have to be considered on two different systems. Since they are made of multiple terms, we will first define the reactor temperature heat balance:
∑ Cp
reac j
j
⎛ ⎞ * F j * Treac − Tref + ⎜ ∑ −Δ r H i (Treac ) * ηi * ri * ρcata ⎟ * AC * L * dξ ⎝ i ⎠
(
)
(
(
) (
)
)
(
)
= ∑ Cp reac * F j + dF j * Treac + dTreac − Tref + hmemb * 2π * ( r0 + δ ) * L * Treac − Tperm * dξ j j
Let’s develop the first exit term:
∑ Cp
€
reac j
j
(
) (
)
(
* F j + dF j * Treac + dTreac − Tref = ∑ Cp reac * F j * Treac − Tref j j
+ ∑ Cp reac * j j
)
dF j * Treac − Tref * dξ dξ
(
)
+ ∑ Cp reac * F j * dTreac j j
+ ∑ Cp reac * dF j * dTreac j j
As we can see, the first term is equal to the energy flux-entrance term, so it can be simplified. € The second and third terms are first order differentials, respectively in temperature and in dimensionless length (since the dF are functions of dksi). The fourth and last term is a second order differential, thus it can be neglected from the equation. This gives us: ⎛ ⎞ ⎜∑ −Δ r H i (Treac ) * ηi * ri * ρcata ⎟ * AC * L * dξ ⎝ i ⎠ dF j = ∑ Cp reac * * Treac − Tref * dξ + ∑ Cp reac * F j * dTreac + hmemb * 2π * ( r0 + δ ) * L * Treac − Tperm * dξ j j d ξ j j
(
)
(
€
)
(
)
It is important to notice that the heat capacities involved in this heat balance have to be calculated at a specific temperature. We already wrote the heat capacity general expression, which is a non-linear function of the temperature. As a result, the heat capacities differential 55
variations due to the differential change of temperature are difficult to express, and only a part of it is a first order differential. In order to simplify our model, we will consider this variation negligible. Therefore, our first energy differential equation can be written:
We can now write the heat balance for the tube side:
[Cp
perm sweep
] (
)
(
) )
(
)
* Fsweep + Cp Hperm * FH 2 * Tperm − Tref + Cp reac H 2 * df H 2 * Treac − Tref + h memb * 2π * ( r0 + δ ) * L * Treac − Tperm * dξ 2
[
(
)] (
perm = Cpsweep * Fsweep + Cp Hperm * FH 2 + df H 2 * Tperm + dTperm − Tref 2
The only two species involved are the sweep gas and the hydrogen. As a reminder, the sweep gas is used to pull the hydrogen out of the tube faster. It can be made of pure hydrogen (use of a by-pass) or steam (easy to condensate). The drawback of using hydrogen is that the partial pressure of hydrogen in the permeate side is equal to the total pressure (fixed at 1 bar in our case). When our system is working at 30 bars, it won’t have an important effect on our system, but at low temperature, it will decrease our reactor’s performances. The drawback of using steam is the need of a condenser after the reactor in order to separate water from hydrogen. However, a simple cooler as a separation unit is needed to reach a temperature below one hundred Celsius degrees. Also, the steam used has to be without impurities (Nitrogen for example). In our case, we will consider the sweep gas to be steam (data easier to find).
The heat capacities involved in this heat balance have to be calculated at a specific temperature, which means that since the hydrogen flux through the membrane is at the reaction temperature, its heat capacity is calculated for the same temperature.
56
Let’s develop the exit term:
[Cp
perm sweep
(
)] (
) [
] (
perm * Fsweep + Cp Hperm * FH 2 + df H 2 * Tperm + dTperm − Tref = Cpsweep * Fsweep + Cp Hperm * FH 2 * Tperm − Tref 2 2
[
)
]
perm + Cpsweep * Fsweep + Cp Hperm * FH 2 * dTperm 2 0 + FCH * 4
0 + FCH * 4
dYH 2 dξ dYH 2 dξ
(
)
* Tperm − Tref * dξ * dTperm * dξ
As we saw in the first heat balance, the first term can be simplified with the entrance and the last term is a second order differential, so it can be neglected.
[Cp
perm sweep
]
0 * Fsweep + Cp Hperm * FH 2 * dTperm = Cp reac H 2 * FCH 4 * 2
0 − Cp reac H 2 * FCH 4 *
dYH 2
dξ dYH 2 dξ
(
)
(
)
* Treac − Tref * dξ * Tperm − Tref * dξ
(
)
+ hmemb * 2π * ( r0 + δ ) * L * Treac − Tperm * dξ
Therefore, the second heat balance can be written:
€
dYH 2 reac 0 dTperm Cp H 2 * FCH 4 * dξ * Treac − Tperm + hmemb * 2π * ( r0 + δ ) * L * Treac − Tperm = perm dξ Cpsweep * Fsweep + Cp Hperm * FH 2 2
(
)
(
)
The final differential system for the adiabatic membrane reactor is the following:
€
57
Remark: Our situation is considering the absence of energy exchange between outside and the shell. In order to obtain the non-isothermal case, we will have to add a new convection term in entrance of the heat balance (of the reaction side of course): conv dQshell = UOR * 2π * Re * L * (Tout − Treac ) * dξ
€
58
MATLAB model
The adiabatic MATLAB model is based on the isothermal model we built in the previous chapter. However, we have to be very careful concerning the terms depending on temperature. They have to be integrated in the functions since T is now a variable of the system (and not a simple parameter).
Changes from isothermal model
Even if the heart of the model will remain the same, multiple changes have to be considered. First of all, the parameters depending on temperature have to be moved from the main program to the functions. The constants, defined in the first part of the program, are considered as global variables so the parameters can be calculated in the PBR and MR functions. Secondly, since the temperature is the main characteristic of our system, its precision will have to be important. A one-degree step will be the maximum considered. However, since the integration is taking much more time than in the isothermal case, a temperature step of 0.1 degree will be the minimum chosen (and only for a small range of temperature). Finally, as we will see in the next paragraph, the resolution for the temperatures can diverge for some points. This will cause the system to give impossible values (negative or infinite). In order to obtain accurate results, an “if” statement is written after the use of the functions but before the loops end. The condition is on the sign of the temperatures and the range of the conversion (between 0 and 1). If these conditions are not fulfilled, the previous value for this parameter (at the previous Temperature step) is replacing them. Indeed, considering the precision of our system, this is a minor error and it doesn’t have incidence on the functions.
59
Functions and changes of volumetric flow rate
Multiple statements have to be added to the functions in the adiabatic case compared to the isothermal one. First of all, the following parameters have to be calculated for the temperature of the system. In the case of the membrane reactor, heat capacities for water and hydrogen have to be calculated for both temperatures, due to the difference of system considered (shell and tube). Moreover, new elements have to be defined in the functions. The derivative of flowrate in report to the dimensionless length is an element needed in the heat balances and have to be calculated at each step. Each flowrate is needed too, and are determined thanks to the partial pressures, the temperature (of the system considered) and the volumetric flowrate. This last parameter needs to be initiated and calculated at each step. Since it is not a differential variable, we can’t initiate it in the main program. Therefore, an “if” statement is used on the Temperature, since this parameter will vary all over the reactor and won’t be equal to its initial value at any length of the reactor. The initial volumetric flowrate is calculated at the beginning of the function (depending on the inlet and the initial state of the system) and the volumetric flowrate for each step is calculated after the total flowrate is determined, so at the end of the function. The condition is given just after calculating the initial volumetric flowrate. If the temperature is equal to the initial temperature, the initial volumetric flowrate is used. Otherwise, the calculated flowrate is used. While simulating the MR adiabatic model, curious divergence appeared in our results while the paced bed reactor gives us accurate graphic representations. These perturbations may be explained by the precision used (1e-4), not small enough to take in account the multiple variations of our system, or by reasons not determined yet.
60
Results PBR adiabatic
The first results to show are the results found for the Packed Bed reactor working in adiabatic. In the following results, the dotted line (up for conversion, down for temperature) represents a pressure of 2bars and then, each line has the respective pressure [5 ; 10 ; 20 ; 30 ; 40 ; 50 ; 100] bars:
Figure 34:PBR Adiabatic, CH4 conv. vs inlet Temperature
61
Figure 35:PBR adiabatic, Exit Temperature vs Inlet temperature
One of the major results appearing from these graphic representations is the Temperature reached by the outlet. The temperature is decreasing along the reactor, confirming the endothermic behavior of our process. This was the result we were expecting, proving that our heat balances are correct.
We plotted through these two graphs the effect of pressure on our system. As we saw in isothermal, an increase of pressure in the PBR provokes the decrease of methane conversion due to the system equilibrium. Moreover, we can observe something that we didn’t notice previously. Even if the methane conversion obtained for the outlet is decreasing with the increase of pressure, this evolution is reaching a limit and its decrease is not constant. If between 2 and 50 bars, the conversion is decreasing from 0.3 to 0.175, it is only decreasing to 0.15 at 100 bars. The same behavior is observed for the outlet temperature, which is increasing with the pressure but not in a linear way.
Another interesting thing we can notice is the evolution of the temperature outlet with the pressure. The reactor temperature decrease is less important for when the pressure is increased. The reason for this behavior cannot be understood directly; therefore we need more results.
62
Let’s represent the evolution of temperature along the reactor for different pressures at 850 Celsius:
Figure 36: Effect of pressure, Temperature vs ksi
The temperature used is high in order to compensate the energy loss. We can try to represent the conversion and the reactor temperature for a lower temperature, which means lower conversion. The blue line is referring to the temperature, the green line to the conversion; for a pressure of 15 bars and a temperature of 700 Celsius, and for a pressure of 30 bars and a temperature of 700 Celsius:
63
Figure 37: PBR adiabatic, CH4 conv and Temp vs ksi, Pt=15 bars, T = 700C
Figure 38: PBR adiabatic, CH4 conv and Temp vs ksi, Pt=30 bars, T = 700C
64
It appears that most of the reaction is happening in the first part of the reactor. Final methane conversion and Temperature are obtained at a dimensionless factor of 0.1 in both cases. From these results, the reason why the outlet temperature is higher for high pressure is due to the PBR behavior. Indeed, the conversion reached at high pressure is lower than at low pressure, which means fewer molecules are reacting in an endothermic behavior. The energy consumed by the reactions is therefore smaller. This conclusion is also showing how the energy consumed by the reaction is the major heat exchange in our heat balance. Increasing the pressure will decrease the heat loss by conduction since it is linked to the pressure, but this will also decrease the conversion of methane. Let’s consider the reaction temperature versus the dimensionless length, for multiple inlet temperature. The pressure is a constant fixed at 30 bars and 850 Celsius:
Figure 39: Temperature vs ksi, at different inlet temperature
As for the previous case with the different pressures, the temperature is decreasing mostly in the first part of the reactor. And the behavior of these temperatures with an increase of the temperature inlet is similar to the increase of pressure we saw previously. The temperature keeps on decreasing along the reactor at low inlet temperature while it is reaching its outlet value sooner for high inlet temperature.
65
Again, it appears that most of the heat loss is due to the reaction, which happens at the beginning of the reactor. The other heat loss, due to conduction, is decreased while working at high pressure and high temperature.
Comparison PBR isothermal vs Adiabatic
Let’s remind that an adiabatic reactor will not be built for a methane steam reforming producing, since this reaction is endothermic and needs energy to produce good results (energy brought by furnaces usually). However, in order to build the non-isothermal model, this adiabatic model is necessary in order to take in account the energy loss due to the reaction and the conduction. This is why comparing Isothermal and Adiabatic results, even if it won’t give any realistic conclusions, is interesting as a mathematical problem.
Let’s consider the methane conversion depending on the pressure and inlet temperature, for both models:
Figure 40: iso vs adia, effect of pressure in a PBR
The results are what we were expecting. The performances of the adiabatic reactor are really low compared to the isothermal reactor ones. 66
It is interesting however to notice that at low temperature (400Celsius), the results are close and around the same values, whatever the pressure (methane conversion between 0.3 and 0.1). This is showing again how our system needs an energy source.
MR Adiabatic
As we said earlier, the results for the adiabatic membrane reactor are not as knit as for the packed bed reactor, due to convergence issues. The only interesting temperature, the reactor’s one, is plotted. As we did for the packed bed reactor, the methane conversion and the outlet temperature are graphically represented against the inlet temperature, and so for multiple pressures. The red line is at 2 bars while the other lines are respectively at 5, 10, 15 and 30 bars:
Figure 41: CH4 conv. vs inlet temperature, effect of pressure for MR
67
Figure 42: Outlet Temperature vs inlet temperature, effect of pressure for MR
These results are showing multiple points. First of all, the general behavior of our adiabatic membrane reactor is the same than our isothermal membrane reactor. An increase of temperature or pressure provokes an increase of the methane conversion. Secondly, the outlet temperature depends on the conversion, in the same way we saw in the adiabatic packed bed model. A low pressure will provoke a low conversion, so a lesser need in energy and a global temperature higher than for high pressure. Finally, we can see as for the adiabatic PBR that the parameter evolution is not linear but curvy. High temperature is provoking an increase in the hydrogen production (endothermic kinetic behavior) and therefore, the system shows this typical exponential behavior, even if it is dampen compared to the isothermal case.
We saw in the isothermal simulations that high pressures were giving interesting results, therefore we will consider now simulations for the values defined previously (pressure superior to 16 bars). This is allowing us to do simulations taking less time to be completed. Also, the divergences issues were encountered for low pressures; therefore part of our code can be simplified like this. Let’s now see the how the temperature and conversion behave inside the reactor: 68
Figure 43: MR adiabatic, CH4 conv and Temp vs ksi, Pt=30 bars, T = 700C
Contrary to the PBR, the membrane reactor is showing a constant evolution of the conversion (and so, the temperature) at high pressure. This evolution has the form of an exponential increase. And it seems the maximum conversion is not obtained at ksi=1. The temperature behavior is similar to what we saw in the PBR case, with a decrease parallel to the conversion increase.
It is also interesting to represent the permeate temperature and the reactor temperature in the reactor:
69
Figure 44: MR adiabatic, permeate and reactor temperatures vs ksi, T=700C and Pt=30bars
The permeate temperature is slightly superior to the reactor temperature. From our previous results, we saw how the reactions demand in energy was the main cause of the temperature drop and not the conduction of the gases. From this plot, we can understand the energy transfer through the membrane as an important element in our energy balance. Indeed, if this transfer were minimal, we would have seen a large difference between both temperatures. Since they appear to be similar, it means the convection between both sides is extremely high, bringing energy to the system. This is an interesting result since it will be identical to the non-isothermal case. Indeed, the convection through the membrane depends only on the temperature difference between sides, so it is not linked to the external temperature.
Comparison adiabatic MR vs adiabatic PBR
As we said earlier, the adiabatic model is not an industrial practical model. It is only a mathematical representation of the methane steam reforming while considering the system intrinsic energy balances.
70
However, the difference between adiabatic and non-isothermal is only a convection term added to both models. Therefore, a comparison between the packed bed reactor and the membrane reactor will in fact represent the energy effect of the membrane on the system. This won’t be comparable to the non-isothermal case however, since the convection term to add depends on multiple variables and is not linear. But the results are still interesting
The first comparison to do is on the outlet temperature against the inlet temperature, for different pressure (2, 5, 10, 20, 30, and 50 bars). The red curves represent the PBR while the blue curves are for the MR:
Figure 45: PBR vs MR, T out vs T in, effect of pressure
It appears that the energy loss in the membrane reactor is much more important than in the PBR. From the conversions observed in the previous part, this seems normal. The conversion graph should show us a larger conversion for the MR than for the PBR. An interesting thing is the behavior of the curves. They are much more “horizontal” for the MR than the PBR. The temperature difference between the systems is therefore increasing with the temperature inlet.
71
Since we saw that temperature and conversion were linked to each other, we will now represent the conversion on a similar graph. The following plot is defined for the same parameters, however the lower pressure has been taken at 1.5 bars:
Figure 46: MR vs PBR, CH4 conv. vs T, effect of pressure
This graphic representation is extremely interesting. In the isothermal case, the conversion difference between PBR and MR was large, especially at medium temperature. Here, we still observe a high methane conversion for the MR but not as important. At 350 Celsius, this difference is 0.02 at 1.5 bars and 0.06 at 30 bars. At 650 Celsius, it becomes 0.15 at 30 bars and 0 at 1.5 bars. The conversion is even better for the PBR at 1.5 bars and at temperatures higher than 700 Celsius. The choice of a minimal pressure at 1.5 bars has been made to show the energy effect on our system. This may be due to the fact that the MR is having higher energy loss due to the tube side. However at a pressure superior to 2 bars, the MR is more efficient in terms of methane conversion than the PBR. And we saw in the isothermal case (which is our ideal case) the range of temperature and pressure we are should operate, and high pressures were preferred. And at high pressure, the divergence problems are barely appearing, as we can see in the following plot:
72
Figure 47:PBR vs MR, high pressures, conversion vs temperature
Figure 48:PBR vs MR, high pressures, out temp vs temp
73
The parameter Delta has finally been plotted to show this:
Figure 49: Delta vs T in, effect of pressure (2 and 30 bars)
Figure 50:Delta vs T, effect of pressure (15 and 30 bars)
74
As we can see on this plot, Delta is barely reaching 0.08 at 2 bars for an inlet temperature of 900 Celsius, while it is almost 0.28 at 30 bars for the same temperature. The irregularity observed at 705 Celsius is one of the divergences we talked about in the previous chapter. Compared to the values of Delta in the isothermal case, the performances of the adiabatic model are very poor. However, this is a mathematical observation and not a physical one.
As we did in the isothermal case, we can see the effect of the catalyst deactivation on our system. The parameter alpha has been introduced in the same way:
Figure 51:Effect of catalyst deactivation on Delta, Pt=30 bars
75
Figure 52:Effect of catalyst deactivation on methane conversion, Pt=30 bars
Figure 53:Effect of catalyst deactivation on Temperature, Pt=30 bars
76
We can notice that the effect of alpha is similar to what it was in isothermal. The conversion decreased by only a few percents due to the deactivation. The outlet temperature is therefore higher since lesser reactants are consumed, meaning less energy used by the system.
One drawback of our results is the impossibility to compare them to sources. Indeed, the mathematical models developed which considered the heat balances wrote results for the nonisothermal case, not the adiabatic one. The comparison between these models and ours will be pointless.
77
Non-isothermal The 1D non-isothermal steady state model has been studied intensively. It is defined from the adiabatic model by including the convection term corresponding to the energy transfer between the exterior to the reaction side through the wall of the reactor. The major difficulty is the heat transfer coefficient expression since the heat transfer between the reactor wall, catalyst and the exterior temperature bring multiple parameters to define in our model. It is usually agreed to use an overall heat transfer coefficient in order to take in account all of these parameters.
Differential system The differential systems are extremely close to the adiabatic ones. We just have to add the following convection term:
For the membrane reactor, the temperature involved is the reaction side one.
Therefore, the non-isothermal packed bed reactor differential equation set is:
And the non-isothermal membrane reactor differential equation set is:
78
U overall
Different equations can be found for this overall heat transfer coefficient and among all works, Dixon (1996) published results that have been recognized and used in papers. The expression he defined is the following:
The parameters used in this expression are defined as:
79
Let’s first define the wall heat transfer coefficient. The expression we use has been reported by Tsotsas and Schlünder (1990) in the line of Dixon’s work:
The parameter dt is the equivalent diameter of a cylindrical tube with the same section of the annular tube. In our case, the ratio is equal to 11, allowing us to use the expression. As a reminder, the dimensionless number Rep is the Reynolds particular, and Pr is the Prandtl number, defined by: ; The expression of the mixture viscosity has been already defined when we calculated the pressure drop in the reactor, the mixture density is the sum of the densities at the correspondent ratio (flux of the specie over the total flux), Dp is the particle diameter and Gx is mass specific gas flow rate, defined with the densities, the area (linked to the void fraction) and the volumetric flow rate.
Finally, the mixture thermal conductivity is defined by Wassiljewa (1904) expression:
with
80
The effective radial thermal bed conductivity is defined by the following expression:
With
As we can see, these parameters depend on the mixture composition and the reaction temperature. Therefore, they have to be calculated at each step of our model. Even for one inlet temperature and one pressure, in conditions that were not diverging in the adiabatic case, no result was obtained. Better precisions have been tried but the same problem appeared, and the time of calculation was extremely important. In order to obtain results, simplifications in the code should be done. Limit values should be defined for the parameters provoking the divergence, in order to avoid obtaining imaginary parts in the final results.
U overall constant A simplification that can be made is to consider the heat transfer coefficient as a constant. Considering the expression we have just defined, this simplification is putting aside a lot of different parameters and therefore, the results will be far from the actual behavior of the reactor. However by considering the flux as a constant, we will obtain a representation of our system closer to the physical non-isothermal reactor than what the adiabatic model was giving us. 81
The value defined by Madia et al. (1999) will be used as the overall heat transfer coefficient in a packed bed reactor. Associated to this, we will use their reactor wall thickness: ⎧UOR = 817200 J.m −2 .h −1 .K −1 ⎨ ⎩ l = 2.54 cm
In a first time, we will consider an outside temperature of 500 Celsius, and we will take the same temperature for the inlet. € Temperature along the reactor, Tw=T0=500Celsius
Figure 54:Non-‐isothermal MR and PBR reactors, Temperature vs ksi, Pt=30bars
From these results and our previous observations, we can understand the temperature profile: the initial decrease is due to the major consumption of reactants at the beginning of the
82
reactor (until ksi = 0.1), then the heat due to convection is compensating it and provokes the increase of temperature, until the heat equilibrium is obtained (Treac=Twall=500C=773K). The packed bed reactor is showing an equilibrium at the exit while the membrane reactor is having a 45 degrees difference.
Figure 55:Non-‐isothermal MR and PBR reactors, CH4 conversion vs ksi, Pt=30bars
The conversion results are showing a much higher conversion for the membrane reactor than for the packed bed reactor, as we saw for the isothermal case. We can notice that the maximum conversion obtained for the MR is only 0.62 while it was 0.9 in isothermal, while the PFR conversion is pretty much the same than the one obtained in isothermal (for the same temperature). This is explaining the temperature profile. The PBR is reaching is equilibrium concentration and therefore, the temperature is reaching its stable value, equal to the exterior temperature. The MR is experiencing a continuous methane conversion and didn’t reach its maximum by the end of the reactor, the temperature profiles is the direct result of this. For a lower pressure (16 bars), we found the following results: 83
Figure 56:Non-‐isothermal MR and PBR, Temperature vs ksi, Pt=16bars
Figure 57:Non-‐isothermal MR and PBR, CH4 conversion vs ksi, Pt=16bars
84
The results for a pressure of 16bars are similar to what we found at 30 bars. Again, the conversion obtained with the membrane reactor is lower than what we found in the isothermal case (0.58 instead of 0.75) while it is unchanged for the packed bed reactor case. And again, the PBR reached an almost maximum methane conversion at the exit of the reactor whereas for the membrane reactor, the conversion is still in evolution. It is interesting to notice that in the high-pressure case, the same conversion is obtained with a difference 50 degrees whereas for 16 bars, this difference is only 35 degrees.
Temperature along the reactor, Tw=T0=600Celsius
Figure 58:Non-‐Isothermal MR and PBR, CH4 conv vs ksi
85
Figure 59:Non-‐Isothermal MR and PBR, Temperature vs ksi
The increase of the inlet and exterior temperatures has provoked the increase of the exit temperature and conversion for both reactors, whatever the pressure. Similar to the previous case, the PBR is showing system stability at the reactor exit whereas the membrane reactor is showing a constant increase of conversion without reaching a maximum. At this temperature, the high-pressure membrane reactor is reaching its stable state at the very exit of the reactor. Indeed, the methane conversion is approaching 1 (0.988 from the numerical results) and the temperature is increasing more than for the other pressures. We can see that the temperature at high pressure is even passing the temperatures for the other two pressures. The same behavior can be observed for the 20 bars situation, since the exit temperature is equal to the 16bars case while the conversion is higher (0.955 > 0.932).
Another thing to look at is the difference between reactions and permeate temperatures in the membrane reactor case:
86
Figure 60:Non-‐Isothermal MR, Reaction and Permeate Temperature vs ksi, Tw=T0=600C, Pt=30bars
We can notice how close the two temperatures are to each other. The permeate temperature is higher at the beginning of the reactor and lower at the end. This is confirming our hypothesis saying that at the reactor beginning, the heat transfer is more the most part the reaction. This is consuming the system energy and the temperature drops abruptly, until the temperature difference is so important that the convection term becomes major. Therefore, the reaction temperature increases and the reaction can be continued.
This can be seen as followed: while the membrane is countering the thermodynamic limitations, a constant gain of energy by convection is countering the kinetic limitations.
87
Conclusion
We can conclude this thesis by summarizing the work that has been done. First of all, an isothermal mathematical model has been created and showed extremely accurate results. These results have been compared with multiple sources, allowing us to validate our model. The main conclusion was the highest performances obtained for the Membrane reactor in a specific set of parameters. It has been reminded that the isothermal model is the less accurate model possible, since the energy balances are not considered. Therefore, the adiabatic model has been written and presented interesting results in the same set of parameters. The influence of the reactions energy demand was indeed revealed to be the major heat transfer in our system. We also saw how the convection through the membrane was extremely important. This model however is not representative of reality since the external convection is not considered. This is why we tried to develop the non-isothermal model. However, the complexity of the heat transfer coefficient made the system diverge and no results were obtained. A constant value of this parameter has therefore been chosen and results show an atypical behavior of the system. High methane conversion was reached in the membrane reactor case for a set of parameters similar to the isothermal case. The validity of this reactor over the traditional packed bed reactor has been confirmed. Nevertheless, further work has to be conducted in order to diminish the simplifications made, such as the overall heat transfer value, the number of our model dimensions or the steady state.
88
References
[1] Lee F. Brown ; International Journal of Hydrogen Energy, 26 (2001) 381-397 [2] Savvas Vasileiadis, Zoe Ziaka-Vasileiadou ; Separation and Purification Technology, 34 (2004) 213-225 [3] Julien Godat, Francois Marechal ; Journal of Power Sources 118 (2003) 411-423 [4] Jianguo Xu, Gilbert F. Froment ; American Institute of Chemical Engineering Journal, 35 (1989) 88-103 [5] Engin Ayturk, Nikolas K. Kazantzis, Yi Hua Ma ; Energy & Environmental Science, 2 (2009) 430-438 [6] Jun Shu, Bernard P.A. Grandjean, Serge Kaliaguine ; Applied Catalysis A: General, 119 (1994) 305-325 [7] Giuseppe Barbieri, Fransesco P. Di Maio ; Industrial Engineering Chemical Research, 36 (1997) 2121-2127 [8] Fausto Galluci, Luca Paturzo, Angelo Basile ; International Journal of Hydrogen Energy, 29 (2004) 611-617 [9] M. De Falco, L. Di Paola, L. Marrelli ; International Journal of Hydrogen Energy, 32 (2007) 2902-2913 [10] Byron Smith R J, Muruganandam Loganathan, Murthy Shekhar Shantha ; International Journal of Chemical Reactor Engineering, 8 (2010) Review R4 [11] G. Madia, G. Barbieri, E. Drioli ; Canadian Journal of Chemical Engineering, 77 (1999) 698-706 [12] A. Bottino, A. Comite, G. Capannelli, R. Di Felice, P. Pinacci ; Catalysis Today, 118 (2006) 214-222 [13] Fabiano A.N. Fernandes, Aldo B. Soares Jr ; Fuel, 85 (2006) 569-573 [14] Anthony Dixon ; Chemical Engineering and Processing, 35 (1996) 323-331 [15] D.L. Hoang, S.H. Chan ; Applied Catalysis A: General, 268 (2004) 207-216 [16] J.M. Smith, H.C. Van Ness, M.M. Abbott ; 7th Edition (2005) [17] Estéban Saatgjian ; 3rd Edition (1998) [18] H. Scott Fogler ; 3rd Edition (1999) [19] Bruce E. Poling, John M. Prausnitz, John P. O’Connell ; 5th Edition (2007) 89
Appendices MATLAB Code – Isothermal Main program % % % % % % % %
(1) CH4 + H2O = (2) CO + H2O = (3) CH4 + 2*H2O = (1-a) CH4 (2-b) H2O (3-c) CO (4-d) CO2 (5-e) H2
CO + CO2 + CO2 +
3*H2 H2 4*H2
global theta sweep Pp Ka Ke kr C1 C2 Pt Eta kr=zeros(3,1);krr=zeros(3,1);E=zeros(3,1);Kar=zeros(5,1);Ka=zeros(5,1); DH=zeros(5,1);Ke=zeros(3,1);G0=zeros(5,1); Tr=ones(5,1);theta=zeros(5,1); % CONSTANTES % inchangeables divers r0=1*2.54e-2; % Radius of the tube without membrane, [m] Ri=2.5*2.54e-2; % Inner radius of the reactor, [m] PhiB=0.5; % Bed porosity rhocat=2100; % catalyst density [kg/m3], 10%Ni/Al2O3, should be 4900 but Engin said 2100 Eta(1)=0.01; % Oveall effectiveness factor, r(1) Eta(2)=0.3; % Oveall effectiveness factor, r(2) Eta(3)=0.01; % Oveall effectiveness factor, r(3) Pp=1; % Permeate pressure [bar] Pi=3.14159; R=8.314; % [J/K.mol] % inchangeables cinetiques krr(1)=1.8e-4; % [kmol.bar0.5/kgcat.h] krr(2)=7.6; % [kmol/kgcat.h.bar] krr(3)=2.2e-5; % [kmol.bar0.5/kgcat.h] E(1)=240100; % [J/mol] E(2)=67100; % [J/mol] E(3)=243900; % [J/mol] Kar(3)=40.9; % [1/bar] Kar(5)=2.9e-2; % [1/bar] Kar(1)=1.8e-1; % [1/bar], Tref=823K Kar(2)=0.4; % Tref=823K DH(3)=-70700; % [J/mol] DH(5)=-82900; % [J/mol] DH(1)=-38300; % [J/mol] DH(2)=88700; % [J/mol]
90
% Gibbs Energy coefficients Gc=[-75.3 7.6e-2 1.9e-5 % CH4 -241.7 4.2e-2 7.4e-6 % H2O -109.9 -9.2e-2 1.5e-6 % CO -393.4 3.8e-3 1.3e-6]; % CO2 % Rk: G(H2)=0 because pure specie % stochiometric % r1 r2 r3 mu=[-1 0 -1 -1 -1 -2 1 -1 0 0 1 1 3 1 4];
coefficients % % % % %
CH4 H2O CO CO2 H2
% Reference Temperature for kinetics parameters, different from Tref0=298K Tr=648*Tr; for i=1:2 Tr(i)=823; % [K] end % reactor & membrane L=7; delta=5e-6; rit=r0+delta; % inlet Ptot=[2;5;10;20;30]; kmax=length(Ptot); %Pt=30; Fa0=1; theta(2)=3; theta(3)=0; theta(4)=0; theta(5)=0.01; sweep=10; T0=350+273.15;
% Reactor Length, [m] % membrane thickness [meter] % inner tube radius [m]
% % % % % % % %
Total pressure [bar] CH4 initial feed, [kmol/h] F0(H2O)/F0(CH4) F0(CO)/F0(CH4) F0(CO2)/F0(CH4) F0(H2)/F0(CH4) sweep factor F0(sweep)/F0(CH4) Inlet temperature [K]
nmax=401; X1=zeros(nmax,kmax);X2=zeros(nmax,kmax); X3=zeros(nmax,kmax);X4=zeros(nmax,kmax); X5=zeros(nmax,kmax); for k=1:kmax Pt=Ptot(k); %Pt=Pt0+(k-1)*1; %L=Lsize(k); %theta(2)=H2Oratio+(k-1)*1; for n=1:nmax %delta=delta0+(n-1)*1e-6; %rit=r0+delta; %sweep=s0+(n-1); %T=T0;
% % % %
If If If If
% % % %
Pt is the parameter/Matrix Pt is an axis L is the parameter/Matrix m is the parameter/Matrix
If delta is the parameter/Matrix inner tube radius [m] If sweep is the parameter/Matrix If T is a constant
91
T=T0+(n-1)*1; % kinetic values for i=1:3 kr(i)=krr(i)*exp(E(i)*(1/648-1/T)/R); % same units than kr end for i=1:5 Ka(i)=Kar(i)*exp(DH(i)*(1/Tr(i)-1/T)/R); % same units than Kr end for i=1:4 G0(i)=Gc(i,1)+Gc(i,2)*T+Gc(i,3)*T^2; end DGrxn=zeros(3,1); % is changing for each T for i=1:3 for j=1:4 DGrxn(i)=DGrxn(i)+mu(j,i)*G0(j); end Ke(i)=exp(-DGrxn(i)/(R*T*1e-3)); end % other parameters Ac=Pi*(Ri^2-rit^2); Qpd=6.3227*(1e-3)*exp(-15630/(R*T)); [m3.um.h.atm0.5/m2] rhoB=rhocat*(1-PhiB); C1=rhoB*Ac*L/Fa0; and MR C2=Qpd*2*Pi*rit*L/(Fa0*delta*22.4);
% area, [m2] % H2 permeability in Pd, % Bed density [kg/m3] % Constant for diff. eq. 1 and 2, PBR % Constant for diff. eq. 3, MR
% INTEGRATION u0=[0 0]; % initial values ksispan=(0:0.0001:1); % precision [~,S1]=ode23s('f2',ksispan,u0); % PBR reactor X1(n,k)=real(S1(end,1)); % methane finale conversion X2(n,k)=real(S1(end,2)); % carbon dioxide finale conversion u1(:,k)=S1(:,1); % methane conversion along ksi u0=[0 0 0]; [ksi,S2]=ode23s('f3',ksispan,u0); X3(n,k)=real(S2(end,1)); X4(n,k)=real(S2(end,2)); X5(n,k)=real(S2(end,3)); u1(:,k)=S2(:,1);
% % % %
% MR reactor methane finale conversion carbon dioxide finale conversion hydrogene finale conversion methane conversion along ksi
end end T0=T0-273.15;T=T-273.15; % Temperature expressed in Celsius Tspan=(T0:1:T); Pspan=(Pt0:1:Pt); %Dspan=(delta0*1e6:1:delta*1e6); %Sspan=(s0:1:sweep); D(:,:)=X3(:,:)-X1(:,:); % Delta = XCH4(MR) - XCH4(PBR)
92
clf %figure(100) %mesh(Pspan,Tspan,D) %title('Effect of pressure and temperature on Methane conversion - L=7m, d=5um, m=3, Fa0=1kmol/h') figure(1) plot(Tspan,X1(:,1),'--',Tspan,X1(:,2),'--',Tspan,X1(:,3),'-',Tspan,X1(:,4),'--',Tspan,X1(:,5),'--') hold on,plot(Tspan,X3(:,1),Tspan,X3(:,2),Tspan,X3(:,3),Tspan,X2(:,4),Tspan,X3(:,5) ) legend('PBR 2','PBR 5','PBR 10','PBR 20','PBR 30','MR 2','MR 5','MR 10','MR 20','MR 30') xlabel('Temperature (Celsius)') ylabel('Methane Conversion') axis([T0 T 0 1]) title('Membrane VS Packed Bed - Conversion of Methane (L=7m, d=5um, m=3, Fa0=1kmol/h)') figure(2) plot(Tspan,X2(:,1),'--',Tspan,X2(:,2),'--',Tspan,X2(:,3),'-',Tspan,X2(:,4),'--',Tspan,X2(:,5),'--') hold on,plot(Tspan,X4(:,1),Tspan,X4(:,2),Tspan,X4(:,3),Tspan,X4(:,4),Tspan,X4(:,5) ) legend('PBR 2','PBR 5','PBR 10','PBR 20','PBR 30','MR 2','MR 5','MR 10','MR 20','MR 30') xlabel('Temperature (Celsius)') ylabel('CO2 Conversion') axis([T0 T 0 1]) title('Membrane VS Packed Bed - Conversion of carbon Dioxide (L=7m, d=5um, m=3, Fa0=1kmol/h)') figure(3) plot(Tspan,X5(:,1),Tspan,X5(:,2),Tspan,X5(:,3),Tspan,X5(:,4),Tspan,X5(:,5)) legend('MR 2','MR 5','MR 10','MR 20','MR 30') xlabel('Temperature (Celsius)') ylabel('Hydrogen production') axis([T0 T 0 4]) title('Membrane VS Packed Bed - Hydrogen production (L=7m, d=5um, m=3, Fa0=1kmol/h)') figure(4) plot(ksi,u1(:,1),'--',ksi,u1(:,2),'--',ksi,u1(:,3),'--',ksi,u1(:,4),'-',ksi,u1(:,5),'--',ksi,u1(:,6),'--') hold on,plot(ksi,u2(:,1),ksi,u2(:,2),ksi,u2(:,3),ksi,u2(:,4),ksi,u2(:,5),ksi,u2(:, 6)) legend('PBR 2','PBR 5','PBR 10','PBR 20','PBR 30','MR 2','MR 5','MR 10','MR 20','MR 30') xlabel('ksi') ylabel('conversion CH4')
93
title('effect of reactor length - T=750 deg. C') figure(5) plot(Tspan,D(:,1),Tspan,D(:,2),Tspan,D(:,3),Tspan,D(:,4),Tspan,D(:,5)) legend('2 bars','5 bars','10 bars','20 bars','30 bars') xlabel('Temperature (Celsius)') ylabel('Delta = XCH4(MR) - XCH4(PBR)') axis([T0 T 0 1]) title('Membrane VS Packed Bed - Delta representation (m=3, L=7m, d=5um, Fa0=1kmol/h)')
MATLAB Code – Isothermal Functions % Packed bed reactor function, Isothermal % the variable Y=u(3) is not considered function diff=f2(ksi,u) global theta Ka Ke kr C1 Pt Eta % Expression of the partial pressures, functions of methane and carbon % dioxide conversions, theta and Pt. S=Pt/(1+theta(2)+theta(3)+theta(4)+theta(5)+2*u(1)); P(1)=(1-u(1))*S; P(2)=(theta(2)-u(1)-u(2))*S; P(3)=(theta(3)+u(1)-u(2))*S; P(4)=(theta(4)+u(2))*S; P(5)=(theta(5)+3*u(1)+u(2))*S; % Expression of the reaction rates, functions of Kinetic coefficients, % partial pressures and conversions. DEN2=(1+Ka(1)*P(1)+Ka(3)*P(3)+Ka(5)*P(5)+Ka(2)*P(2)/P(5))^2; r(1)=Eta(1)*(kr(1)/(P(5)^2.5))*(P(1)*P(2)-(P(5)^3)*P(3)/Ke(1))/DEN2; r(2)=Eta(2)*(kr(2)/P(5))*(P(3)*P(2)-P(5)*P(4)/Ke(2))/DEN2; r(3)=Eta(3)*(kr(3)/(P(5)^3.5))*(P(1)*P(2)^2-(P(5)^4)*P(4)/Ke(3))/DEN2; % Differential System diff1=C1*(r(1)+r(3)); diff2=C1*(r(2)+r(3)); diff=[diff1 diff2]';
% Membrane reactor function, Isothermal function diff=f3(ksi,u) global theta Ka Ke kr C1 C2 Pt Pp sweep Eta % Expression of the partial pressures, functions of methane and carbon % dioxide conversions, hydrogen permeability, theta and Pt. S=Pt/(1+theta(2)+theta(3)+theta(4)+theta(5)+2*u(1)-u(3)); P(1)=(1-u(1))*S; P(2)=(theta(2)-u(1)-u(2))*S; P(3)=(theta(3)+u(1)-u(2))*S; P(4)=(theta(4)+u(2))*S; P(5)=(theta(5)+3*u(1)+u(2)-u(3))*S;
94
% Expression of the reaction rates, functions of Kinetic coefficients, % partial pressures and conversions. DEN2=(1+Ka(1)*P(1)+Ka(3)*P(3)+Ka(5)*P(5)+Ka(2)*P(2)/P(5))^2; r(1)=Eta(1)*(kr(1)/(P(5)^2.5))*(P(1)*P(2)-(P(5)^3)*P(3)/Ke(1))/DEN2; r(2)=Eta(2)*(kr(2)/P(5))*(P(3)*P(2)-P(5)*P(4)/Ke(2))/DEN2; r(3)=Eta(3)*(kr(3)/(P(5)^3.5))*(P(1)*P(2)^2-(P(5)^4)*P(4)/Ke(3))/DEN2; % Pressure of hydrogen in the permeate side Pmem=Pp*u(3)/(u(3)+sweep); % Differential system diff1=C1*(r(1)+r(3)); diff2=C1*(r(2)+r(3)); diff3=C2*(P(5)^0.5-Pmem^0.5); diff=[diff1 diff2 diff3]';
MATLAB Code – Adiabatic Main Code global theta Kar krr E DH Gc C1 Pt Pp sweep R mu Pi rit L Fa0 delta DeltarH0 MC Tref0 Tref Ac rhoB Eta T0 kr=zeros(3,1);krr=zeros(3,1);E=zeros(3,1);Kar=zeros(5,1);Ka=zeros(5,1); DH=zeros(5,1);Ke=zeros(3,1);G0=zeros(5,1); Tref=ones(5,1);theta=zeros(5,1);DeltarH0=zeros(3,1);Eta=zeros(3,1); % CONSTANTES % inchangeables divers r0=1*2.54e-2; % Outer radius of the membrane, [m] Ri=2.5*2.54e-2;% Inner radius of the reactor, [m] PhiB=0.5; % Bed porosity rhocat=2100; % catalyst density [kg/m3], 10%Ni/Al2O3, should be 4900 but Engin said 2100 Eta(1)=0.01; % Oveall effectiveness factor for r(1) Eta(2)=0.3; % Oveall effectiveness factor for r(2) Eta(3)=0.01; % Oveall effectiveness factor for r(3) Pp=1; % Permeate pressure [bar] Pi=3.14159; R=8.314; % [J/K.mol] Tref0=298.15; % Ref temperature for Delta r H calculations [K] % inchangeables cinetiques krr(1)=1.8e-4; % [kmol.bar0.5/kgcat.h] krr(2)=7.6; % [kmol/kgcat.h.bar] krr(3)=2.2e-5; % [kmol.bar0.5/kgcat.h] E(1)=240100; % [J/mol] E(2)=67100; % [J/mol] E(3)=243900; % [J/mol] Kar(3)=40.9; % [1/bar] Kar(5)=2.9e-2; % [1/bar] Kar(1)=1.8e-1; % [1/bar], Tref=823K Kar(2)=0.4; % Tref=823K DH(3)=-70700; % [J/mol] DH(5)=-82900; % [J/mol] DH(1)=-38300; % [J/mol] DH(2)=88700; % [J/mol]
95
% Reactions enthalpies at 298K DeltarH0(1)=206310; DeltarH0(3)=165110; DeltarH0(2)=DeltarH0(3)-DeltarH0(1); % Gibbs Energy coefficients Gc=[-75.3 7.6e-2 1.9e-5 -241.7 4.2e-2 7.4e-6 -109.9 -9.2e-2 1.5e-6 -393.4 3.8e-3 1.3e-6]; % Rk: G(H2)=0 because pure specie % CP coefficients % A B*e3 C*e6 D*e-5 MC=[1.702 9.081 -2.164 0 % CH4 3.470 1.450 0 0.121 % H2O 3.376 0.557 0 -0.031 % CO 5.457 1.045 0 -1.157 % CO2 3.249 0.422 0 0.083]; % H2 % stochiometric coefficients % 1 2 3 mu=[-1 0 -1 % CH4 -1 -1 -2 % H2O 1 -1 0 % CO 0 1 1 % CO2 3 1 4]; % H2 % Reference Temperature for kinetics Tref=648*Tref; for i=1:2 Tref(i)=823; % [K] end
% [J/mol] % [J/mol] % [J/mol]
parameters, different from Tref0=298K
% reactor & membrane L=7; % Reactor Length, [m] delta=5e-6; % membrane thickness [meter] rit=r0+delta; % inner tube radius [m] % inlet Ptot=[2;5;10;15;30]; % Total pressure [bar] %Ptot=[30]; kmax=length(Ptot); %Pt=30; % [bar] Fa0=1; % CH4 initial feed, [kmol/h] theta(2)=3; % F0(H2O)/F0(CH4) theta(3)=0; % F0(CO)/F0(CH4) theta(4)=0; % F0(CO2)/F0(CH4) theta(5)=0.001; % F0(H2)/F0(CH4) sweep=20; % sweep factor F0(sweep)/F0(CH4) Tini=400+273.15; % Inlet temperature [K] nmax=401; X=zeros(5,nmax,kmax,4); DT=zeros(nmax,kmax); for k=1:kmax Pt=Ptot(k); for n=1:nmax T0=Tini+(n-1);
96
Ac=Pi*(Ri^2-rit^2); % area, [m2] rhoB=rhocat*(1-PhiB); % Bed density [kg/m3] C1=rhoB*Ac*L/Fa0; %________________________________________________________________________ % % INTEGRATION %________________________________________________________________________ % Packed Bed reactor adia u0=[0 0 T0]; ksispan=(0:0.0001:1); [ksia,u1]=ode15s('fadia1',ksispan,u0); X1(n,k)=u1(end,1); X3(n,k)=u1(end,3)-273.15; SolC(:,n,k)=u1(:,1); SolT(:,n,k)=u1(:,3)-273.15;
% % % %
methane conversion Outlet Temperature System conversion System temperature
[-] [K] [-] [K]
DT(n,k)=T0-(u1(end,3));
% Inlet - outlet [degrees]
% Packed Bed reactor iso u0=[0 0 T0]; ksispan=(0:0.001:1); [ksib,u2]=ode15s('f2',ksispan,u0); X2(n,k)=u2(end,1);
% methane conversion [-]
% Due to convergence problem, a retro-analysis of the results suppress the % imaginary parts and give real results. for u=1:4 for n=3:nmax-1 for k=1:kmax for i=1:5 if imag(X(i,n,k,u))==0 X(i,n,k,u)=X(i,n,k,u); else X(i,n,k,u)=2*X(i,n-1,k,u)-X(i,n-2,k,u); end end end end end end end
MATLAB Code – Adiabatic Functions % Packed bed reactor function, Adiabatic % the variable Y=u(3) is not considered, Temperature is the third variable function diff=fadia1(ksi,u) global theta Kar krr E DH Gc C1 Pt R mu DeltarH0 MC Tref0 Tref Ac rhoB Fa0 L Eta Ftot T0
97
kr=zeros(3,1);Ka=zeros(5,1);Ke=zeros(3,1); G0=zeros(5,1);Cp=zeros(5,1);F=zeros(5,1);dFdksi=zeros(5,1); K1=0;K2=0;Q=0; % third variable defined as Temperature T=u(3); % Expression of volumetric flowrate Q Q0=Fa0*(1+theta(2)+theta(3)+theta(4)+theta(5))*R*T/Pt; % Initial Q if T==T0 Q=Q0; % Initial situation else Q=Ftot*R*T/Pt*1e-3; % Other cases, "1e-3" for units end % Expression of heat capacities % - general expression (for heat balance) for j=1:5 Cp(j)=R*(MC(j,1)+MC(j,2)*T*1e-3+MC(j,3)*T^2*1e-6+MC(j,4)*(1/T^2)*1e5); %[J/(mol.K)] end % - integrated values (for Enthalpy expressions) DCp=zeros(3,1); for i=1:3 for j=1:5 DCp(i)=DCp(i)+mu(j,i)*R*(MC(j,1)*(T-Tref0)+MC(j,2)/2*(T^2Tref0^2)*1e-3+MC(j,3)/3*(T^3-Tref0^3)*1e-6-MC(j,4)*(1/T-1/Tref0)*1e5); %[J/(mol.K)] end end % Expression of Enthalpies DeltarH=zeros(3,1); DeltarH(1)=DeltarH0(1)+DCp(1); % [J/mol] DeltarH(3)=DeltarH0(3)+DCp(3); % [J/mol] DeltarH(2)=DeltarH(3)-DeltarH(1); % [J/mol]
% Expression of Kinetic parameters % - reaction rate constants for i=1:3 kr(i)=krr(i)*exp(E(i)*(1/648-1/T)/R); % same units than krr end % - Adsorption constants for i=1:5 Ka(i)=Kar(i)*exp(DH(i)*(1/Tref(i)-1/T)/R); % same units than Kar end % - Gibbs energie constants for i=1:4 G0(i)=Gc(i,1)+Gc(i,2)*T+Gc(i,3)*T^2; end DGrxn=zeros(3,1); for i=1:3 for j=1:4
98
DGrxn(i)=DGrxn(i)+mu(j,i)*G0(j); end % Equilibrium constants Ke(i)=exp(-DGrxn(i)/(R*T*1e-3)); end
% [-] no unit
% Expression of the partial pressures, functions of methane and carbon % dioxide conversions, theta and Pt. S=Pt/(1+theta(2)+theta(3)+theta(4)+theta(5)+2*u(1)); P(1)=(1-u(1))*S; % [bar] P(2)=(theta(2)-u(1)-u(2))*S; % [bar] P(3)=(theta(3)+u(1)-u(2))*S; % [bar] P(4)=(theta(4)+u(2))*S; % [bar] P(5)=(theta(5)+3*u(1)+u(2))*S; % [bar] % Expression of the reaction rates, functions of Kinetic coefficients, % partial pressures and conversions. DEN2=(1+Ka(1)*P(1)+Ka(3)*P(3)+Ka(5)*P(5)+Ka(2)*P(2)/P(5))^2; r(1)=Eta(1)*(kr(1)/(P(5)^2.5))*(P(1)*P(2)-(P(5)^3)*P(3)/Ke(1))/DEN2; % [kmol/(kg.h)] r(2)=Eta(2)*(kr(2)/P(5))*(P(3)*P(2)-P(5)*P(4)/Ke(2))/DEN2; % [kmol/(kg.h)] r(3)=Eta(3)*(kr(3)/(P(5)^3.5))*(P(1)*P(2)^2-(P(5)^4)*P(4)/Ke(3))/DEN2; % [kmol/(kg.h)] % Expression of differential flows dF/dksi, from mass balances dFdksi(1)=-rhoB*Ac*L*(r(1)+r(3))*1e3; %[mol/h] dFdksi(2)=-rhoB*Ac*L*(r(1)+r(2)+2*r(3))*1e3; %[mol/h] dFdksi(3)=rhoB*Ac*L*(r(1)-r(2))*1e3; %[mol/h] dFdksi(4)=rhoB*Ac*L*(r(2)+r(3))*1e3; %[mol/h] dFdksi(5)=rhoB*Ac*L*(3*r(1)+r(2)+4*r(3))*1e3; %[mol/h] % Expression of general flowrates, function of temperature, pressure and Q Ftot=0; for j=1:5 F(j)=P(j)*Q/(R*T)*1e3; % [mol/h] Ftot=Ftot+F(j); % [mol/h] end % Rk: needs of "1e3" in order to have same units in heat balance % Expression of the heat balance % - denominator for j=1:5 K2=K2+F(j)*Cp(j); %[J/(K.h)] end % - numerator for i=1:3 K1=K1+(-DeltarH(i))*r(i)*rhoB*Ac*L*1e3; for j=1:5 K1=K1-Cp(j)*dFdksi(j)*(T-Tref0); %[J/h] end end % Differential system
99
diff1=C1*(r(1)+r(3)); % methane conversion [-] diff2=C1*(r(2)+r(3)); % Carbon dioxide conversion [-] diff3=K1/K2; % Temperature evolution [K] diff=[diff1 diff2 diff3]';
% membrane reactor function, Adiabatic function diff=fadia2(ksi,u) global theta Kar krr E DH Gc C1 Pt R mu DeltarH0 MC T0 Tref0 Tref Ac rhoB Fa0 L Eta delta Pi rit Pp sweep Umemb Ftot kr=zeros(3,1);Ka=zeros(5,1);Ke=zeros(3,1);G0=zeros(5,1); Cpr=zeros(5,1);Cpp=zeros(5,1);DCpr=zeros(3,1);F=zeros(5,1);dFdksi=zeros(5,1); K1=0;K2=0;Q=0; % Fourth and Fifthe variables are the temperatures Tr=abs(real(u(4))); % Reaction Temperature, index "r" Tp=abs(real(u(5))); % Permeate Temperature, index "p" % Expression of volumetric flowrate Q, [J/(h.bar)] Q0=Fa0*1e3*(1+theta(2)+theta(3)+theta(4)+theta(5))*R*T/Pt; % Initial Q if T==T0 Q=Q0; % Initial situation else Q=Ftot*R*T/Pt; % Other cases end % Expression of heat capacities % - general expression (for heat balance) for j=1:5 Cpr(j)=R*(MC(j,1)+MC(j,2)*Tr*1e-3+MC(j,3)*Tr^2*1e6+MC(j,4)*(1/Tr^2)*1e5); % Reactor side [J/(mol.K)] Cpp(j)=R*(MC(j,1)+MC(j,2)*Tp*1e-3+MC(j,3)*Tp^2*1e6+MC(j,4)*(1/Tp^2)*1e5); % Permeate side [J/(mol.K)] end % - integrated values (for Enthalpy expressions), Reactor side for i=1:3 for j=1:5 DCpr(i)=DCpr(i)+mu(j,i)*R*(MC(j,1)*(Tr-Tref0)+MC(j,2)/2*(Tr^2Tref0^2)*1e-3+MC(j,3)/3*(Tr^3-Tref0^3)*1e-6-MC(j,4)*(1/Tr-1/Tref0)*1e5); %[J/(mol.K)] end end % Expression of Enthalpies DeltarH=zeros(3,1); DeltarH(1)=DeltarH0(1)+DCpr(1); % [J/mol] DeltarH(3)=DeltarH0(3)+DCpr(3); % [J/mol] DeltarH(2)=DeltarH(3)-DeltarH(1); % [J/mol] % Expression of Kinetic parameters for i=1:3 kr(i)=krr(i)*exp(E(i)*(1/648-1/Tr)/R); % same units than krr end for i=1:5 Ka(i)=Kar(i)*exp(DH(i)*(1/Tref(i)-1/Tr)/R); % same units than Kar end for i=1:4
100
G0(i)=Gc(i,1)+Gc(i,2)*Tr+Gc(i,3)*Tr^2; end DGrxn=zeros(3,1); for i=1:3 for j=1:4 DGrxn(i)=DGrxn(i)+mu(j,i)*G0(j); end Ke(i)=exp(-DGrxn(i)/(R*Tr*1e-3)); end % Definition of Partial Pressures S=Pt/(1+theta(2)+theta(3)+theta(4)+theta(5)+2*u(1)-u(3)); P(1)=(1-u(1))*S; % [bar] P(2)=(theta(2)-u(1)-u(2))*S; % [bar] P(3)=(theta(3)+u(1)-u(2))*S; % [bar] P(4)=(theta(4)+u(2))*S; % [bar] P(5)=(theta(5)+3*u(1)+u(2)-u(3))*S; % [bar] % Condition on Hrydrogen fraction (it can't be inferior to 0) if P(5)
