EngOpt 2008 - International Conference on Engineering Optimization Rio de Janeiro, Brazil, 01 - 05 June 2008.

Optimization of a Membrane Reactor for Hydrogen Production Through Methane Steam Reforming Using Experimental Design Techniques and NPSOL L. C. Silva; V. V. Murata; C. E. Hori; A. J. Assis School of Chemical Engineering – Federal University of Uberlandia – Uberlandia – MG – Brazil, Phone: (+55)34 3239 4189 - E-mail: [email protected]

1. Abstract In this work two mathematical models for methane steam reforming in membrane reactors were developed. The first one is a steady state, non isothermal, non isobaric and one dimensional model derived from material and energy balances and validated using experimental data from the literature. It is referred as full model. The influence of inlet reactor pressure (PR0), methane feed flow rate (FCH40), sweep gas flow rate (FI), external reactor temperature (TW) and steam to methane feed flow ratio (M) on methane conversion (XCH4 ) and hydrogen recovery (YH2 ) were evaluated through simulations. The second model, referred as simplified model, was obtained though the response surface technique considering a central composite planning and results from 41 simulations using the full model. This simplified model was included into a constrained optimization problem whose objective function was defined as the summation of both XCH4 and YH2 . The optimized codified values related with PR0, FI and TW reached the maximum value while the FCH4 and M reached the minimum value. The resulting XCH4 and YH2 were 99.94% and 98.19% respectively. These preliminay optimized values, which are higher than any other condition, are consistent with the parametric analysis accomplished with the full model. 2. Keywords: Methane steam reforming, membrane reactors, hydrogen, NPSOL. 3. Introduction Hydrogen is used as fuel and as a feedstock in important processes such as ammonia and methanol production and Fischer-Tropsch synthesis, besides the increasing use in fuel cells. Methane steam reforming is the predominant industrial method used for hydrogen production [5; 7] and involves three reactions (Eq.(1) - (3)), according to Xu and Froment [10], that are highly endothermic and occur with increasing number of moles. Thus, these reactions are favored at high temperatures and low pressures. In the conventional process these reactions are carried out in long tubular reactors filled with nickel based catalysts, that operates at temperatures above 850◦ C, pressures between 3 and 25 atm and steam to methane feed flow ratio between 2 and 4. The maximum conversion achieved is the equilibrium conversion and it is around 70 to 80% [4]. On the other hand, the development of Pd-based (palladium-based) membrane reactors is changing this setting. Pd-based membranes are considered to be 100% permselective to hydrogen, the membrane continuously removes hydrogen produced in the catalytic reaction zone thus pushing the chemical equilibrium and allowing higher methane conversions at lower temperature. CH4 + H2 O CO + 3 H2 CO + H2 O CO2 + H2

− ∆H298K = −206.1 kJ/mol

(1)

− ∆H298K = 41.15 kJ/mol

(2)

− ∆H298K = −165.0 kJ/mol

CH4 + 2 H2 O CO2 + 4 H2

(3)

Several theoretical and experimental works showed methane conversions above 90% in membrane reactors (e.g. [3; 2; 9]), but just few investigated hydrogen recovery. Nevertheless, no authors have reported both high conversions and high recoveries at the same time. Thus, the main objectives of this work are: (i) to carry out the mathematical modeling of a Pd-based membrane reactor, for methane steam reforming; (ii) to validate the model with experimental data from the literature; (iii) to find the best operating conditions that optimize hydrogen recovery and methane conversion. 4. Mathematical modeling A schematic diagram of the Pd-based membrane reactor considered in this work is illustrated in Fig. 1. The inner tube is the Pd-based membrane and the outer one is the shell. The catalyst is packed in the annular region. The inert or sweep gas pass into the inner tube in co-current flow mode. 1

Figure 1: Membrane reactor scheme considered in this work. The mathematical model considered real gases for thermodynamic properties calculations. Additional hypotheses adopted in the model development were: • Isothermal/isobaric and nonisothermal/nonisobaric cases; • Pseudo-homogeneous model; • Steady-state operation; • Plug-flow behavior was adopted for the permeate and reaction zone; • Inexistence of boundary layer on the membrane; • Hydrogen diffusion in the membrane is the rate-determining step for hydrogen permeation; The balance equations are: • Mass balance in the reaction chamber: For i = CH4 , CO, CO2 , H2 O

3 X

dFi =W νij rj dZ j=1 For i = H2

3 X

(4)

dFi =W νij rj − JH2 Am dZ j=1

(5)

where JH2 is given by: JH2 =

E A0 e( R Tm ) 0.5 (pH2 ) − pH2 p )0.5 δ

(6)

• Mass balance in the permeate side: dFH2p dZ

= JH2 Am

(7)

• Energy balance in the reaction chamber: 3 X 1 dTR ∂P dvr + (−∆Hj ) rj −HH2 JH2 Am Q1 − Q2 + W ρr Cvr vr · = TR dZ ∂T V dZ Van j=1 • Energy balance in the permeate side:

(8)

ρp Cvp vp ·

dTp ∂P dvp 1 + H J A = Tp + Q 2 H H m 2 2 dZ ∂T V dZ Vpe

(9)

Q1 is the heat exchanged between the reactor wall and the reaction chamber, and Q2 is the heat exchanged between the reaction chamber and the permeate: Q1 = U1 A1 (Tw − TR )

(10)

Q2 = U2 A2 (TR − TP )

(11)

the overall heat transfer coefficients were taken from Ohmori et al. [5]. • The pressure drop in the catalytic bed was calculated using: dPR ρr u2s = −f dZ dp

(12)

where f is calculated by the Ergun’s equation: 1−ε f= ε3

b(1 − ε) a+ Re

(13)

The reactor dimensions used in this work were taken from Shu et al. [7]. All thermodynamic mixture properties were calculated using the Virial equation of state, with van der Waals mixing rules. The gases and the gaseous mixture viscosities were calculated by Chung’s method [6]. The values of constants a and b of Ergun’s equation are 150 and 1.75 respectively [1]. Other catalyst properties used were taken from Xu and Froment [10], including kinetic expressions and parameters. The following dimensionless variables were used to make the model dimensionless: Pp =

PP PR TR Z TP FCH FH2 O ; Pr = 0 ; T r = ; z = ; Tp = ; fCH4 = 0 4 ; fH2 O = 0 0 PR PR Tw L Tw FCH4 M FCH 4 fH2 =

FH2p FH2 FCO ; fH2p = ; fCO2 = 0 2 0 0 4FCH4 4FCH4 FCH4

The resultant model equations are a set of nine nonlinear coupled ordinary differential equations integrated using DLSODE (Double precision Livermore Solver for Ordinary Differential Equations) package of FORTRAN subroutines. This solver makes use of Adams methods (predictor-corrector) in the nonstiff case, and Backward Differentiation Formula (BDF) methods (the Gear methods) in the stiff case. In this work it was used the BDF method with an initial step of 1 × 10−14 , a relative tolerance of 1 × 10−15 and an absolute tolerance of 1 × 10−10 . The following initial conditions (z=0) were used: fCH4 (0) = Tp (0) = Pr (0) = Tr (0) = fH2 O (0) = 1 fCO (0) = fCO2 (0) = fH2p (0) = 0 fH2 (0) = 1 × 10−6 Methane conversion and hydrogen recovery were defined as: XCH4 =

0 FCH − FCH4 0 4 = 1 − fCH 0 4 FCH4

YH2 =

fH2p

fH2p 0 + fH2 − fH 2

(14) (15)

5. Results and discussion 5.1. Model validation Model validation was carried out comparing the simulated data with the experimental data reported by Shu et al. [7], at the same operational conditions and reactor design parameters. Shu et al. [7] considered two types of membranes in their work: a palladium supported in stainless steel (Pd/SS) and an alloy of palladium and silver supported on stainless steel (Pd-Ag/SS). Additionally, they considered only isothermal and isobaric conditions. Thus, in order to make the comparison in the same basis, the operating temperature and pressure from Shu et al. [7] were considered to be the external reactor temperature and the inlet pressure, respectively. The agreement between the simulated results and the experimental ones was good as it can be seen in the Fig. 2. Fig. 2 (a) shows conversion against the molar steam to methane feed flow ratio, and the best agreement between experimental and simulated data was found at steam to methane feed flow ratio around 3. A good agreement is showed in Fig. 2 (b), since the model tracked the methane conversion behavior in all range of temperature. A standard or reference operational condition for membrane reactors in methane steam reforming is, according to Shu et al. [7]: PR0 = 136000 Pa, Pp = 101325 Pa, TW = 773.15 K, M = 3, FI = 2.75×10−5 mol/s and FCH40 = 2.75×10−5 mol/s. Fig.3 shows a typical composition profile using this standard operational condition, where it can be seen that most of the reaction takes place up to 10% reactor’s length.

(a)

(b)

Figure 2: Model validation. (a) Methane conversion behavior versus steam to methane feed flow ratio and (b) Methane conversion behavior versus temperature.

5.2. Temperature and pressure effects In order to find out the influence of temperature both in methane conversion and hydrogen recovery, the external reactor temperature (TW) was changed from 300 to 700◦ C (or from 573.15 to 973.15 K), in nine steps. The results are plotted in Fig. 4(a). High conversions and recoveries are achieved at high temperatures. According to Sjardin et al. [8], the temperature range commonly used to operate Pd-based membrane reactors is up to 600◦ C (or 873.15 K). Moreover, at temperatures around 300◦ C, and in the presence of hydrogen, palladium has a transition between two different lattice phases. This can cause membrane stress, resulting in microcracks and embrittlement of the membrane. Using the model developed in this work, temperature of 600◦ C and standard conditions for other operating conditions, methane conversion and hydrogen recovery calculated are around 80% and 40%, respectively. The inlet reaction pressure (PR0) effect on methane conversion was experimentally studied by Tonge et al. [9]. These authors adopted a range of 101.325 to 506.625 kPa, which is relatively bigger than what is reported in other studies (e.g. [7; 2]). In this work, a range of 100 to 500 kPa was adopted and the results are in Fig. 4(b). Contrasting with conventional reactors, the increase of reaction pressure has a positive effect on methane conversion and hydrogen recovery. With an inlet pressure of 506.625 kPa and standard condition for other variables, a methane conversion of 65.23% and a hydrogen recovery of 83.24% were achieved.

Figure 3: Composition profiles along the reactor for standard operational condition

(a)

(b)

Figure 4: Effect of the inlet pressure and external temperature on the methane conversion and hydrogen recovery 5.3. The effects of steam to methane feed flow ratio and inert flow rate To study the effect of steam to methane feed flow ratio (M) it was considered a range of 2 to 6, which is intermediate between the ranges used by Shu et al. [7] and Gallucci et al. [2]. Additionally, this range is wide enough to allow reasonable conclusion. The reader should attempt to the fact that M values lower than two implies in carbon formation and M values higher than six reduce excessively the hydrogen partial pressure by dilution. As it can be seen in Fig. 5 (a), an increase in M implies in a raise in methane conversion, and reduction in hydrogen recovery. Therefore, an optimization of both methane conversion and hydrogen recovery from the parametric standpoint is difficult to achieve. Fig. 5 (b) shows the effect of the inert (or sweep gas) molar flow rate (FI). When the inert or sweep gas flow rate is increased, both methane conversion and hydrogen recovery increase, reaching the same value, 74%, in a sweep gas flow rate of 1.4×10−4 mols/s.

(a)

(b)

Figure 5: Effect of steam to methane feed and sweep gas flow rate on methane conversion and hydrogen recovery 5.4. Effect of methane feed flow rate Increasing methane feed flow rate (FCH40), decreases both methane conversion and hydrogen recovery. This effect is more proeminent in hydrogen recovery at lower values of methane feed and exhibit an exponencial behavior, as it can be seen in Fig. 6.

Figure 6: Effect of methane feed flow rate on the methane conversion and hydrogen recovery 6. Influence of Model Parameters by Experimental Design As shown in previous results (Fig. 4, 5, 6) obtained from the full model, several variables strongly affect both methane conversion and hydrogen recovery. From an process standpoint view, the objective of a reforming process is to maximize both methane conversion and hydrogen recovery, subject to bounds in variables and model restrictions. Thus, a rigorous optimization for this problem is a very difficult task. Considering the complex nature of the full model a preliminary optimization step which uses a simplified model was done. Due to the nonlinear nature of the parameters behavior, a central composite experimental design was chosen to fit a second order surface response on methane conversion and hydrogen recovery. Five factors were analyzed, namely: inlet reactor pressure, methane feed flow rate, sweep gas flow rate, external reactor temperature and steam to methane feed flow ratio. In this case no more than one duplicate value in the center is necessary, because there is no experimental error. An α value of 1.596 was chosen to the central composite design be orthogonal. This experimental design led to 43 simulations in different conditions, in the following five

parameters level: Table 1: Central composite variables level Var\ level PR0(X1 ) FCH40 (X2 ) FI(X3 ) TW (X4 ) M (X5 )

-1.596 101.325 1.000×10−5 2.75×10−5 573.150 2.500

-1 177.001 3.38×10−5 4.804×10−5 629.165 3.155

0 303.975 7.375×10−5 8.250×10−5 723.150 4.250

1 430.949 1.137×10−4 1.1696×10−4 817.135 5.346

1.596 506.625 1.3750×10−4 1.3750×10−4 873.150 6.000

The following codification equations were used: F CH40 − 7.375 × 10−5 F I − 8.250 × 10−5 ; X = 3 3.994×−5 3.446 × 10−5 T W − 723.150 M − 4.250 X4 = ; X5 = 93.985 1.096 The following fitted surfaces were found as a function of the codified variables: X1 =

P R0 − 303.975 ; 202.65

X2 =

XCH4 =43.8235 + 2.5037 · X1 − 7.00623 · X2 + 5.88484 · X22 + 3.70743 · X3 + 30.8484 · X4 + + 3.0722 · X5 + 2.553 · X1 · X4 YH2 =52.59577 + 14.41726 · X1 − 3.06291 · X12 − 11.9214 · X2 + 4.01804 · X22 + 8.81369 · X3 + − 4.5154 · X5 + 2.605086 · X1 · X4

(16)

(17)

The first surface response is the methane conversion; some effects were eliminated based on the confidence interval of 95%. The model fitted has a correlation coefficient of the R2 = 0.973. The second one is the surface response for the hydrogen recovery. This equation has a correlation coefficient of R2 = 0.942. As it can be seen on equations above, the mainly quadratic effects are methane feed flow ratio (X2 ) and the inlet pressure (X1 ), and only the linear interaction between the inlet pressure and the steam to methane feed flow rate was relevant. 7. Optimization using NPSOL and surface responses NPSOL is a set of FORTRAN 77 subroutines for minimizing a smooth function subject to constraints, which may include simple bounds on the variables, linear constraints and smooth nonlinear constraints, and may also be used for unconstrained, bound-constrained and linearly constrained optimization. NPSOL uses a sequential quadratic programming (SQP) algorithm, in which each search direction is the solution of a QP sub problem. Bounds, linear constraints and nonlinear constraints are treated separately [11]. In this work NPSOL was used to optimize methane conversion and hydrogen recovery using surface response derived in previous item. Due to the fact that methane conversion and hydrogen recovery have the same order of magnitude, the sum of these variables was defined as the objective function. Therefore the optimization problem was posed as: Maximize: XCH4 + YH2 Subject to: XCH4 =43.8235 + 2.5037 · X1 − 7.00623 · X2 + 5.88484 · X22 + 3.70743 · X3 + 30.8484 · X4 + + 3.0722 · X5 + 2.553 · X1 · X4 YH2 =52.59577 + 14.41726 · X1 − 3.06291 · X12 − 11.9214 · X2 + 4.01804 · X22 + 8.81369 · X3 + − 4.5154 · X5 + 2.605086 · X1 · X4 −1.596 ≤ X1 ≤ 1.596 −1.596 ≤ X2 ≤ 1.596

(18)

(19)

−1.596 ≤ X3 ≤ 1.596 −1.596 ≤ X4 ≤ 1.596 −1.596 ≤ X5 ≤ 1.596 The selection of lower and upper bounds of the decision variables are based on the operating limits of the membrane reactors and Shu’s work [7]. In fact, they are the α value in the central composite design, but they were previously chosen with this purpose. The optimization led to the following values: X2 = X5 = −1.596

X1 = X3 = X4 = 1.596; Which codified back led to: P R0 = 506.625 kP a;

F CH40 = 1.0 × 10−5 mol/s; T W = 873.15 K;

F I = 1.375 × 10−5 mol/s

M = 2.5

The optimized methane conversion and hydrogen recovery calculated through the complete model using the above optimal values were 99.94% and 98.19%, which are higher than any other condition. 8. Conclusions The membrane reactor mathematical model developed showed a good agreement with the experimental data. The influence of five important parameters on methane conversion and hydrogen recovery was analyzed. Steam to methane feed flow ratio has an opposite effect when the objective is to maximize both methane conversion and hydrogen recovery. Optimization procedure was done in two steps. First, experimental design techniques were employed to adjust a surface response that shows how the target variables depend on optimizing parameters. Second, the optimum was found using the optimizer code NPSOL. It was shown that high methane conversion and hydrogen recovery, 99.94 and 98.19%, respectively, can be achieved using the following parameters values: PR0 = 506.625 kPa; FCH40 = 1.0×10−5 mol/s; FI = 1.375× 10−5 mol/s; TW = 873.15 K; M = 2.5. These preliminay optimized values, which are higher than any other condition, are consistent with the parametric analysis accomplished with the full model. 9. List of symbols a - Constant in Ergun’s equation [−] A0 - Pre-exponential factor of the permeability

h

N LT

T 2L M

n i

2 A1 - Heat transfer area between the reaction chamber and the external reactor temperature L 2 A2 - Heat transfer area between the reaction chamber and the permeate side L Am - Permeation area of membrane L2 b - Constant in Ergun’s equation [−] h i 2

Cvp - Specific heat at constant volume of the permeate side mixture TM2 NL Θ h i 2 Cvr - Specific heat at constant volume of the reaction side mixture TM2 NL Θ dp - Equivalent diameter particle [L] h of catalyst i E - Activation energy

M L2 T 2N

f - Friction factor [−] fH2 p - Dimensionless flow rate of hydrogen in the permeate side [−] fi - Dimensionless flow rate of i component in the reaction side [−] F CH40 - Methane feed flow rate N T FH2 p - Hydrogen flow rate in the permeate N T Fi - Molar flow rate of component i in the reaction side N T F I - Molar flow rate of inert gas [mol] h or sweep i

HH2 - Hydrogen enthalpy

M L2 T 2N

JH2 - Hydrogen permeation flux

N L2 T

L - Reactor length [L] M - Steam to methane feed flow H2 O e CH4 [−] ratio P R0 - Inlet reactor pressure TM 2L PR - Reactor pressure TM 2L Pr - Dimensionless reaction pressure [−] Pp - Permeate pressure TM 2L pH2 - Hydrogen partial pressure in the reaction side TM 2L pH2p - Hydrogen partial pressure in the permeate side TM 2L

h i 2 Q1 - Heat exchanged between the reaction side and external environment MT L2 i h 2 Q2 - Heat exchanged between the reaction side and the permeate side MT L2 i h Kmol rj - Rate of the reaction j Kg h cat h2 i R - Universal gas constant TM2 NL Θ

Re - Reynolds’ number [−] Tm - Membrane temperature [Θ] Tp - Dimensionless permeate temperature [−] Tr - Dimensionless temperature of the reaction side [−] TR - Temperature of the reaction side [Θ] Tw - External temperature [Θ] us - Gaseous velocity of the reaction chamber TL U1 - Overall heat transfer coefficient between reaction side and external side TM M3 Θ U2 - Overall heat transfer coefficient between reaction and permeation zone T 3 Θ vp - Mean velocity of the gaseous mixture in the permeate TL vr - Mean velocity of the mixture in the reaction chamber TL 3gaseous Van - Annular volume L Vpe - Permeate volume L3 XCH4 - Methane conversion [−] W - Catalyst mass [M ] YH2 - Hydrogen recovery [−] z - Dimensionless axial position [−] Z - Axial position [−] Greek letters δ - Membrane thickness [L] h i ∆Hj - Heat of the reaction j

M L2 T 2N

ε - Void fraction [−] νij - Stoichiometric coefficient of the component i in the reaction j [−] ρp - Permeate density LM3 ρr - Reaction density LM3 Subscripts i = CH4 , CO, H2 , H2 O, CO2 . j = 1, 2, 3 (reactions 1,2,3) k = CH4 , CO, H2 , H2 O. p = permeate side r = reaction side Acknowledgements This work has been supported by Brazilian funding agencies, CAPES and CNPq (Grant no 475934/2006-7). The authors are grateful for the assistance of M.Sc. Adriene Artiaga Pfeifer and Eng. Alice Medeiros de Lima in the use of DIRCOL and NPSOL codes.

10. References [1] Falco, M. D. T.; Di Paola, L.; Marrelli and Nardella, L. Simulation of large-scale membrane reformers by a two-dimensional model, Chemical Engineering Journal, 2006, 128, 115-125. [1] Froment, K. and Bischoff, K. B. Chemical Reactor Analysis and Design. John Wiley, New York, 1990. [2] Gallucci, F.; Paturzo, L.; Fama, A.; Basile, A. Experimental Study of the Methane Steam Reforming Reaction in a Dense Pd/Ag Membrane Reactor, Industrial and Engineering Chemistry Research , 2004, 43, 928-933. [3] Lin, Y.-M.; Liu, S.-L.; Chuanga, C.-H. and Chub, Y-T.. Effect of incipient removal of hydrogen through palladium membrane on the converion of methane steam reforming: experimental and modeling, Catalysis Today, 2003, 82, 127-139. [4] Ogden, J. M. Review of small stationary reformers for hydrogen production, Report No. IEA/H2/TR-02/002. Princeton, USA: Princeton University. 2001. [5] Ohmori, W. Y. T.; YU, W.; Yamamoto, A.; Endo, A.; Nakaiwa, M.; Hayakawa, T. and Itoh, N. Simulation of a porous ceramic membrane reactor for hydrogen production, International Journal of Hydrogen Energy, 2005, 30(10), 1071-1079, . [6] Poling, B. E.; Prausnits. M. J.; O’Connell, J. P. The Properties of Gases and Liquids, 5th ed., McGraw-Hill, 2004. [7] Shu, J.; Grandjean, P. A. and Kaliaguine, S. Methane steam reforming in asymmetric Pd- and Pd-Ag/porous SS membrane reactors, Applied Catalysis A: General, 1994, 119, 305-325. [8] Sjardin, M.; Damen, K. J. and Faaij, A. P. C. Techno-economic prospects of small-scale membrane reactors in a future hydrogen-fuelled transportation sector, Energy, 2006, 31, 2523-2555. [9] Tong, J.; Matsumura, Y.; Suda, H. and Haraya, K. Experimental Study of Steam Reforming of Methane in a Thin (6 µM) Pd-Based Membrane Reactor, Industrial and Engineering Chemistry Research, 2005, 44, 1454-1465. [10] Xu, J. and Froment, G. F. Methane Steam Reforming, Methanation and Water-Gas Shift: I. Intrinsic Kinetics, AIChE Journal , 1989, 35(1), 88-96. [11] Gill, P. E.; Murray, W.; Saunders, M. A.; Wright, M. H. User’s Guide for NPSOL (Version 4.0): A Fortran package for nonlinear programming, SOL 86-2, 1986.

Optimization of a Membrane Reactor for Hydrogen Production Through Methane Steam Reforming Using Experimental Design Techniques and NPSOL L. C. Silva; V. V. Murata; C. E. Hori; A. J. Assis School of Chemical Engineering – Federal University of Uberlandia – Uberlandia – MG – Brazil, Phone: (+55)34 3239 4189 - E-mail: [email protected]

1. Abstract In this work two mathematical models for methane steam reforming in membrane reactors were developed. The first one is a steady state, non isothermal, non isobaric and one dimensional model derived from material and energy balances and validated using experimental data from the literature. It is referred as full model. The influence of inlet reactor pressure (PR0), methane feed flow rate (FCH40), sweep gas flow rate (FI), external reactor temperature (TW) and steam to methane feed flow ratio (M) on methane conversion (XCH4 ) and hydrogen recovery (YH2 ) were evaluated through simulations. The second model, referred as simplified model, was obtained though the response surface technique considering a central composite planning and results from 41 simulations using the full model. This simplified model was included into a constrained optimization problem whose objective function was defined as the summation of both XCH4 and YH2 . The optimized codified values related with PR0, FI and TW reached the maximum value while the FCH4 and M reached the minimum value. The resulting XCH4 and YH2 were 99.94% and 98.19% respectively. These preliminay optimized values, which are higher than any other condition, are consistent with the parametric analysis accomplished with the full model. 2. Keywords: Methane steam reforming, membrane reactors, hydrogen, NPSOL. 3. Introduction Hydrogen is used as fuel and as a feedstock in important processes such as ammonia and methanol production and Fischer-Tropsch synthesis, besides the increasing use in fuel cells. Methane steam reforming is the predominant industrial method used for hydrogen production [5; 7] and involves three reactions (Eq.(1) - (3)), according to Xu and Froment [10], that are highly endothermic and occur with increasing number of moles. Thus, these reactions are favored at high temperatures and low pressures. In the conventional process these reactions are carried out in long tubular reactors filled with nickel based catalysts, that operates at temperatures above 850◦ C, pressures between 3 and 25 atm and steam to methane feed flow ratio between 2 and 4. The maximum conversion achieved is the equilibrium conversion and it is around 70 to 80% [4]. On the other hand, the development of Pd-based (palladium-based) membrane reactors is changing this setting. Pd-based membranes are considered to be 100% permselective to hydrogen, the membrane continuously removes hydrogen produced in the catalytic reaction zone thus pushing the chemical equilibrium and allowing higher methane conversions at lower temperature. CH4 + H2 O CO + 3 H2 CO + H2 O CO2 + H2

− ∆H298K = −206.1 kJ/mol

(1)

− ∆H298K = 41.15 kJ/mol

(2)

− ∆H298K = −165.0 kJ/mol

CH4 + 2 H2 O CO2 + 4 H2

(3)

Several theoretical and experimental works showed methane conversions above 90% in membrane reactors (e.g. [3; 2; 9]), but just few investigated hydrogen recovery. Nevertheless, no authors have reported both high conversions and high recoveries at the same time. Thus, the main objectives of this work are: (i) to carry out the mathematical modeling of a Pd-based membrane reactor, for methane steam reforming; (ii) to validate the model with experimental data from the literature; (iii) to find the best operating conditions that optimize hydrogen recovery and methane conversion. 4. Mathematical modeling A schematic diagram of the Pd-based membrane reactor considered in this work is illustrated in Fig. 1. The inner tube is the Pd-based membrane and the outer one is the shell. The catalyst is packed in the annular region. The inert or sweep gas pass into the inner tube in co-current flow mode. 1

Figure 1: Membrane reactor scheme considered in this work. The mathematical model considered real gases for thermodynamic properties calculations. Additional hypotheses adopted in the model development were: • Isothermal/isobaric and nonisothermal/nonisobaric cases; • Pseudo-homogeneous model; • Steady-state operation; • Plug-flow behavior was adopted for the permeate and reaction zone; • Inexistence of boundary layer on the membrane; • Hydrogen diffusion in the membrane is the rate-determining step for hydrogen permeation; The balance equations are: • Mass balance in the reaction chamber: For i = CH4 , CO, CO2 , H2 O

3 X

dFi =W νij rj dZ j=1 For i = H2

3 X

(4)

dFi =W νij rj − JH2 Am dZ j=1

(5)

where JH2 is given by: JH2 =

E A0 e( R Tm ) 0.5 (pH2 ) − pH2 p )0.5 δ

(6)

• Mass balance in the permeate side: dFH2p dZ

= JH2 Am

(7)

• Energy balance in the reaction chamber: 3 X 1 dTR ∂P dvr + (−∆Hj ) rj −HH2 JH2 Am Q1 − Q2 + W ρr Cvr vr · = TR dZ ∂T V dZ Van j=1 • Energy balance in the permeate side:

(8)

ρp Cvp vp ·

dTp ∂P dvp 1 + H J A = Tp + Q 2 H H m 2 2 dZ ∂T V dZ Vpe

(9)

Q1 is the heat exchanged between the reactor wall and the reaction chamber, and Q2 is the heat exchanged between the reaction chamber and the permeate: Q1 = U1 A1 (Tw − TR )

(10)

Q2 = U2 A2 (TR − TP )

(11)

the overall heat transfer coefficients were taken from Ohmori et al. [5]. • The pressure drop in the catalytic bed was calculated using: dPR ρr u2s = −f dZ dp

(12)

where f is calculated by the Ergun’s equation: 1−ε f= ε3

b(1 − ε) a+ Re

(13)

The reactor dimensions used in this work were taken from Shu et al. [7]. All thermodynamic mixture properties were calculated using the Virial equation of state, with van der Waals mixing rules. The gases and the gaseous mixture viscosities were calculated by Chung’s method [6]. The values of constants a and b of Ergun’s equation are 150 and 1.75 respectively [1]. Other catalyst properties used were taken from Xu and Froment [10], including kinetic expressions and parameters. The following dimensionless variables were used to make the model dimensionless: Pp =

PP PR TR Z TP FCH FH2 O ; Pr = 0 ; T r = ; z = ; Tp = ; fCH4 = 0 4 ; fH2 O = 0 0 PR PR Tw L Tw FCH4 M FCH 4 fH2 =

FH2p FH2 FCO ; fH2p = ; fCO2 = 0 2 0 0 4FCH4 4FCH4 FCH4

The resultant model equations are a set of nine nonlinear coupled ordinary differential equations integrated using DLSODE (Double precision Livermore Solver for Ordinary Differential Equations) package of FORTRAN subroutines. This solver makes use of Adams methods (predictor-corrector) in the nonstiff case, and Backward Differentiation Formula (BDF) methods (the Gear methods) in the stiff case. In this work it was used the BDF method with an initial step of 1 × 10−14 , a relative tolerance of 1 × 10−15 and an absolute tolerance of 1 × 10−10 . The following initial conditions (z=0) were used: fCH4 (0) = Tp (0) = Pr (0) = Tr (0) = fH2 O (0) = 1 fCO (0) = fCO2 (0) = fH2p (0) = 0 fH2 (0) = 1 × 10−6 Methane conversion and hydrogen recovery were defined as: XCH4 =

0 FCH − FCH4 0 4 = 1 − fCH 0 4 FCH4

YH2 =

fH2p

fH2p 0 + fH2 − fH 2

(14) (15)

5. Results and discussion 5.1. Model validation Model validation was carried out comparing the simulated data with the experimental data reported by Shu et al. [7], at the same operational conditions and reactor design parameters. Shu et al. [7] considered two types of membranes in their work: a palladium supported in stainless steel (Pd/SS) and an alloy of palladium and silver supported on stainless steel (Pd-Ag/SS). Additionally, they considered only isothermal and isobaric conditions. Thus, in order to make the comparison in the same basis, the operating temperature and pressure from Shu et al. [7] were considered to be the external reactor temperature and the inlet pressure, respectively. The agreement between the simulated results and the experimental ones was good as it can be seen in the Fig. 2. Fig. 2 (a) shows conversion against the molar steam to methane feed flow ratio, and the best agreement between experimental and simulated data was found at steam to methane feed flow ratio around 3. A good agreement is showed in Fig. 2 (b), since the model tracked the methane conversion behavior in all range of temperature. A standard or reference operational condition for membrane reactors in methane steam reforming is, according to Shu et al. [7]: PR0 = 136000 Pa, Pp = 101325 Pa, TW = 773.15 K, M = 3, FI = 2.75×10−5 mol/s and FCH40 = 2.75×10−5 mol/s. Fig.3 shows a typical composition profile using this standard operational condition, where it can be seen that most of the reaction takes place up to 10% reactor’s length.

(a)

(b)

Figure 2: Model validation. (a) Methane conversion behavior versus steam to methane feed flow ratio and (b) Methane conversion behavior versus temperature.

5.2. Temperature and pressure effects In order to find out the influence of temperature both in methane conversion and hydrogen recovery, the external reactor temperature (TW) was changed from 300 to 700◦ C (or from 573.15 to 973.15 K), in nine steps. The results are plotted in Fig. 4(a). High conversions and recoveries are achieved at high temperatures. According to Sjardin et al. [8], the temperature range commonly used to operate Pd-based membrane reactors is up to 600◦ C (or 873.15 K). Moreover, at temperatures around 300◦ C, and in the presence of hydrogen, palladium has a transition between two different lattice phases. This can cause membrane stress, resulting in microcracks and embrittlement of the membrane. Using the model developed in this work, temperature of 600◦ C and standard conditions for other operating conditions, methane conversion and hydrogen recovery calculated are around 80% and 40%, respectively. The inlet reaction pressure (PR0) effect on methane conversion was experimentally studied by Tonge et al. [9]. These authors adopted a range of 101.325 to 506.625 kPa, which is relatively bigger than what is reported in other studies (e.g. [7; 2]). In this work, a range of 100 to 500 kPa was adopted and the results are in Fig. 4(b). Contrasting with conventional reactors, the increase of reaction pressure has a positive effect on methane conversion and hydrogen recovery. With an inlet pressure of 506.625 kPa and standard condition for other variables, a methane conversion of 65.23% and a hydrogen recovery of 83.24% were achieved.

Figure 3: Composition profiles along the reactor for standard operational condition

(a)

(b)

Figure 4: Effect of the inlet pressure and external temperature on the methane conversion and hydrogen recovery 5.3. The effects of steam to methane feed flow ratio and inert flow rate To study the effect of steam to methane feed flow ratio (M) it was considered a range of 2 to 6, which is intermediate between the ranges used by Shu et al. [7] and Gallucci et al. [2]. Additionally, this range is wide enough to allow reasonable conclusion. The reader should attempt to the fact that M values lower than two implies in carbon formation and M values higher than six reduce excessively the hydrogen partial pressure by dilution. As it can be seen in Fig. 5 (a), an increase in M implies in a raise in methane conversion, and reduction in hydrogen recovery. Therefore, an optimization of both methane conversion and hydrogen recovery from the parametric standpoint is difficult to achieve. Fig. 5 (b) shows the effect of the inert (or sweep gas) molar flow rate (FI). When the inert or sweep gas flow rate is increased, both methane conversion and hydrogen recovery increase, reaching the same value, 74%, in a sweep gas flow rate of 1.4×10−4 mols/s.

(a)

(b)

Figure 5: Effect of steam to methane feed and sweep gas flow rate on methane conversion and hydrogen recovery 5.4. Effect of methane feed flow rate Increasing methane feed flow rate (FCH40), decreases both methane conversion and hydrogen recovery. This effect is more proeminent in hydrogen recovery at lower values of methane feed and exhibit an exponencial behavior, as it can be seen in Fig. 6.

Figure 6: Effect of methane feed flow rate on the methane conversion and hydrogen recovery 6. Influence of Model Parameters by Experimental Design As shown in previous results (Fig. 4, 5, 6) obtained from the full model, several variables strongly affect both methane conversion and hydrogen recovery. From an process standpoint view, the objective of a reforming process is to maximize both methane conversion and hydrogen recovery, subject to bounds in variables and model restrictions. Thus, a rigorous optimization for this problem is a very difficult task. Considering the complex nature of the full model a preliminary optimization step which uses a simplified model was done. Due to the nonlinear nature of the parameters behavior, a central composite experimental design was chosen to fit a second order surface response on methane conversion and hydrogen recovery. Five factors were analyzed, namely: inlet reactor pressure, methane feed flow rate, sweep gas flow rate, external reactor temperature and steam to methane feed flow ratio. In this case no more than one duplicate value in the center is necessary, because there is no experimental error. An α value of 1.596 was chosen to the central composite design be orthogonal. This experimental design led to 43 simulations in different conditions, in the following five

parameters level: Table 1: Central composite variables level Var\ level PR0(X1 ) FCH40 (X2 ) FI(X3 ) TW (X4 ) M (X5 )

-1.596 101.325 1.000×10−5 2.75×10−5 573.150 2.500

-1 177.001 3.38×10−5 4.804×10−5 629.165 3.155

0 303.975 7.375×10−5 8.250×10−5 723.150 4.250

1 430.949 1.137×10−4 1.1696×10−4 817.135 5.346

1.596 506.625 1.3750×10−4 1.3750×10−4 873.150 6.000

The following codification equations were used: F CH40 − 7.375 × 10−5 F I − 8.250 × 10−5 ; X = 3 3.994×−5 3.446 × 10−5 T W − 723.150 M − 4.250 X4 = ; X5 = 93.985 1.096 The following fitted surfaces were found as a function of the codified variables: X1 =

P R0 − 303.975 ; 202.65

X2 =

XCH4 =43.8235 + 2.5037 · X1 − 7.00623 · X2 + 5.88484 · X22 + 3.70743 · X3 + 30.8484 · X4 + + 3.0722 · X5 + 2.553 · X1 · X4 YH2 =52.59577 + 14.41726 · X1 − 3.06291 · X12 − 11.9214 · X2 + 4.01804 · X22 + 8.81369 · X3 + − 4.5154 · X5 + 2.605086 · X1 · X4

(16)

(17)

The first surface response is the methane conversion; some effects were eliminated based on the confidence interval of 95%. The model fitted has a correlation coefficient of the R2 = 0.973. The second one is the surface response for the hydrogen recovery. This equation has a correlation coefficient of R2 = 0.942. As it can be seen on equations above, the mainly quadratic effects are methane feed flow ratio (X2 ) and the inlet pressure (X1 ), and only the linear interaction between the inlet pressure and the steam to methane feed flow rate was relevant. 7. Optimization using NPSOL and surface responses NPSOL is a set of FORTRAN 77 subroutines for minimizing a smooth function subject to constraints, which may include simple bounds on the variables, linear constraints and smooth nonlinear constraints, and may also be used for unconstrained, bound-constrained and linearly constrained optimization. NPSOL uses a sequential quadratic programming (SQP) algorithm, in which each search direction is the solution of a QP sub problem. Bounds, linear constraints and nonlinear constraints are treated separately [11]. In this work NPSOL was used to optimize methane conversion and hydrogen recovery using surface response derived in previous item. Due to the fact that methane conversion and hydrogen recovery have the same order of magnitude, the sum of these variables was defined as the objective function. Therefore the optimization problem was posed as: Maximize: XCH4 + YH2 Subject to: XCH4 =43.8235 + 2.5037 · X1 − 7.00623 · X2 + 5.88484 · X22 + 3.70743 · X3 + 30.8484 · X4 + + 3.0722 · X5 + 2.553 · X1 · X4 YH2 =52.59577 + 14.41726 · X1 − 3.06291 · X12 − 11.9214 · X2 + 4.01804 · X22 + 8.81369 · X3 + − 4.5154 · X5 + 2.605086 · X1 · X4 −1.596 ≤ X1 ≤ 1.596 −1.596 ≤ X2 ≤ 1.596

(18)

(19)

−1.596 ≤ X3 ≤ 1.596 −1.596 ≤ X4 ≤ 1.596 −1.596 ≤ X5 ≤ 1.596 The selection of lower and upper bounds of the decision variables are based on the operating limits of the membrane reactors and Shu’s work [7]. In fact, they are the α value in the central composite design, but they were previously chosen with this purpose. The optimization led to the following values: X2 = X5 = −1.596

X1 = X3 = X4 = 1.596; Which codified back led to: P R0 = 506.625 kP a;

F CH40 = 1.0 × 10−5 mol/s; T W = 873.15 K;

F I = 1.375 × 10−5 mol/s

M = 2.5

The optimized methane conversion and hydrogen recovery calculated through the complete model using the above optimal values were 99.94% and 98.19%, which are higher than any other condition. 8. Conclusions The membrane reactor mathematical model developed showed a good agreement with the experimental data. The influence of five important parameters on methane conversion and hydrogen recovery was analyzed. Steam to methane feed flow ratio has an opposite effect when the objective is to maximize both methane conversion and hydrogen recovery. Optimization procedure was done in two steps. First, experimental design techniques were employed to adjust a surface response that shows how the target variables depend on optimizing parameters. Second, the optimum was found using the optimizer code NPSOL. It was shown that high methane conversion and hydrogen recovery, 99.94 and 98.19%, respectively, can be achieved using the following parameters values: PR0 = 506.625 kPa; FCH40 = 1.0×10−5 mol/s; FI = 1.375× 10−5 mol/s; TW = 873.15 K; M = 2.5. These preliminay optimized values, which are higher than any other condition, are consistent with the parametric analysis accomplished with the full model. 9. List of symbols a - Constant in Ergun’s equation [−] A0 - Pre-exponential factor of the permeability

h

N LT

T 2L M

n i

2 A1 - Heat transfer area between the reaction chamber and the external reactor temperature L 2 A2 - Heat transfer area between the reaction chamber and the permeate side L Am - Permeation area of membrane L2 b - Constant in Ergun’s equation [−] h i 2

Cvp - Specific heat at constant volume of the permeate side mixture TM2 NL Θ h i 2 Cvr - Specific heat at constant volume of the reaction side mixture TM2 NL Θ dp - Equivalent diameter particle [L] h of catalyst i E - Activation energy

M L2 T 2N

f - Friction factor [−] fH2 p - Dimensionless flow rate of hydrogen in the permeate side [−] fi - Dimensionless flow rate of i component in the reaction side [−] F CH40 - Methane feed flow rate N T FH2 p - Hydrogen flow rate in the permeate N T Fi - Molar flow rate of component i in the reaction side N T F I - Molar flow rate of inert gas [mol] h or sweep i

HH2 - Hydrogen enthalpy

M L2 T 2N

JH2 - Hydrogen permeation flux

N L2 T

L - Reactor length [L] M - Steam to methane feed flow H2 O e CH4 [−] ratio P R0 - Inlet reactor pressure TM 2L PR - Reactor pressure TM 2L Pr - Dimensionless reaction pressure [−] Pp - Permeate pressure TM 2L pH2 - Hydrogen partial pressure in the reaction side TM 2L pH2p - Hydrogen partial pressure in the permeate side TM 2L

h i 2 Q1 - Heat exchanged between the reaction side and external environment MT L2 i h 2 Q2 - Heat exchanged between the reaction side and the permeate side MT L2 i h Kmol rj - Rate of the reaction j Kg h cat h2 i R - Universal gas constant TM2 NL Θ

Re - Reynolds’ number [−] Tm - Membrane temperature [Θ] Tp - Dimensionless permeate temperature [−] Tr - Dimensionless temperature of the reaction side [−] TR - Temperature of the reaction side [Θ] Tw - External temperature [Θ] us - Gaseous velocity of the reaction chamber TL U1 - Overall heat transfer coefficient between reaction side and external side TM M3 Θ U2 - Overall heat transfer coefficient between reaction and permeation zone T 3 Θ vp - Mean velocity of the gaseous mixture in the permeate TL vr - Mean velocity of the mixture in the reaction chamber TL 3gaseous Van - Annular volume L Vpe - Permeate volume L3 XCH4 - Methane conversion [−] W - Catalyst mass [M ] YH2 - Hydrogen recovery [−] z - Dimensionless axial position [−] Z - Axial position [−] Greek letters δ - Membrane thickness [L] h i ∆Hj - Heat of the reaction j

M L2 T 2N

ε - Void fraction [−] νij - Stoichiometric coefficient of the component i in the reaction j [−] ρp - Permeate density LM3 ρr - Reaction density LM3 Subscripts i = CH4 , CO, H2 , H2 O, CO2 . j = 1, 2, 3 (reactions 1,2,3) k = CH4 , CO, H2 , H2 O. p = permeate side r = reaction side Acknowledgements This work has been supported by Brazilian funding agencies, CAPES and CNPq (Grant no 475934/2006-7). The authors are grateful for the assistance of M.Sc. Adriene Artiaga Pfeifer and Eng. Alice Medeiros de Lima in the use of DIRCOL and NPSOL codes.

10. References [1] Falco, M. D. T.; Di Paola, L.; Marrelli and Nardella, L. Simulation of large-scale membrane reformers by a two-dimensional model, Chemical Engineering Journal, 2006, 128, 115-125. [1] Froment, K. and Bischoff, K. B. Chemical Reactor Analysis and Design. John Wiley, New York, 1990. [2] Gallucci, F.; Paturzo, L.; Fama, A.; Basile, A. Experimental Study of the Methane Steam Reforming Reaction in a Dense Pd/Ag Membrane Reactor, Industrial and Engineering Chemistry Research , 2004, 43, 928-933. [3] Lin, Y.-M.; Liu, S.-L.; Chuanga, C.-H. and Chub, Y-T.. Effect of incipient removal of hydrogen through palladium membrane on the converion of methane steam reforming: experimental and modeling, Catalysis Today, 2003, 82, 127-139. [4] Ogden, J. M. Review of small stationary reformers for hydrogen production, Report No. IEA/H2/TR-02/002. Princeton, USA: Princeton University. 2001. [5] Ohmori, W. Y. T.; YU, W.; Yamamoto, A.; Endo, A.; Nakaiwa, M.; Hayakawa, T. and Itoh, N. Simulation of a porous ceramic membrane reactor for hydrogen production, International Journal of Hydrogen Energy, 2005, 30(10), 1071-1079, . [6] Poling, B. E.; Prausnits. M. J.; O’Connell, J. P. The Properties of Gases and Liquids, 5th ed., McGraw-Hill, 2004. [7] Shu, J.; Grandjean, P. A. and Kaliaguine, S. Methane steam reforming in asymmetric Pd- and Pd-Ag/porous SS membrane reactors, Applied Catalysis A: General, 1994, 119, 305-325. [8] Sjardin, M.; Damen, K. J. and Faaij, A. P. C. Techno-economic prospects of small-scale membrane reactors in a future hydrogen-fuelled transportation sector, Energy, 2006, 31, 2523-2555. [9] Tong, J.; Matsumura, Y.; Suda, H. and Haraya, K. Experimental Study of Steam Reforming of Methane in a Thin (6 µM) Pd-Based Membrane Reactor, Industrial and Engineering Chemistry Research, 2005, 44, 1454-1465. [10] Xu, J. and Froment, G. F. Methane Steam Reforming, Methanation and Water-Gas Shift: I. Intrinsic Kinetics, AIChE Journal , 1989, 35(1), 88-96. [11] Gill, P. E.; Murray, W.; Saunders, M. A.; Wright, M. H. User’s Guide for NPSOL (Version 4.0): A Fortran package for nonlinear programming, SOL 86-2, 1986.