Multiobjective Optimization of Industrial Autothermal Reformer for

Jul 15, 2009 - (RIPI), Tehran 14665-1998, Iran. The multiobjective optimization of an industrial autothermal reformer consisting of a noncatalytic par...
0 downloads 0 Views 2MB Size
Ind. Eng. Chem. Res. 2009, 48, 7529–7539

7529

Multiobjective Optimization of Industrial Autothermal Reformer for Syngas Production Using Nonsorting Genetic Algorithm II Alireza Behroozsarand,*,† Hadi Ebrahimi,‡ and Akbar Zamaniyan‡ Department of Chemical Engineering, Sahand UniVersity of Technology, Tabriz, 651335-1996, Iran, and Department of Natural Gas ConVersion, Gas Research DiVision, Research Institute of Petroleum Industry (RIPI), Tehran 14665-1998, Iran

The multiobjective optimization of an industrial autothermal reformer consisting of a noncatalytic partial oxidation (POX) chamber and a two-section catalytic steam reformer was studied in order to produce syngas with H2/CO ratio near 1 for application in Fischer-Tropsch and oxo processes. Both thermodynamic modeling of the POX section and mathematical modeling of fixed beds were used for the optimization of the process for maximizing of methane conversion and CO selectivity and minimizing of CO2 makeup flow. The nonsorting genetic algorithm II was employed for the optimization problem. The adjustable parameters are feed temperature, O2/CH4 ratio, steam to methane (S/C) ratio, and the first and second catalytic bed lengths of the reforming part. The optimized case proposes a syngas product with H2/CO of 1 by a feed including O2/CH4 and S/C of 0.59 and 2.012, respectively. The optimized lengths of two catalytic beds are 1.79 m for the first bed and 5.51 m for the second. 1. Introduction Synthesis gas, consisting of hydrogen and carbon monoxide, a useful gas, is employed in several chemical industries such as synthesis of Fischer-Tropsch (FT) for producing higher hydrocarbons (eq 1), methanol (eq 2), and ammonia (eq 3) and hydrogen production. Pure hydrogen can be produced from synthesis gas and the reaction between carbon monoxide and steam in the water gas shift reaction (eq 4) followed by undesired species removal. nCO + 2nH2 f CnH2n + nH2O

(1)

CO + 2H2 f CH3OH

(2)

3H2 + N2 f 3NH3

(3)

CO + H2O f CO2 + H2

(4)

In practice, synthesis gas is mostly produced by natural gas or methane steam reforming in tubular catalytic reactors, although this process is very expensive in large scale due to the high heat demand.1-3 Alternative routes to convert natural gas into synthesis gas have been studied by some authors including methane and CO2 reforming (eqs 5 and 6),4-6 partial methane oxidation (eq 7),7-9 and combination of them.10-14 The reaction network is very complex in all cases and includes several dominated reactions, presented in eqs 5-14 in Table 1.15 Steam reforming of methane is a preferred commercial process for production of syngas with H2/CO molar ratio of about 3.0;2 the process is highly endothermic.16 It is the most conventional and economical process for hydrogen and syngas production. Depending on the usage, the H2/CO ratio of the product is varied; it is mainly changed by CO2 injection into the feed. In this case, CO2 is converted to CO in the water gas * To whom correspondence should be addressed. Tel.: +98 (21) 44739551/9 2397. Fax: +98 (21) 44739537. E-mail: behroozsarandar@ ripi.ir. † Sahand University of Technology. ‡ Research Institute of Petroleum Industry.

shift reaction (eq 9). The other effective parameters are temperature and pressure of the system. So by controlling these parameters, the produced H2/CO can be adjusted. Besides conventional steam reforming, researchers have recently focused on noncatalytic and catalytic partial oxidation of methane, which would directly give the desired H2/CO ratio of 2.0: CH4 + 0.5O2 f CO + 2H2

∆H ) -22 kJ/mol, 1000 K (15)

CO2 reforming (eq 6), one of the alternatives to steam reforming of CH4, has received considerable attention recently for several potential advantages, including a lower H2/CO ratio (near 1), recovery of CO2, and a possible means of energy storage/transmission. The major problem encountered in CO2 reforming is rapid deactivation of the catalysts, mainly due to carbon deposition. Autothermal reforming (ATR) is a combination of noncatalytic partial oxidation and steam reforming developed by Haldor Topsoe in the late 1950s with the aim of reforming in a single reactor. In this reactor the preheated feed streams (CH4, H2O, and O2) are mixed in a burner placed at the top where the methane partial oxidation reactions occur (Figure 1). The steam reforming and water gas shift reactions take place in the catalyst bed below the burner. The catalyst bed, including two sub-beds, is designed for additional conversion, adjusting of H2/CO ratio, and controlling the equilibrium condition. The first bed has more fresh high-active catalyst than that of the second bed. Typically, the ATR operates at high temperatures, about 2200 K in the combustion zone and 1200-1400 K in the catalytic zone. The conditions result in low oxygen consumption (O2/CH4 ) 0.55-0.60); however, with a certain amount of steam is added to the feedstock to eliminate carbon formation. In all cases, the H2/CO ratio at the outlet of the reactor can be precisely adjusted by varying the H2O/CH4 and/or O2/CH4 molar ratio in the feed.1 For instance, addition of steam to the feed is caused an increase in H2/CO; however, it causes the decrease carbon formation. For the case of maximum H2 yield, some studies have been already performed for limited ATR reactor conditions for fuel cell applications; they are catalytic ATR and consist of a

10.1021/ie900259n CCC: $40.75  2009 American Chemical Society Published on Web 07/15/2009

7530

Ind. Eng. Chem. Res., Vol. 48, No. 16, 2009

Table 1. Reaction Network (Eqs 5-14) CH4 + H2O T CO + 3H2 CH4 + CO2 T 2CO + 2H2 CH4 + 0.5O2 f CO + 2H2 CH4 + 2O2 f CO2 + 2H2O CO + H2O T CO2 + H2 2CO f C + CO2 CH4 T C + 2H2 C + H2O T CO + H2 CO + 3H2 T CH4 + H2O CO2 + 4H2 T CH4 + 2H2O

∆H298K ∆H298K ∆H298K ∆H298K ∆H298K ∆H298K ∆H298K ∆H298K ∆H298K ∆H298K

catalytic POX.17-21 On the contrary the present work focuses on noncatalytic partial oxidation (NC-POX). It is assumed that the amine or PSA unit is a black box which separates carbon dioxide from sour syngas stream (Figure 1). Since 1968, various models have been developed for modeling of fixed bed reformer reactors: the one-dimensional homogeneous model,22-25 one-dimensional heterogeneous model,26,27 and two-dimensional heterogeneous model.28,29 In this work, the one-dimensional heterogeneous model, developed for reactor side by Xu and Froment,27 and already approved by some authors,29,30 has been linked to the model of POX section. In this paper, a typical industrial ATR is simulated using two models: an NC-POX thermodynamic model and a onedimensional (1D) pseudoheterogeneous model for the fixed bed catalytic sections (two fixed beds). The exothermic combustion reactions in the NC-POX section are very fast so that all oxygen is consumed by the reaction with methane. Methane combustion reaction occurred through many radical reactions. “A perfectly stirred reactor model (PSR)”31 is used for the NC-POX model. The present study considers a typical industrial autothermal reformer for analyzing the optimum product based on the detailed modeling of both POX and fixed bed sections in accordance with the industrial purposes (Figure 1). In the noncatalytic partial oxidation section, thermodynamic equilibrium analysis is performed by means of the Gibbs free energy minimization method. The possible species that might be formed in the final product gas are assumed to be CH4, O2, CO2, CO, H2O, and H2. The product yields are normalized on the basis of one initial unit mole of CH4.

Figure 1. Industrial autothermal reformer with CO2 separation scheme.

) ) ) ) ) ) ) ) ) )

206 kJ/mol (steam reforming) 247 kJ/mol (CO2 reforming) -36 kJ/mol (partial oxidation) -802 kJ/mol (total oxidation) -41 kJ/mol (water gas shift) -172 kJ/mol (Boudouard reaction) 75 kJ/mol (CH4 cracking) 130 kJ/mol (coal gasification) -206 kJ/mol (methanation) -160 kJ/mol (CO2 methanation)

(5) (6) (7) (8) (9) (10) (11) (12) (13) (14)

After modeling of the system, parameter analysis is employed for optimizing of the important parameters: maximizing of methane conversion and carbon monoxide selectivity besides the minimizing of makeup CO2 flow for a low H2/CO. Two possible methods may be applied for obtaining a low H2/CO ratio: • decreasing hydrogen production flow rate; • increasing carbon monoxide production flow rate. The first method may be performed by decreasing of feed steam to carbon ratio (S/C) that is not suitable, because of probability of carbon formation in low S/C. The second method, however, is more realistic. On the other hand, maximizing of methane conversion and carbon monoxide selectivity increases the flow rates of both synthesis gas and carbon monoxide. Some authors have already studied the optimization of autothermal reformers in various applications. Hagh20 studied the effect of important parameters (such as O2/CH4, steam/ carbon, heat loss) on the autothermal reformer for maximum hydrogen production. Babu and Rakesh32 considered optimal design of an autothermal reformer for ammonia synthesis. In the optimization section, they focused on the optimal length of the reformer for different top temperatures with the constraints of energy, mass balance of reaction, feed gas temperature, and mass flow rate of nitrogen for ammonia production. Ariane et al.15 focused on the optimization of reactions taking place in both primary and autothermal reformers. They have studied the optimization of the combined carbon dioxide reforming and partial methane oxidation over a 1% Pt/Al2O3 catalyst in order to produce synthesis gas with hydrogen/carbon monoxide ratio close to 1, for applications in metallurgical and polycarbonate pro-

Ind. Eng. Chem. Res., Vol. 48, No. 16, 2009

7531

Table 2. Coefficients of Eq 21 coefficients

A

CH4 CO CO2 H2 H2O

B

C -5

-8

-4.87 × 10 1.77 × 10-5 6.59 × 10-5 9.61 × 10-5 -4.13 × 10-5

1.0247 -0.0247 0.0118 -0.0590 0.0493

4.21 × 10 5.00 × 10-8 -2.51 × 10-7 -8.56 × 10-8 -4.07 × 10-8

cesses, production of oxygenated compounds, and hydrocarbons. It must be noted that most of these studies have been focused on catalytic autothermal reformers. In the present study a real parameter nonsorting genetic algorithm II (NSGAII) with blend crossover is used for optimizing the autothermal reformer parameters. The aim of study is maximizing of methane conversion and carbon monoxide selectivity and minimizing of CO2 makeup molar flow, while the hydrogen to carbon monoxide ratio is fixed near 1 for application in the Fischer-Tropsch and oxo syntheses. 2. Noncatalytic Partial Oxidation Model 2.1. Thermodynamic Equilibrium. Some authors have considered a simplified model of one molecular reaction, i.e., the highly exothermic combustion of CH4 to CO and H2O with an O2/CH4 ratio of 1.5 for the first part of the ATR: a noncatalytic POX chamber: CH4 + γO2 f CO + 2H2O

γ ) 1.5

(16)

33

Zhu et al. considered the reforming of the methane by partial oxidation for syngas production by means of thermodynamic analysis and kinetic simulation. They used the CHEMKIN package incorporating the GRI 1.2 mechanisms to predict the syngas compositions. The concept is employed here. In this work, the thermodynamic equilibrium analysis is performed by means of Gibbs free energy minimization method.34 Here, methane and oxygen are mixed and combustion occurred in the POX chamber. The products including CO, CO2, H2, and H2O with the relevant value (fi) is obtained to be entered to the catalytic fixed beds as the following equation: CH4 + γO2 f fCOCO + fCO2CO2 + fH2OH2O + fH2H2, O2 0.0 < γ ) e 2.0 CH4

(17)

where the fi (i ) CO, CO2, H2O, and H2) are constant coefficients. The equilibrium compositions (in an equilibrium temperature) are obtained using the direct minimization of the Gibbs function of the system, given by N

nG )

∑ n G¯i ) ∑ n G i

i

0 i

+ RT

i)1

fi

∑ n ln f i

(18)

0 i

For gas-phase reactions, ˆfi is equal to φˆ iyiP. Since the standard state is taken as the pure ideal gas state at the pressure of 101.32 kPa or f 0i ) 101.32 kPa, and since G0i is set equal to zero for each chemical element in its standard state, G0i ) ∆Gfi0 for each component. Rephrasing the equation (eq 18) gives nG(nis, T, P) )

∑ n ∆G i

0 fi

∑ n RT ln P + ∑ n RT ln y + ∑ n RT ln φˆ +

i

i

i

i

i

(19)

D

E

F

G

H

-4.2691 1.1017 0.1115 2.0684 0.3550

6.9971 -0.7122 -0.4693 0.0005 -2.3540

-5.5217 0.0873 0.7353 -3.4697 5.1032

2.1051 0.0735 -0.3877 2.6236 -3.2459

-0.3112 -0.0215 0.0690 -0.5636 0.6574

The set of nis which minimizes nG at constant T and P are solved by the constraints of elemental balances:35 N

∑na

i ji

) bj

j ) 1, ..., K

(20)

i)1

The objective function (eq 19) is minimized, and the constraints (eq 20) are solved by the Newton-Raphson and Lagrange multiplier methods.15,35 After obtaining of partial oxidation products by Gibbs free energy minimizing (solving of eqs 19 and 20, simultaneously), data is imported to Table Curve v1.0 software, and the following polynomial is obtained based on the best equation preferred by the software (therefore, the normalized product yields (fi) are calculated as a function of feed temperature (T) and O2/CH4 ratio (γ)): fi ) A + BTFeed + CTFeed2 + Dγ + Eγ2 + Fγ3 + Gγ4 + Hγ5 (21) Use of feed temperature as an input variable is more convenient than that of equilibrium or output temperature since they are not accessible in a simulation problem. Constant coefficients A, B, C, D, E, F, G, and H are prepared in Table 2 for all product components. The results of thermodynamic analysis according to the above representation are depicted in Figure 2. In this figure the normalized yields of compositions versus oxygen to methane ratio and temperature are presented for adiabatic conditions (Figure 2, parts a and b). Previous authors have mostly focused on the isothermal conditions; however, in the real industrial cases the adiabatic conditions may be more favorable, because the reformer output temperature is not adjustable. There are also some authors that use the equilibrium temperature and pressure for presenting of the product yields.33,36,37 But the inlet stream temperature is more controllable in the case of adiabatic conditions. Figure 2a shows normalized yields of the components and as a function of feed O2/CH4 ratio. Figure 2b shows the effect of pressure on normalized H2 yield as a function of temperature. According to Figure 2a, maximum hydrogen production is obtained at an O2/CH4 ratio of 0.5. Beside this, in higher ratio of O2/CH4 (>0.8), the methane is converted completely. The normalized CO yield is almost constant at O2/CH4 > 1, but the normalized H2O yield is significantly increased. Figure 2b presents that increasing of the reformer pressure in low feed temperatures decreases hydrogen production, but there is not considerable variation in high feed temperature that is common in industrial cases. Therefore, pressure effect on stoichiometric coefficients in high temperature may be negligible. 2.2. Mass and Energy Balances. The noncatalytic partial oxidation process was simulated by a “perfectly stirred reactor model” (PSR).38 The PSR model corresponds to the perfect mixing between reactants and products. On the other hands, PSR represents a reactor where perfect mixing of reactants and products is assumed. It is characterized by spatial homogeneity and steady-state operation. The mixing is instantaneous, and

7532

Ind. Eng. Chem. Res., Vol. 48, No. 16, 2009

Figure 2. (a) Normalized components yield for adiabatic condition (P ) 36 bar, T ) 670 °C); (b) effect of the total pressure on H2 yield as a function of feed temperature for adiabatic condition.

the conditions in the outlet of the reactor are equal to the conditions in the reactor. The steady-state governing equations for the gas-phase species include31 ω ˙ iWi 1 (Y - Yi) + ) 0, τ i,in F

i ) 1, ..., N

(22)

and the gas energy equation

∑Y

(I) SR:

(II) WGS:

N

m ˙

presented in this section, are applied for both fixed beds. Xu and Froment27 considered three chemical reactions of steam reforming (SR), water gas shift (WGS), and methanation (MR) reactions:

N

i,in(hi,in

- hi) - V

i

∑ h ω˙ W i

i

i

)0

(23)

i)1

In eqs 22 and 23, τ and V are the nominal residence time in the reactor and the reactor volume, respectively. The subscript “in” indicates inlet stream quantities. The nominal residence time in the reactor is given by FV m ˙

τ)

j ) W

j pW RT

(25)

(∑ ) N

i)1

Yi Wi

r1

(28)

CO + H2O T CO2 + H2

r2

(29)

CH4 + 2H2O T CO2 + 4H2

-1

(26)

r1 )

(k1 /PH22.5)(PCH4PH2O - PH23PCO /K1) DEN2

r2 )

r3 )

(k2 /PH22)(PCOPH2O - PH2PCO2 /K2) DEN2

(k3 /PH23.5)(PCH4PH2O2 - PH24PCO2 /K3) DEN2



T

T0

Cp,i(T) dT

)Ω

(27)

3. Catalytic Fixed Bed Model 3.1. Rate Equations. After consuming all O2 by the combustion in the POX section, the stream mainly consisting of CO, CO2, H2, H2O, and the reminded (unreacted) CH4 is entered to the first fixed bed of the catalytic section, followed by the second one, under the first bed (Figure 1). The governing equations,

(31)

(32)

(33)

3.2. Material Balances. Mass transfer occurs in two sections, the gas phase and catalyst particles. Mass balance equations in the gas phase for two components, CH4 and CO2, are13

dy hi ) hi,o +

(30)

DEN ) 1 + KCOPCO + KH2PH2 + KCH4PCH4 + KH2OPH2O/PH2 (34)

dxCH4

In eq 23 hi is

r3

with the corresponding intrinsic rate equations27 of

(24)

The eqs 22 and 23 are solved simultaneously to obtain the temperature and compositions of the product. The density is calculated from the ideal gas law, and the enthalpy is related to the temperature through the constitutive equation. These two relations are given in eqs 25 and 26. The residence time, the pressure, and the inlet species mass fractions, and temperature are used as the input data. F)

(III) MR:

CH4 + H2O T CO + 3H2

dxCO2 dy

)Ω

FBηCH4rCH4 FCH4◦ FBηCO2rCO2 FCH4◦

(35)

(36)

These equations can be calculated by the finite difference method along a differential element of tube length. The effectiveness factor η is calculated from the heterogeneous model.27,39

Ind. Eng. Chem. Res., Vol. 48, No. 16, 2009

3.3. Energy Balance. The energy balance for an axial element of dy of the reactor is39 FBAcross

∑ [(∆H )r η ] dy + [F Cp u A j

j j

g

j

g s crossTipg]net

+

UaAside(Ts,k - Tipg) ) 0

(37)

Subscript “ipg” refers to the element i of process gas in the tube. This equation includes three terms of heat generated by chemical reactions, convection energy transferred by fluid, and energy transferred from tube skin, respectively. Ua is overall heat transfer coefficient which includes both tube material and catalyst bed resistances:

()

dte dte 1 1 ) ln + Ua 2λt dti R

(38)

R is the convective heat transfer coefficient in the catalyst bed which is calculated using correlations developed by Xu and Froment.27 More details of the fixed bed catalytic part have been considered previously in the literature.27,30,39 3.4. Momentum Equation. The momentum equation is40 -

FguS2 dP )f dy gdp

(39)

And friction factor is f)

1 - ε 1.74 + 150(1 - ε) ε Redp

(40)

4. Genetic Algorithm for Multiobjective Optimization 4.1. Genetic Algorithm (GA). As the new methods, genetic algorithms have been successfully applied to nonlinear optimization problems in many dimensions, where traditional methods are often found to fail.41 Moreover, more traditional ones such as deterministic and gradient-based optimization methods do not search the parameter space and can tend to converge toward a local extreme of the fitness function. It is clearly unsatisfied for problems where the fitness varies nonmonotonously with parameters. On the other hand, genetic algorithms are able to depart from local optima due to the variability of the parameters within the “gene pool” and the element of randomness inherent within the methods. Furthermore, genetic algorithms do not require knowledge of the gradient of the fitness functions, which makes them particularly suited to optimization problems for which an analytical expression is not known for the fitness functions. 4.2. Multiobjective Optimization. There are two general approaches to multiple-objective optimization. One is to combine the individual objective functions into a single composite function or move all but one objective to the constraint set. In the former case, determination of a single objective is possible with methods such as utility theory, weighted sum method. But the problem lies in the proper selection of the weights or utility functions to characterize the decision-maker’s preferences.42 The second general approach is to determine an entire Pareto-optimal solution set or a representative subset. A Pareto-optimal set is a set of solutions that are nondominated with respect to each other. While moving from one Pareto solution to another, there is always a certain amount of sacrifice in one objective(s) to achieve a certain amount of gain in the other(s). Pareto-optimal solution sets are often preferred to single solutions because they can be practical when considering real-life problems since the

7533

final solution of the decision-maker is always a trade-off. Paretooptimal sets can be of varied sizes, but the size of the Pareto set usually increases with the increase in the number of objectives. The ultimate goal of a multiobjective optimization algorithm is to identify solutions in the Pareto-optimal set. However, identifying the entire Pareto-optimal set, for many multiobjective problems, is practically impossible due to its size. In addition, for many problems, especially for combinatorial optimization problems, proof of solution optimality is computationally infeasible. Therefore, a practical approach to multiobjective optimization is to investigate a set of solutions (the best-known Pareto set) that represent the Pareto-optimal set as well as possible. With these concerns in mind, a multiobjective optimization approach should achieve the following three conflicting goals:43 1. The best-known Pareto front should be as close as possible to the true Pareto front. Ideally, the best-known Pareto set should be a subset of the Pareto-optimal set. 2. Solutions in the best-known Pareto set should be uniformly distributed and diverse over of the Pareto front in order to provide the decision-maker a true picture of trade-offs. 3. The best-known Pareto front should capture the whole spectrum of the Pareto front. This requires investigating solutions at the extreme ends of the objective function space.42 4.3. Multiobjective GA. Being a population-based approach, GA is well-suited to solve multiobjective optimization problems. A generic single-objective GA can be modified to find a set of multiple nondominated solutions in a single run. The ability of GA to simultaneously search different regions of a solution space makes it possible to find a diverse set of solutions for difficult problems with nonconvex, discontinuous, and multimodal solutions spaces. 4.4. Model Equations for Process Parameters Optimization. The optimization problem is considered for an industrial retrofit plant. In this optimization, methane flow rate is constant. In the previous section, at O2/CH4 ratio of 0.55, the methane conversion and CO selectivity were maximized and makeup CO2 molar flow was minimized. It can logically search for operating scenarios that will maximize methane conversion and CO selectivity and minimize makeup CO2 molar flow simultaneously. Performing a constrained optimization with both of them as objectives can identify such scenarios. The optimization problem can be expressed mathematically as the following: Maximize object(1) ) XCH4 (maximize methane conversion) (41) object(2) ) SCO (maximize CO selectivity)

(42)

Minimize object(3) ) FCO2 (minimize makeup CO2 molar flow) (43) that XCH4 )

molar flow HCin - molar flow HCout molar flow HCin

(44)

7534

Ind. Eng. Chem. Res., Vol. 48, No. 16, 2009

SCO )

molar flow COout

∑ (molar flow of j)

j ) CH4, H2, CO2, H2O

out

(45)

subject to 600 °C e TFeed e 900 °C

(46)

0.4 e O2 /CH4 e 0.6

(47)

2.0 e S/C e 3.0

(48)

0.1 m < LB1 < 2.5 m

(49)

2.5 m < LB2 < 6.0 m

(50)

FCH4 ) constant

(51)

Table 3. Feed Characteristics of Industrial ATR temperature (°C) pressure (kPa) molar flow, kmol/h compositions, mol % CO2 CO H2 CH4 N2 H2O CH3OH C2H6 C3H8 O2

2.73 1.44 14.66 26.79 1.38 37.96 0.01 0.08 0.04 14.91

Table 4. Comparison of Experimental and Calculated Results for Output Syngas of Industrial ATR variable

LB1 and LB2 are the first and second catalytic fixed bed, respectively. The wide range of the feed temperature constraint of the autothermal reformer has been selected based on the reformer wall resistance temperature and lifetime of Ni catalyst at catalytic fixed bed sections. According to Figure 3 at low temperature (below 600 °C) both two first objective functions (CH4 conversion and CO selectivity) are decreased; therefore, the lower limit, 600 °C, is proper. At an O2/CH4 ratio between 0.4 and 0.6 the objectives have maximum value if other parameters such as S/C, TFeed, and FCH4 are constant. And, according to the literature and industrial steam reformers, the lower limit of S/C is set at 2.0 to avoid carbon formation on the catalyst,1 which occurs at low value of S/C. Very high S/C ratio affects the process economics adversely because of the energy requirement for heating up the feed. A fixed feed flow rate of methane is assumed to optimize the retrofit syngas production unit. The limits of LB1 and LB2 were defined based on design information.

670 3600 23850

T (°C) P (kPa) flow (kmol/h) compositions, mol % CO2 CO H2 CH4 H2O N2

calculated 973.3 3500 32569.8 0.063 0.156 0.46 0.01 0.301 0.01

experimental 975 3500 32420 0.057 0.161 0.483 0.013 0.276 0.01

relativea error %

absoluteb error

0.18 0 0.46 9.86 3.4 5.13 0.23 8.17

0.006 0.005 0.024 0.002 0.025

a Relative error (%) ) ((experimental - calculated)/experimental) × 100. b Absolute error ) (experimental - calculated).

5. Results and Discussion 5.1. Model Verification. Table 3 presents the feed data of industrial autothermal reformer. Feed compositions include feed stream, steam, and oxygen. Although the natural gas is used for the simulation as an industrial point of view, it is assumed that methane is a dominant component, and only methane is used for the calculation. For simplification, all other hydrocarbons such as heavier ones and methanol are converted to

Figure 4. Effect of oxygen/methane ratio at feed on CH4 conversion and CO selectivity; TFeed ) 850 °C, FCH4 ) 4000 kmol/h, S/C ) 2.5, H2/CO ) 2.0.

Figure 3. Effect of feed temperature on CH4 conversion and CO selectivity; FCH4 ) 4000 kmol/h, O2/CH4 ) 0.55, S/C ) 2.5, H2/CO ) 2.0.

methane. On the other hand, the equaled methane is named methane, consisting of all rich-methane hydrocarbons, existed in the natural gas feed. Table 4 shows syngas product obtained from both an industrial plant and simulated results. Both relative and absolute errors are reported. The relative error for CO2 mole fraction is 9.86 due to its low order of magnitude; therefore, the absolute error should be considered that it is 0.006, an acceptable value. It shows that the model is able to predict the industrial reformer performance with high accuracy. Temperature values of main points in the ATR reformer are depicted in Figure 7. This figure shows that a sharp increase in the temperature (almost 60%) due to the high exothermic oxidation reactions of the combustion chamber occurred in the

Ind. Eng. Chem. Res., Vol. 48, No. 16, 2009

Figure 5. Effect of steam/methane ratio of feed on CH4 conversion and CO selectivity; TFeed ) 850 °C, FCH4 ) 4000 kmol/h, O2/CH4 ) 0.55, H2/ CO ) 2.0.

POX section. The temperature is increased again from the bottom of the POX to the outlet of the first catalytic bed. Then, a 0.2% decrease in temperature occurred along the second bed of the catalytic section. For more considering, the temperature profile along the catalytic fixed beds is depicted in Figure 8. The figure shows

7535

that temperature is increased sharply at the entrance of the first catalytic bed. This is caused by high reaction rates when the amount of reactants is very much. After reaching a temperature pick (almost 968 °C), the profile is decreased immediately in the first bed. Because of the catalyst characteristic, the rate of decrease in the second bed is lower than in the first bed. The fresh catalyst with high activity is located in the first bed. 5.2. Analysis of Effective Parameters. In this section the effect of process variables on the two objective functions (only XCH4 and SCO) are considered and analyzed. 5.2.a. Feed Temperature. Figure 3 presents the effect of feed temperature on methane conversion and carbon monoxide selectivity. It is obvious that both methane conversion and CO selectivity are increased with increasing of the feed temperature as an adjustable input variable. The maximum values are obtained in a high feed temperature. For instance, the value of 1200 °C for feed temperature should be applied in the special industrial case according to Figure 3. There is a fact that preparation of high heat supply for feed stream is undesirable, because of economic aspects, and probably catalyst destruction may be occurred. It is noted that other variables such as O2/ CH4 and S/C ratios are constant during the simulation. Therefore, very high feed temperature is not sufficient parameter for the process optimization.

Figure 6. Effect of O2/CH4 ratio on H2/CO ratio at syngas (a) TFeed ) 850 °C, FCH4 ) 4000 kmol/h, S/C ) 2.5; effect of feed temperature ratio on H2/CO ratio at syngas (b) FCH4 ) 4000 kmol/h, S/C ) 2.5, O2/CH4 ) 0.55.

Figure 7. Temperature values of main points in the ATR reformer.

7536

Ind. Eng. Chem. Res., Vol. 48, No. 16, 2009

Figure 8. Temperature profile of the catalytic fixed bed of the ATR. Table 7. Nondominated Pareto-Optimal Solutions after 50 Generations

Table 5. Effect of an Increase in the Decision Variable on the Objective Functions effect of increase in objective function

TFeed

O2/CH4

steam/CH4

LB1

LB2

XCH4(%) SCO (%) FCO2(%)

v v v

v v V

v V Vv

V V constant

v constant constant

Table 6. Governing Parameters of NSGAII parameter name

method and value

no. of decision variables no. of objectives population size maximum generations replace proportion tournament pool size crossover method crossover probability mutation method mutation probability

5 3 50 50 0.9 2 two points 0.9 selective 0.1

5.2.b. O2/CH4 Ratio of Feed. Figure 4 shows the effect of oxygen/methane ratio of feed on CH4 conversion and CO selectivity. By increasing molar flow of oxygen in reformer feedstock, the required oxygen for methane or natural gas combustion is increased; therefore, methane conversion will be increased. The trend of CO selectivity consists of two sections. By increasing the O2/CH4 ratio between 0 and 0.55, CO selectivity reaches a maximum point and then is decreased. Therefore, an oxygen/methane ratio of 0.55 is the best value for obtaining maximum value of CO selectivity. According to Figure 2, oxygen/methane ratio has significant effect on normalized H2O yield. In high oxygen/methane ratio, it is increased rapidly. Increasing H2O molar flow in the output NC-POX section causes that the SR and WGS reactions in the catalytic section occurred in the forward direction. Therefore, the molar flow of hydrogen is increased more than CO molar flow and CO selectivity is decreased subsequently. 5.2.c. Steam/CH4 Ratio. Figure 5 presents the effect of S/C of feed on CH4 conversion and CO selectivity. It is clear that the effect of S/C on methane conversion is very low. Therefore, the S/C parameter is not suitable for maximizing methane conversion. According to Figure 6, increasing the S/C ratio decreases significantly CO selectivity because of increasing steam molar flow in the feed; the SR and WGS reactions in the catalytic bed tend to follow the forward direction. Therefore, H2, CO, and CO2 molar flows in the syngas product are

no.

TFeed (°C)

O2/CH4

S/C

LB1 (m)

LB2 (m)

1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49 50

898.33 898.33 898.33 898.33 898.33 898.33 898.33 898.33 898.33 898.33 898.33 898.33 898.33 898.33 898.33 898.33 898.33 898.33 898.33 898.33 898.33 898.33 898.33 898.33 898.33 898.33

0.59739 0.57767 0.58092 0.58092 0.59739 0.59739 0.59739 0.58448 0.59739 0.59739 0.59739 0.59739 0.58448 0.58448 0.59739 0.59739 0.59739 0.59739 0.59739 0.59739 0.59739 0.59739 0.59739 0.59739 0.59739 0.59739

2.0198 2.1299 2.0244 2.0244 2.0198 2.0198 2.0198 2.0229 2.1299 2.0198 2.0198 2.0198 2.0244 2.0229 2.0198 2.0198 2.0198 2.0198 2.0198 2.0198 2.0198 2.0198 2.0198 2.0198 2.0198 2.0198

1.7889 0.99667 2.0233 2.4566 2.4566 1.7889 1.7889 0.6149 1.7889 1.7889 1.7889 1.7889 1.7889 0.6149 1.7889 2.4566 1.7889 1.7889 1.7889 1.7889 1.7889 1.7889 2.4566 1.7889 1.7889 1.7889

5.5057 3.351 2.9482 5.448 5.448 5.5057 5.5057 3.377 5.5057 5.5057 5.5057 5.5057 5.5057 3.377 5.5057 5.448 5.5057 5.5057 5.5057 5.5057 5.5057 5.5057 5.448 5.5057 5.5057 5.5057

makeup XCH4 SCO CO2 (kmol/h) 98.1 96.9 97.2 97.3 97.6 97.8 98 97.7 97.6 98 98 98.1 97.6 97.7 97.6 98 97.9 98 98 98.1 97.6 98.1 98.1 98.1 98.1 98.1

25 22 24 24 24 25 25 25 23 25 25 25 25 25 24 25 25 25 25 25 24 25 25 25 25 25

82686 30767 50893 59315 46248 69221 79629 76563 46248 81018 80189 82686 74397 76801 46248 81018 74726 76807 81018 81994 46248 81994 81994 81994 81994 81994

increased. However, because of the higher stoichiometric coefficient of hydrogen in the SR reaction, hydrogen molar flow increases more than CO in the syngas product. 5.2.d. H2/CO Ratio of Syngas. Figure 6 presents the effect of two important parameters on H2/CO ratio in the syngas stream. It is found that feed temperature and O2/CH4 ratio have significant effect on the H2/CO ratio. According to Figure 6, low H2/CO is obtained at low feed temperature and high O2/ CH4 ratio. In fact the fixed bed CH4 entrance would be decreased by increasing of feed O2/CH4 ratio in the feed. Also, the POX product molar flows such as CO, CO2, and H2 are increased in this case. Therefore, both CO and H2 molar flows are decreased, in contrast to H2O, based on the SR reaction. Also CO and H2O molar flow are increased by contrary with the molar flow of H2 according to the WGS reaction. Finally the H2/CO ratio is decreased after catalytic fixed bed section. 5.3. Optimization Results. The optimization problem involves the three mentioned objective functions, which have different behaviors.

Ind. Eng. Chem. Res., Vol. 48, No. 16, 2009

7537

Table 8. Comparison of Operating Parameters for Three Optimal Cases (FCH4 ) 5000 kmol/h)

Figure 9. Pareto-optimal set of solutions obtained for the simultaneous optimization of XCH4, SCO, and makeup CO2 molar flow.

Table 5 shows the effect of variation in the decision parameters on the objective functions. It has been generated based on a parametric sensitivity analysis of the simulated mathematical model. For instance, increasing of feed temperature (TFeed) increases methane conversion, CO selectivity, and CO2 makeup molar flow (CO2 makeup is adjusted by H2/CO ) 1). The objective functions were optimized fulfilling the constraints and bounds given in eqs 46-50. As stated earlier the NSGAII algorithm was used for obtaining the Pareto-optimal solutions. A MATLAB code for real-parameter NSGAII described in an earlier section has been used to optimize the parameters. A population size of 50 was chosen with crossover of 0.9 and mutation probability of 0.1. The input parameters of NSGAII are given in Table 6.

parameters

chromosome A

chromosome B

chromosome C

TFeed (°C) O2/CH4 S/C LB1 LB2 XCH4 (%) SCO (%) FCO2 (kmol/h) H2/CO

898.33 0.59 2.012 1.79 5.51 98.1 25 82686 1.00

898.33 0.58 2.024 2.02 2.95 97.2 24 50893 1.00

898.33 0.57 2.130 1.00 3.35 96.9 22 30767 1.00

The different operations were performed for 50 generations to obtain the nondominated Pareto-optimal solutions. The Pareto-optimal solution sets after 50 generations are shown in Table 7. Figure 9 shows the Pareto set of optimal solutions obtained for the problem formulated above. The values of XCH4, SCO, and FCO2 are plotted. The points in the Pareto set (Figure 9) indicate the maximum possible of XCH4, SCO or minimum possible of FCO2 with the given operating constraints. The benefit of a multiobjective optimization is evident upon observing the wide choice of operating points available in the Pareto-optimal set. The CPU time taken to generate one set of Pareto-optimal solutions (such as those in Figure 9) is 300 min on the PCPentium (R) 4 CPU 3.00 GHz, 512 MB of RAM. Each point (referred to as a chromosome) on the Pareto set (Figure 9) is associated with a set of decision variables. Figure 10 presents the plot of decision variables corresponding to each of the points on the Pareto set. Table 8 lists values of the decision variables corresponding to three chromosomes A, B, and C selected from the Pareto set in Figure 9. In fact, a tradeoff between maximizing methane conversion, CO selectivity, and minimizing carbon dioxide makeup molar flow was employed to select the three solutions (chromosomes), prepared in Table 8. H2/CO ratio near 1 of the syngas product is an

Figure 10. Decision variables corresponding to each of the Pareto-optimal solutions shown in Figure 9.

7538

Ind. Eng. Chem. Res., Vol. 48, No. 16, 2009

important parameter which has been considered in all cases. The best choice can be selected based on the other considerations such as mechanical and especially economical aspects. According to Table 8 case A seems the best based on the proper H2/ CO (equal 1). In this case, the O2/CH4 of 0.59 and S/C of 2.012 are favorable in contrast with the high feed temperature that is 898.33 °C (it needs more heat load). 6. Conclusion An autothermal reformer was modeled and optimized in this paper. In partial oxidation and fixed catalytic sections, “a perfectly stirred reactor model (PSR)” and one-dimensional pseudoheterogeneous model are used, respectively. For kinetics of partial oxidation, a correlation between normalized yield of products and feed parameters is derived based on the Gibbs free energy minimization theory. The simulation results have good agreement with experimental data. By analyzing process parameters, the trend and constraint of adjustable variables are obtained. The process is optimized for obtaining H2/CO ratio near 1, maximizing methane conversion, and CO selectivity. Results show that feed temperatures between 850 and 950 °C, O2/CH4 between 0.5 and 0.6, S/C between 2.0 and 2.5, first catalytic bed length between 1.5 and 2.0 m, and second catalytic bed length between 5.0 and 6.0 m are optimum values for industrial cases. The benefit of a multiobjective optimization is evident upon observing the wide choice of operating points available in the Pareto-optimal set. Between four cases, case A has maximum methane conversion and carbon monoxide selectivity at low H2/CO ratio (near 1). It is concluded that the unified modeling of two parts of the autothermal reformer is useful for optimization of any industrial autothermal reformer. Nomenclature Acronyms 1D ) one dimensional ATR ) autothermal reformer FT ) Fischer-Tropsch NC-POX ) noncatalytic partial oxidation NSGA ) nonsorting genetic algorithm POX ) partial oxidation PSR ) perfectly stirred reactor syngas ) synthesis gas WGS ) water gas shift Latin across ) cross section of the fixed bed section, m2 A, B, ..., H ) constant of eq 7 Cp,i ) constant pressure heat capacity of species i, kJ · kg-1 · K-1 dP ) catalyst particle diameter, m dte ) tube inside diameter, m dto ) tube outside diameter, m fCO ) carbon monoxide stoichiometric coefficient fCO2 ) carbon dioxide stoichiometric coefficient fH2O ) steam stoichiometric coefficient fH2 ) hydrogen stoichiometric coefficient fCH04 ) CH4 molar flow of feed fCO20 ) CO2 molar flow of feed FCO2 ) makeup CO2 molar flow g ) gravity acceleration, m/s2 Gi ) Gibbs free energy of pure species i at operating conditions j i ) average Gibbs free energy of pure species i at operating G conditions

∆Gfi0 ) standard Gibbs free energy of formation of species i hi ) enthalpy of species i, kJ · kg-1 K ) total number of atomic elements Kij ) binary interaction parameters LB1 ) length of first section catalytic bed, m LB2 ) length of second section catalytic bed, m N ) total number of species in the reaction mixture ni ) number of moles of species i nG ) total Gibbs free energy of the system m ˙ ) mass flow rate, kg · s-1 P ) pressure, kPa R ) universal gas constant, kg · kmol-1 · K-1 ri ) rate of reaction for component j, kmol · kg-1 · s-1 SCO ) CO selectivity T ) temperature, K Ua ) overall heat transfer coefficient Us ) gas velocity, m/s V ) volume of PSR, m3 Wi ) molecular weight of species i, kg · kmol-1 j ) average molecular weight, kg · kmol-1 W xi ) conversion (i ) CH4 and CO2) Xi ) process variables (optimization) XCH4 ) methane conversion Y ) mass fraction of species i Z ) compressibility factor Greek R ) convective heat transfer coefficient γ ) oxygen/methane ratio 0 ∆H298 ) molar enthalpy, kJ · mol-1 ε ) catalyst porosity ηi ) effectiveness factor of reaction i λt ) tube conductivity coefficient F ) mixture density, kg · m-3 FB ) bulk density of catalytic bed, kgcat · mbed-3 Fg ) density of process gas, kg · m-3 τ ) nominal residence time in the reactor, s φˆ i ) fugacity coefficient of species i in the gas mixture Ω ) fixed bed reactor cross section, m2 ω ˙ i ) molar rate of production of species i, kmol · m-3 · s-1 ω ) acentric factor

Literature Cited (1) Pena, M. A.; Go´mez, J. P.; Fierro, J. L. G. New catalytic routes for syngas and hydrogen production. Appl. Catal., A 1996, 144, 7. (2) Yong, L. u.; Jinzhen, X. e.; Changchun, Y.; Yu, L.; Shikong, S. Mechanistic investigations on the partial oxidation of methane to synthesis gas overa nickel-on-alumina catalyst. Appl. Catal., A 1998, 174, 121. (3) Haitham, a. Q. Effect of ageing on a steam reforming catalyst. Chem. Eng. J. 1997, 66, 51. (4) Bitter, J. H.; Seshan, K.; Lercher, J. A. Mono and bifunctional pathways of CO2/CH4 reforming over Pt and Rh based catalysts. J. Catal. 1998, 176, 93. (5) Bradford, M. C. J.; Vannice, M. A. CO2 reforming of CH4. Catal. ReV. Sci. Eng. 1999, 41, 1. (6) Michael, C. J.; Bradford, M. C. J.; Vannice, M. A. CO2 reforming of CH4 over supported Ru catalysts. J. Catal. 1999, 183, 69. (7) Yong, L.; Jinzhen, X.; Changchun, Y.; Yu, L.; Shikong, S. Mechanistic investigations on the partial oxidation of methane to synthesis gas over a nickel-on-alumina catalyst. Appl. Catal., A 1998, 174, 121. (8) Kiyoharu, N.; Kengo, A.; Naoko, M.; Naoki, I.; Toshimitsu, S.; Yonghong, T.; Tetsuhiko, K.; Masatake, H. Effect of support on the conversion of methane to synthesis gas over supported iridium catalysts. Catal. Lett. 1998, 51, 163. (9) Aisling, M. O.; Julian, R. H. R. The effect of O2 addition on the carbon dioxide reforming of methane over Pt/ZrO2 catalysts. Catal. Today 1998, 46, 203.

Ind. Eng. Chem. Res., Vol. 48, No. 16, 2009 (10) Vernon, P.; Green, M.; Cheetham, A.; Ashcroft, A. T. Partial oxidation of methane to synthesis gas, and carbon dioxide as an oxidising agent for methane conversion. Catal. Today 1992, 13, 417. (11) Matsumura, Y.; Moffat, J. B. Oxidative methane conversion to carbon monoxide and hydrogen at low reactor wall temperatures over ruthenium supported on silica. Catal. Lett. 1994, 24, 59. (12) Choudhary, V. R.; Rajput, A. M.; Prabhakar, B. Energy efficient methane-to-syngas conversion with low H2/CO with carbon dioxide and oxygen. Catal. Lett. 1995, 32, 391. (13) Choudhary, V. R.; Uphade, B. S.; Mamman, A. S. Simultaneous steam and CO2 reforming of methane to syngas over NiO/MgO/SA-5205 in presence and absence of oxygen. Appl. Catal., A 1998, 168, 33. (14) Ruckenstein, E.; Hu, Y. H. Combination of CO2 reforming and partial oxidation of methane over Ni/MgO catalyst. Ind. Eng. Chem. Res. 1998, 37, 1744. (15) Ariane, L. L.; Neuman, S. d. R.; Vera, M. M. S.; Jose´, C. P. Modeling and optimization of the combined carbon dioxide reforming and partial oxidation of natural gas. Appl. Catal., A 2001, 215, 211. (16) Renesme, G.; Saint, J. J.; Muller, Y. Transportation fuels and chemicals directly from natural gas: how expensive? Catal. Today 1992, 13, 371. (17) Chan, S. H.; Wang, H. M. Thermodynamic analysis of natural-gas fuel processing for fuel cell applications. Int. J. Hydrogen Energy 2000, 25, 441. (18) Doss, E.; Kumar, R.; Ahluwalia, R. K.; Krumpelt, M. Fuel processors for automotive fuel cell systems: a parametric analysis. J. Power Sources 2001, 102, 1. (19) Ahmed, S.; Krumpelt, M. Hydrogen from hydrocarbon fuels for fuel cells. Int. J. Hydrogen Energy 2001, 26, 291. (20) Hagh, B. F. Optimization of autothermal reactor for maximum hydrogen production. Int. J. Hydrogen Energy 2003, 28, 1369. (21) Chan, S. H.; Hoang, D. L.; Ding, O. L. Transient performance of an autothermal reformersA 2-D modeling approach. Int. J. Heat Mass Transfer 2005, 48, 4205. (22) Davies, J.; Lihou, D. Optimal design of methane steam reformer. Chem. Proc. Eng. 1971, 52, 71. (23) Golebiowski, A.; Walas, T. Thermal process in catalytic reforming of methane with water vapor. Int. Chem. Eng. 1973, 13, 133. (24) Hyman, M. Simulate methane reformer reactions. Hydrocarbon Process. 1968, 47, 131. (25) Singh, C. P. P.; Saraf, D. N. Simulation of side fired steam hydrocarbon reformers. Ind. Eng. Chem. Process Des. DeV. 1979, 18, 30. (26) Elnashaie, S.; Adris, A. M.; Soliman, M. A.; Al-Ubaid, A. S. Digital simulation of industrial steam reformers for natural gas using heterogeneous models. Can. J. Chem. Eng. 1992, 70, 786. (27) Xu, J.; Froment, G. F. Methane steam reforming: IIsDiffusional limitations and reactor simulation. AIChE J. 1998, 35, 97.

7539

(28) Ferreira, R. M. Q.; Marques, M. M.; Babo, M. F.; Rodrigues, A. E. Modeling of the methane steam reforming reactor with large pore catalyst. Chem. Eng. Sci. 1992, 47, 2909. (29) Pedernera, M. N.; Pina, J.; Borio, D. O.; Bucala, V. Use of a heterogeneous two dimensional model to improve the primary steam reformer performance. Chem. Eng. J. 2003, 94, 29. (30) Soltan, M. J. S.; Zamaniyan, A. Simulation of terraced wall methane steam reforming reactors. Iran. J. Sci. Technol. 2002, 26, 249. (31) Boersma, J. M. Modeling of NOx emission from natural gas fired gas turbine combustors. Ph.D. Thesis, University of Twente, The Netherlands, 1993. (32) Babu, B. V.; Rakesh, A. Optimal design of an auto-thermal ammonia synthesis reactor. Comput. Chem. Eng. 2005, 29, 1041. (33) Zhu, J.; Zhang, D.; King, K. D. Reforming of CH4 by partial oxidation: Thermodynamic and kinetic analyses. Fuel 2001, 80, 899. (34) Wark, K. W. Thermodynamics, 4th Edition; McGraw-Hill: New York, 1983. (35) Lwin, Y. Chemical equilibrium by Gibbs energy minimization on spreadsheets. Int. J. Eng. Educ. 2000, 16 (4), 335. (36) Chan, S. H.; Ding, O. L.; Hoang, D. L. A thermodynamic view of partial oxidation, steam reforming, and autothermal reforming on methane. Int. J. Green Energy 2004, 1 (2), 265–278. (37) Seo, Y. S.; Shirley, A.; Kolaczkowski, S. T. Evaluation of thermodynamically favorable operating conditions for product of hydrogen in three different reforming technologies. J. Power Sources 2002, 108, 213– 225. (38) Albrecht, B. Reactor modeling and process analysis for partial oxidation of natural gas. Ph.D. Thesis, University of Twente, The Netherlands, 2004. (39) Zamaniayan, A.; Ebrahimi, H.; Soltan, M. J. S. A unified model for top fired methane steam reformers using three-dimensional zonal analysis. Chem. Eng. Process. 2008, 47, 946. (40) Ergun, S. Fluid flow through packed columns. Chem. Eng. Prog. 1952, 48, 89. (41) Harris, S. D. The optimization of reaction rate parameters for chemical kinetic modeling of combustion using genetic algorithms. Comput. Methods Appl. Mech. Eng. 2000, 190, 1065. (42) Konak, A. Multi-objective optimization using genetic algorithms: A tutorial. Reliab. Eng. Syst. Saf. 2006, 91, 992. (43) Zitzler, E.; Deb, K.; Thiele, L. Comparison of multiobjective evolutionary algorithms: empirical results. EVol. Comput. 2000, 8 (2), 173.

ReceiVed for reView February 17, 2009 ReVised manuscript receiVed June 3, 2009 Accepted June 23, 2009 IE900259N