Autothermal Reforming of Gasoline for Fuel Cell Applications: A

Ind. Eng. Chem. Res. 2008, 47, 9437–9446

9437

Autothermal Reforming of Gasoline for Fuel Cell Applications: A Control-Oriented Dynamic Model Yongyou Hu and Donald J. Chmielewski* Department of Chemical and Biological Engineering, Illinois Institute of Technology, Chicago, Illinois 60616

Dennis Papadias Chemical Engineering DiVision, Argonne National Laboratory, 9700 South Cass AVenue, Argonne, Illinois 60439

In this work, we develop a control-oriented, reduced order dynamic model of an autothermal reforming (ATR) reactor. The targeted application is within the on-board fuel-processing unit of a fuel cell vehicle. A previous effort has illustrated that a predictive-type controller may be required to achieve desired performance within this application. The objective of the current effort is to determine the existence of a reduced order model with enough speed and accuracy to meet the online computational demands of a predictive controller. Central to the model development is an approximation of reaction rates that achieve reasonable accuracy near the inlet while preserving the overall energy balance. The resulting scheme converts a partial differential equation model into a set of ordinary differential/algebraic equations and achieves nearly a 4 orders of magnitude improvement in computational speed while preserving most of the nonlinear characteristics of the original system. Such results give clear indication that the hurdle of computational viability can be overcome and opens the door for further development of a predictive controller for the ATR application. 1. Introduction Fuel cells have been touted as a clean alternative to traditional combustion-based power systems. However, in transportation applications the challenges associated with hydrogen storage have motivated the investigation of on-board fuel processing as an alternative. In this work we study just one component of a proposed on-board fuel processor system,1 namely, the autothermal reforming (ATR) reactor. Since the target application is vehicle propulsion, such a reactor will be expected to respond to frequent changes in hydrogen demand as well as numerous start-up/shut-down cycles. As such, throughput control and temperature regulation will be key aspects of reactor operation. Although the use of traditional feedback control is an option, such a scheme is expected to result in limited dynamic performance. Specifically, the highly nonlinear nature of the reactor coupled with severe penalties for exceeding temperature bounds will require conservative tuning within the controller and thus result in sluggish responses during load changes and start-up.2 The alternative is to use a predictive-type controller,3-5 which can incorporate a nonlinear dynamic model and explicitly enforce hard limits on reactor temperature. However, the penalty within the controller for this expected performance improvement is a massive increase in computational burden. This is due to the online optimization routine that resides at the heart of a predictive controller. Qualitatively, this means that hundreds of simulations will likely be required to occur within the sample time of the controller (expected to be on the order of seconds). If we now contrast this computational burden with the simulation time reported in ref 6 (approximately 2 min for a simulation horizon of 180 s), it is easily concluded that the online optimization suggested by predictive control is currently out of reach. In the case of an off-line optimization framework (which includes parameter identification, as well as the look-up table approaches to predictive control7,8), the resulting computational effort is characterized as barely tractable. Thus, the first step * To whom correspondence should be addressed. Tel.: 312-567-3537. Fax: 312-567-8874. E-mail: [email protected].

toward the development of any optimization-based algorithm for the ATR reactor is to construct a sufficiently fast dynamic model. The paper is organized as follows. First we give a brief description of our existing distributed parameter model for the ATR reactor, from ref 6. In section 3, the reduced order model is developed and presented, while simulation results and discussion are given in section 4. 2. Description of the ATR Reactor As discussed above, the ATR reactor is just one component of an overall fuel-processing unit (see Figure 1). The ATR reactor, shown in Figure 2, is composed of a solid catalyst (Rh/ Gd-CeO2) residing within an adiabatic tubular vessel, through which reactant gases are passed at high space velocity (∼50 000/ h). Generally speaking, all reactants are fed in the vapor phase and are preheated to around 250 °C via the recuperating heat exchanger and/or inlet heating rod. In ref 6, a kinetic model for the following mechanism was proposed and validated, the results of which are summarized in Tables 1-3. total oxidation : steam reforming :

CmHn + (m + n/4)O2 f mCO2 + n/2H2O (1) CmHn + mH2O f mCO + (m + n/2)H2 (2)

water gas shift :

CO + H2O T CO2 + H2

(3)

The reactor is usually operated within one of two operational modes. The first is the catalytic partial oxidation (CPOX) mode where the feed consists of fuel and air. Although this mode will quickly heat up the reactor, it suffers from a small H2/CO ratio. The second mode (denoted as the ATR mode) is similar to CPOX but includes steam in the feed, so as to increase the H2/CO ratio. The expected start-up scheme starts with CPOX operation followed by a transition to the ATR mode, where the desired H2/CO ratio will be achieved.

10.1021/ie800209k CCC: $40.75  2008 American Chemical Society Published on Web 11/05/2008

9438 Ind. Eng. Chem. Res., Vol. 47, No. 23, 2008

Figure 1. Flow diagram of the FASTER fuel processor.

Figure 2. Geometrical domain for the adiabatic reforming unit (ref 6). Table 1. Rate Expressions

Table 4. Definitions

reaction

rate equation

1 2 3

r1 ) k1kcyO2 /(k1 + kc) (s) (s) (s) r2 ) k2yCmHnyH2O/(1 + kayCmHn)2 (s) (s) (s) (s) r3 ) k3(yCOyH2O - yCO2yH2 /Ke)

(g)

qr ) (-∆H1)r1 + (-∆H2)r2 + (-∆H3)r3 jqj ) γsqr + qg qg ) βgs(T(g) - T(s)) Rs ) λe,ax/(F cˆp)s(1 - ε) βws ) σwshw/(Fcˆp)s(1 - ε) γs ) 1/(Fcˆp)s(1 - ε) βww ) σwwhw/(Fcˆp)w βgs ) σcshc/(Fcˆp)s(1 - ε) βgg ) Acσcshc/(m ˙ cˆp)g

unit mol/m3/s mol/m3/s mol/m3/s

Table 2. Kinetic Model k1 f ∞ j (m + 0.25n) kc ) σcsF(g)km,O2/M km,O2 ) ShDm,O2/ds Sh ) 0.78Re0.45Sc0.33 Dm,O2 ) 0.196(T(g)/T0)1.75 k2 ) σcsδk02 exp(-E2/RT(s)) ka ) k0a exp(Ea/RT(s)) (s) k3 ) σcsδk03 exp(-E 3/RT ) (s) Ke ) 10((2073/T ) - 2.03)

mol/m3/s mol/m3/s m/s m2/s mol/m3/s mol/m3/s -

Table 3. Kinetic Parameters δ ) 1.4 × 10-5 k02 ) 3 × 1011 k0a ) 105 E2 ) 1.05 × 105 Ea ) 3000 k03 ) 3.4 × 1012 E3 ) 1.3 × 105

M mol/m3/s J/mol J/mol mol/m3/s J/mol

A first-principles model of the reactor should have three independent variables (axial, radial, and time) and should indicate temperatures in the gas, catalyst support, and vessel wall as well as species concentrations in the gas. Toward reducing the computation load of such a simulation, ref 6 employed a number of standard assumptions. The first was a pseudo-steady-state assumption with respect to the gas phase, owing in part to the reactor’s high space velocity but, more importantly, the expected slow temperature response of the solid phase. The second was to assume negligible conduction and diffusion within the gas, again owing to the high space velocity

J/m3/s K/s K/s m2/K/s 1/s m3 K/J 1/s 1/s 1/m

(i.e., dominated by convection). The last standard assumption was with respect to the metal wall of the reactor vessel where conduction is neglected, due to the small annulus of area available for conduction as compared to the large surface area available for heat transfer from the catalyst. Although each of these removed stiffness from the model and improved simulation time, the model still contained three independent variables. Thus, a key contribution of ref 6 was to develop a novel heat transfer correlation that allowed for the elimination of radial direction details, while preserving the heat transfer relationship between the catalyst support and the metal wall. The resulting dynamic model is as follows. (i) Material balances: n

m ˙ (g)

a dwj ) AcMj (νijri), dz i)1

∑

j ) CmHn, CO2, H2O, CO, H2, O2 (4)

(ii) Energy balances: (m ˙ cˆp)g

dT(g) ) Acσcshc(T(s) - T(g)) dz

(5)

Ind. Eng. Chem. Res., Vol. 47, No. 23, 2008 9439

(Fcˆp)s(1 - ε)

(s)

we find the following finite dimensional relationship, indicating the time evolution of the temperature profile coefficients:

2 (s)

∂T ∂T ) λe,ax 2 + σcshc(T(g) - T(s)) + ∂t ∂z σwshw(T(w) - T(s)) + qr (6)

dx ) Ax + qˆ(t) dt

(w)

∂T ) σwwhw(T(s) - T(w)) (7) ∂t where standard no heat flux boundary conditions are applied to eq 6 and qr is defined in Table 4. (Fcˆp)w

(15)

where (s) (s) (w) (w) (w) T x ) [x(s) 1 x2 · · · xN x1 x2 · · · xN ] ,

A)

[

Ass Aws Asw Aww

]

Ass ) diag{-Rsωi - βws}, Aws ) diag{βws}, Asw ) diag{βww}, Aww ) diag{-βww} and 2

3. Model Reduction Although the above model can be considered of reduced order, its solution via finite element methods (FEM) is still too slow for use within an optimization framework. In this section, we illustrate our proposed reduction procedure, specifically tailored to the ATR reactor. The reduction procedure will take place in two steps. The first step focuses on the thermal aspects of the process. By exploiting the linear nature of eqs 5-7, we will arrive at a semianalytic solution procedure, denoted below as reduced order model 1 (ROM1). The second step will focus on the highly nonlinear portion of the model, essentially eq 4. Through application of some reasonable assumptions concerning the reaction rates, we will arrive at a set of analytic functions for the heat generation profiles within the reactor. These functions combined with ROM1 will generate the proposed controloriented model, denoted below as ROM2. Toward clarity of presentation we will stay within the continuous-time domain until the very end of the development, where a discrete-time version will be presented. 3.1. The Solid Material. We begin with a spatial discretization of the solid components via a standard Galerkin procedure.9,10 A rewriting of eqs 6 and 7 (see Table 4) yields ∂2T(s) ∂T(s) ) Rs 2 + βws(T(w) - T(s)) + qj ∂t ∂z

(8)

∂T(w) ) βww(T(s) - T(w)) (9) ∂t where qj contains the heat generated by the reactions as well as heat exchanged with the flowing gas. Next, assume that the catalyst support and vessel wall temperature profiles can be approximated by a linear combination of basis function Fj(z): N

T(s)(z, t) ≈ T(s) N )

∑x

(s) j (t)Fj(z)

(10)

j)1 N

T(w)(z, t) ≈ T(w) N )

∑x

(w) j (t)Fj(z)

(11)

j)1

Given the classic results of Sturm-Liouville theory, the eigenfunctions of eq 8 yield a complete orthonormal sequence over the inner product:

〈Fj, Fi 〉 ) ∫0 Fj(z)Fi(z) dz L

(12)

where Fj(z) ) Hj cos(ωjz), ωj ) (j - 1)π/L, H1 ) 1/L and Hj ) 2/L, j ) 2, 3,.... If we now enforce the conditions 〈FiRN(s)〉 ) 0 and 〈FiRN(w)〉 ) 0, where R(s) N )-

∂T(s) ∂2T(s) N N (s) + Rs 2 + βws(T(w) j N - TN ) + q ∂t ∂z

(13)

∂T(w) N (w) + βww(T(s) N - TN ) ∂t

(14)

R(w) N )-

qˆ(t) ) [〈qj, F1 〉

〈qj, F2 〉 · · · 〈qj, FN 〉

0 · · · 0]T

(16)

3.2. Heat Transfer with the Gas Phase. Next we address the gas-phase energy balance, in an effort to determine the qg contribution to qˆ. We begin by defining a new variable T˜ ( T(g) - T(s) and substitute into eq 5, which yields dT˜ ) -βggT˜ + g(z, t) dz

(17)

N where g(z, t) ) -dT(s)/dz ≈ ∑j)1 xj(s)ωjHj sin (ωjz). Then exploitation of the linear nature of eq 17 yields

N

T˜(z, t) ) e-βggzT˜(0, t) +

∑ x k (z) j j

(18)

j)1

where kj(z) ) Hj{ωjβgg sin(ωjz) - ωj2 cos(ωjz) + ωj2e-βggz}/ 2 (βgg + ωj2). Then, insertion into eq 16 yields qˆg(t) ) βgs[〈T˜, F1 〉 · · · 〈T˜, FN 〉 ]T ) M1x + M2T˜(0, t) (19) where M1, M2 are matrix functions of βgg (see Appendix A for detailed expressions). Through eq 19 we find that the heat transfer term between the gas phase and the catalysts support is a linear function of the inlet gas temperature as well as the coefficients x(s) j of eq 10. However, we also find this heat transfer term to be a nonlinear function of the mass flow rate of the gas (see βgg of Table 4 as well as Appendix A). Since mass flow is expected to be among the manipulated variables of the ATR controller, this will be a nonlinear input to the system. 3.3. Reaction-Generated Heat. Toward calculating the reaction-generated heat, it is clear that we must first solve the material balances indicated by eq 4. However, due to the heterogeneous nature of these reactions, we need not simultaneously solve the gas-phase energy balance. Rather, it is more appropriate to simply apply the instantaneous solid phase temperature profile, indicated by eq 10, to the kinetic expressions found in eq 4. Once we have the resulting concentration profiles, it is a simple matter to complete the model by determining the profile of heat evolved from the three reactions and apply to eq 16. As an initial model reduction step, denoted ROM1, we solve the six spatial ODEs of eq 4 using a standard ODE solver (ode15 from Matlab). As seen in Figures 3 and 4, ROM1 compares quite favorably with the original CFD simulations of ref 6. An additional note of the ROM1 simulations is that we have assumed zero mass transfer resistance between the bulk and (g) surface concentrations (y(s) j ) yj for all species except oxygen, since mass transfer is the basis of r1, see Table 1). In sum, we conclude that the Galerkin approximation of the previous subsections (along with the negligible mass transfer resistance assumption) provides little degradation with regard to model accuracy. Unfortunately, these assumptions also provide little improvement regarding simulation speed. Specifically, for the 180 s simulation horizon of Figure 3 (time starts from 20 s), ROM1 required 30 s of CPU time. (It should be noted that the


Figure 3. Operating conditions (left) and temperature response (right).

Figure 4. Concentration and temperature profiles at steady state (t ) 135 s of Figure 3 for concentrations and t ) 155 s for temperature).

FEM used an adaptive step-size routine with regard to the evolution of time. In contrast, ROM1 assumed a uniform sample time over the entire simulation horizon. Although adaptive step sizing could have been applied to ROM1 to achieve faster simulations, such a scheme does not fit within the framework of typical predictive controllers.) In the remainder of this section, we indicate a procedure aimed at eliminating the need for a numeric-based solution to eq 4. This purely analytic model will be denoted as ROM2. We start with reaction 1. A quick inspection of the oxygen fraction profile of Figure 4 indicates that this reaction goes to completion within the first millimeter of the reactor. It is also noted that, due to the simplicity of the rate term, an analytic expression for rate 1 can be logically determined. r1(z) ) r◦1 e-a1z

(20)

where a1 ) Acσckm,O2F(g)/m ˙ , r1° ) σckm,O2F(g)wO° 2/(m + n/4)MO2, ° and wO2 is the inlet mass fraction of oxygen. Then we use the analytic expression 20 to calculate the oxidation portion of qr: qr1(z) ) (-∆H1)r1(z). Since reaction 1 is virtually complete in the first millimeter of the reactor, it is also reasonable to assume that the species associated with this reaction are completely converted at the inlet plane of the reactor. Said another way, ROM2 assumes that the material balance portion of combustion

reaction will complete before the reforming and WGS reactions start. As such the inlet concentrations for the reforming and WGS reaction portion of ROM2 are as follows: wC1 mHn ) wC◦ mHn 1 ◦ wCO ) wCO + 2 2

wH1 2O ) wH◦ 2O +

MCmHn MO2(m + n/4) mMCO2

MO2(m + n/4) 0.5nMH2O MO2(m + n/4)

wO◦ 2

(21)

wO◦ 2

(22)

wO◦ 2

(23)

The inlet concentrations of hydrogen and CO are unchanged 1 ° (wH1 2 ) wH° 2 and wCO ) wCO ) and wO1 2 ) 0. Next we turn to reaction 2. As shown in Figure 4, the fuel and hydrogen concentration profiles show characteristics similar to that of an exponential function. Thus, we approximate the fuel and hydrogen concentration profiles by the following: wCmHn(z) ) w˜CmHn(1 - e-a2z) + wC1 mHn

(24)

wH2(z) ) w˜H2(1 - e-a2z) + wH1 2

(25)


where w˜CmHn ) wCmHn - wCmHn and w˜H2 ) wH2 - wH2 (the exit mass fractions deviated from the postcombustion fractions). As shown in Figure 5 (CPOX mode: fuel, 40 g/min; air, 90 slpm; steam, 0 g/min), wCe mHn and wHe 2 are not constant but respond in a way similar to a first-order system. Thus, we approximate w˜CmHn and w˜H2 as e

∞

1

e

1

w˜CmHn(t) ) w ˜ C∞mHn(1 - e-t/τ)

(26)

˜ H∞2(1 - e-t/τ) w˜H2(t) ) w

(27)

w˜H∞2

where w˜CmHn and are the steady-state versions of w˜CmHn and w˜H2 (see Appendix B for details). The time constant, τ, shows a linear relationship with reactor residence time, tr, as indicated by Figure 6. In this work, we assume τ ) b1tr + b2, where b1 ) 469 and b2 ) -13.45. However, during nominal operation (i.e., after cold start-up), it was observed that the exit concentration transients were very fast and could be ignored. Substitution of eq 24 into the fuel balance portion of eq 4, gives the following analytic expression for reaction rate 2: r2(z) ) -w˜CmHna2e-a2z/R

completes at the inlet plane, the reforming rate at the inlet is calculated as r2◦ ) K2wC1 mHnwH1 2O j2

(s)

a2 )

1

r2◦R -w˜CmHn

(30)

The analytic expression 28 can now be rewritten as r2(z) ) r2◦e-a2z

(31)

and used to calculate the reforming portion of qr:qr2(z) ) (-∆H2)r2(z). Turning to the WGS reaction, substitution of eqs 25 and 31 into the hydrogen balance portion of eq 4 gives the following analytic expression for reaction rate 3: r3(z) ) r3◦e-a2z

(28)

where R ) AcMCmHn/m ˙ . Assuming the combustion reaction

(29)

j wCmHn/MCmHn)2 where K2 ) k2(T0 )M /MCmHnMH2O(1 + ka(T0 )M (s) and T0 is the temperature of the solid at the inlet. Combining eqs 28 and 29 yields the following expression for a2: (s)

°

∞

∞

(32) °

where r3 ) -(MCmHnw˜H2/MH2w˜CmHn + m + 0.5n)r2. The WGS portion of qr is then given as qr3(z) ) (-∆H3)r3(z). Now that we have procedures for calculating each of the three heat generation terms, we can combine them and rewrite qr in the compact form: qr(z) ) c1e-a1z + c2e-a2z

(33)

where c1 and c2 are easily determined from the above development. We can now proceed to incorporate qr into eq 16.

〈c1e-a z, Fi 〉 ) 1

〈c2e-a z, Fi 〉 ) 2

Figure 5. Dynamics of concentrations at exit.

Figure 6. Relation between τ and tr.

Hic1a1 a12 + ωi2 Hic2a2 a22 + ωi2

{1 - (-1)i-1e-a1L}

(34)

{1 - (-1)i-1e-a2L}

(35)

Within the simulations we assumed that e-a1L ) 0 and e-a2L ) 0. The reasonableness of these assumptions is due to the overdesign of the reactor used in this study. In the case of a smaller reactor it is likely that the second term of eq 35 would need to be retained. It should be also noted that the constants a1 and a2 are functions of the inlet conditions and solid temperature and, thus, should be recalculated at every sample of the simulation. Let us now summarize the above development. Basically, two guiding principles were used to help develop the approximations. The first was that the overall conversion rates should be the same as the original model. This has the important benefit of preserving the overall energy balance on the reactor. Of

Figure 7. Operating conditions (left) and temperature response (right) for CPOX mode.


Figure 8. Concentration and temperature profiles at steady state for a CPOX mode.

original kinetic model evaluated at z ) 0. In the case of r2° , we again employ the original kinetics, but rather than evaluate at the inlet conditions we employ the more reasonable postcombustion conditions defined by eqs 21-23. As we will see in the next section, this procedure of approximating the inlet rates while preserving the overall conversion rate provides a good balance between model accuracy and computational simplicity. 4. Model Validation

Figure 9. Temperature response for transition mode (CPOX to ATR).

particular concern was the temperature-dependent equilibrium constant of the WGS reaction, evaluated at the reactor exit temperature (see Appendix B for details). From a mathematical perspective, this assumption constrained the area beneath each function. Toward approximating the shape of each function, we employed the second guidesinlet rates should be as close as possible to the original model. For reaction 1, this was clearly observed, as indicated by the fact that r1° is essentially the


Before proceeding to model validation, we will need to make two final decisions. The first concerns the number of eigenfunctions used to approximate the axial temperature profile within the reactor (i.e., N of eqs 10 and 11). Although a large N will improve accuracy it also increases computational burden. Toward this effort, we moved the no heat flux boundary condition at the exit from 31 to 19 mm. This allowed us to reduce the number of eigenfunctions and maintain sufficient profile fidelity at the inlet region. The second decision regards the sample period of the discrete-time version of the model. Since the intended application is use within an optimization framework, a uniformly sampled discrete-time model is required. Again, appropriate selection of this parameter will significantly reduce computational effort. After a bit of trial and error, we finally settled on N ) 20 and ∆t ) 0.5 s. These values will be


Figure 11. Concentration and temperature profiles at steady state (t ) 135 s of Figure 3 for concentrations and t ) 155 s for temperature).


Figure 13. Concentration and temperature profiles at steady state for an ATR mode.

used in the remainder of the analysis. (It should be noted that the ROM1 simulations of Figures 3-6 and those to follow used N ) 80, ∆t ) 0.25 s, and a reactor length of 31 mm.)

4.1. Start-Up Scenarios. We begin by comparing the ROM1 simulations with those of ROM2 for CPOX mode (fuel, 40 g/min; air, 90 slpm; steam, 0 g/min). On the basis of Figure 7,


we see that the temperature dynamics at 0.8, 7, and 19 mm compare quite well. The mass fraction profiles, shown in Figure 8, also match very well despite some discrepancy at the inlet, due to our assumption of complete combustion at the inlet plane. In the steady-state temperature profile, an error of about 30 °C is observed in the 2-7 mm interval, which we attribute the slight errors resulting from our assumptions regarding the reforming and WGS reaction rates, highlighted by the slight overconversion of fuel observed in Figure 8. The transition from CPOX mode to ATR mode is an important scenario for reactor start-up. Figure 9 shows simulations for a transition from CPOX mode (fuel, 40 g/min; air, 90 slpm; steam, 0 g/min) to ATR mode (fuel, 40 g/min; air, 165 slpm; steam, 90 g/min). Clearly, ROM2 captures the temperature dynamics quite well with less than a 15 °C error as compared with ROM1. In our next comparison (Figure 10) we applied the inlet flows of Figure 3. From 50 to 135 s, the solid temperature profiles show significant discrepancy (30 and 10 °C for the 7 and 19 mm locations, respectively). However, once the air flow increases from 130 to 140 dm3/min the trajectories seem to match quite well. This sudden change in accuracy seems to stem from the slight inaccuracies resulting from the reforming model of ROM2 (see Figure 11), similar to our analysis of Figures 7 and 8. 4.2. Nominal Operation. Once the reactor has been startedup (i.e., after the transition from CPOX to ATR mode), the dominate disturbance will be load (or throughput) changes while operating in ATR mode. Figure 12 is a simulation of start-up under ATR operating conditions (fuel, 40 g/min; air, 165 slpm; steam, 90 g/min). Figure 13 shows the steady-state profiles for this case, and will serve as the baseline for the step responses to follow. It is interesting to note that the significant discrepancy observed in the steady-state temperature profile (around 40 °C) at some locations due to the small number of eigenfunction we used in ROM2. If the number of eignfunctions is increased from 20 to 40, the performance will improve significantly, but the computational cost will also double. Figure 14 shows comparisons of temperature responses at 7 mm to changes in air flow, steam flow, and inlet temperature. We see that ROM2 captures the responses quite well, about a 3 °C difference as compared with the CFD simulations. 4.3. Computational Effort. As shown in Table 5, the computational effort has been significantly reduced. For a time horizon of 180 s, approximately 2 min are required for the FEM calculations, about 30 s for ROM1, and only 0.04 s for ROM2. It should be noted that during start-up and load changes the time constant of the system is about 30 s. If we were to use this as the simulation horizon for a predictive controller, then the computational cost for the ROM2 would be less than 0.007 s. Assuming the sample time of the controller is 0.5 s this would allow for approximately 70 simulations per control move. 5. Conclusion The objective of this work was to show that a reduced order dynamic model of the ATR reactor with accuracy and speed sufficient for predictive control applications could be achieved. The presented results clearly indicate that such a model is possible and gives credence to the further investigation and the genuine possibility of applying predictive control to this dynamically fast process. Future work will focus on developing such controllers under the now realistic assumption that computational tractability can be achieved.

Figure 14. Step response (top, airflow 165-180 slpm; middle, steam flow rate 90-110 g/min; bottom, inlet temperature 250-275 °C). Table 5. Comparison of Computational Effort

a

numeric method

FEM

ROM1a

ROM2a

CPU time, s

110

30

0.04

Matlab with Core 2 Duo CPU 2.33G Hz and 1G DDR RAM.

Acknowledgment This manuscript was created by the University of Chicago, as an Operator of Argonne National Laboratory, under Contract No. W-31-109-ENG-38, in association with the U.S. Department of Energy. This work was supported by the U.S. Department of Energy Hydrogen, Fuel Cells, and Infrastructure Technologies Division. Special thanks to Nancy Garland, Valri Lightner, Patrick Davis, and Steve Chalk for supporting this study. The authors also thank the participants of the Feasibility of Acceptable Start-Time Experimental Reformer (FASTER) project. This was a collaborative project with contributions from Los Alamos National Laboratory, Oak Ridge National Laboratory, Pacific Northwest National Laboratory, and a number of universities and private organizations providing catalysts, heat exchangers, and hardware build support. The authors gratefully acknowledge


the team at Argonne National Laboratory who participated in the assembly, tests, hardware design, and analysis of the fuel processor: Steven Calderone, Tod Harvey, Sheldon Lee, Hsui-Kai Liao, Daniel Applegate, Vincent Novick, Rajesh Ahluwalia, Steven Lottes, and Shabbir Ahmed.

fractions after the reverse WGS completes are as calculated to be wC2 mHn ) wC1 mHn

(B.2)

2 1 1 ) wCO + wCO MCO/MCO2 wCO 2

(B.3)

2 )0 wCO 2

(B.4)

1 MH2O/MCO2 wH2 2O ) wH1 2O + wCO 2

(B.5)

1 MH2/MCO2 wH2 2 ) wH1 2 - wCO 2

(B.6)

Appendix A: Matrix Functions in the Gas-Phase Model Expanding eq 18 we find T˜(z, t) ) e-βggzT˜(0, t) + N Hj{ωjβgg sin(ωjz) - ωj2 cos(ωjz) + ωj2e-βggz} (A.1) x(s) j βgg2 + ωj2 j)1

∑

Then insertion into eq 16 yields -β z 〈qg, Fi 〉 ) βgs〈T˜, Fi 〉 ) βgsHi∫0 cos(ωiz)e gg dzT˜(0, t) + L

N

βgsHi

∑β j)1

xjHjωj 2

gg

∫

+ ωj

L

2

0

Now that we have the available fuel from eq B.2 and the available steam from eq B.5, we can determine if steam will be completely consumed within an infinite length reactor, given by the following condition: wH2 2O < mwC2 mHnMH2O/MCmHn

cos(ωiz){βgg sin(ωjz) - ωj cos(ωjz) + ωje-βggz} dz (A.2)

Evaluating these integrals for i ) 1, 2,..., N and assuming e-βggL ≈ 0, we find M1 ) [βgs(βggFD - ED + βggDGD) ON×N ] M2 ) βgsβggDP where D ) diag

{

βgg

2

1 + ωi2

}

or equivalently in terms of inlet conditions wH◦ 2O/MH2O + 2w◦ O /MO2 < mwC◦ mHn/MCmHn

In the alterative case fuel will be completely consumed. If B.7 holds, then the resulting mass fractions after the reforming reaction completes are given as wC3 mHn ) wC2 mHn - wH2 2OMCmHn/mMH2O

(B.8)

3 2 ) wCO + wH2 2OMCO/MH2O wCO

(B.9)

3 2 ) wCO wCO 2 2

(B.10)

wH3 2O ) 0

(B.11)

wH3 2 ) wH2 2 + (m + 0.5n)wH2 2OMH2/mMH2O

(B.12)

P ) [H1 · · · HN ]T E)

L diag {Hi2ωi} 2

F ) diag {Hi} · Q · diag {Hiωi}

Then w˜C∞mHn and w˜H∞2 are easily determined as

cos{(ωj - ωi)L} Q(i, j) ) 2 2 2(ωj - ωi) ω -ω ωj

j

w˜C∞mHn ) wC3 mHn - wC1 mHn

i

cos{(ωj + ωi)L} , 2(ωj + ωi)

i, j ) 1, 2, ... , N )-

G ) PP diag {ωi } T

(B.7)

2

2

It is additionally noted that P, E, Q, G, and βgs do not change during the course of the simulation, whereas βgg and D are functions of mass flow rate and need to be recalculated at each time step.

m(4m + n)MO2

w˜H∞2 ) wH3 2 - wH1 2 )

nMH2 mMO2

ww◦ O 2

ww◦ O +

MCmHn mMH2O

wH◦ 2O

(m + 0.5n)MH2

2

mMH2O

(B.13)

ww◦ O

2

(B.14) Returning to the alternative of B.7, where fuel is completely consumed, we assume the reforming reaction will start first:

∞ Appendix B: Calculation of w˜C∞mHn and w˜H 2

A reasonable assumption is that there is still a great deal of fuel remaining after the combustion reaction is complete. wCmHn◦/MCmHn . wO2◦/(m + 0.25n)MO2

(4m + 2n)MCmHn

(B.1)

It is noted that the reforming reaction will continue until either the fuel or steam reaches zero. The available steam for the reforming reaction has three sources: the feed, steam produced by combustion, and steam produced by reverse WGS. For the purpose of this calculation, let us assume the reverse WGS goes to completion in an effort to calculate the maximum amount of H2O available (the resulting hydrogen fraction may become negative, but on the basis of B.1, it will become positive again once the reforming reaction is applied). The resulting mass

wC2 mHn ) wC1 mHn

(B.15)

2 1 ) wCO + mwC1 mHnMCO/MCmHn wCO

(B.16)

2 1 ) wCO wCO 2 2

(B.17)

wH2 2O ) wH1 2O - mwC1 mHnMH2O/MCmHn

(B.18)

wH2 2 ) wH1 2 + (m + 0.5n)wC1 mHnMH2/MCmHn

(B.19)

Next we activate the WGS reaction and assume it reaches equilibrium at the exit. If the conversion of this reaction is ξ, the following will hold:

9446 Ind. Eng. Chem. Res., Vol. 47, No. 23, 2008 2 (wCO /MCO - ξ)(wH2 2O/MH2O - ξ)

1 ) (s) 2 2 (wCO /M Ke(T + ξ)(w /M + ξ) CO2 H2 H2 e ) 2

(B.20)

Then, the energy balance based on eqs 5-7 will lead to (g) cˆpg(T(g) e - Tin ) )

wO2◦(∆H2 - ∆H1) (m + n/4)MO2

+

wCmHn◦(-∆H2)

Greek Symbols

MCmHn + ξ(-∆H3)

or (g) T(s) e ≈ Te )

{

wO2◦(∆H2 - ∆H1) (m + n/4)MO2

+

wCmHn◦(-∆H2) MCmHn

}

+ ξ(-∆H3) /cˆpg + T(g) in (B.21) Clearly, ξ can be determined by substituting eq B.21 into eq B.20 and solving the nonlinear equation. Then w˜C∞mHn, w˜H∞2 is given by w˜C∞mHn ) -wC1 mHn

(B.22)

w˜H∞2 ) (m + 0.5n)wC1 mHnMH2/MCmHn + ξMH2

(B.23)

However, it is noted that solving the resulting nonlinear equation will consume significant time due to the iterations involved. Thus, we do not solve this equation directly but allow it to evolve iteratively during the simulation based on an initial guess for ξ (e.g., ξ ) 0 at t ) 0). In comparison with full iteration at every time step, this simpler and much faster method resulted in very little degradation in accuracy. Nomenclature Ac ) cross-sectional area of catalyst bed, m2 cp ) heat capacity, kJ/(kg · K) ds ) mesh diameter, m Dm ) gas diffusion, m2/s E ) activation energy, J/(mol · K) hc ) heat transfer coefficient in catalyst bed, W/(m2 · K) hw ) heat transfer coefficient to reactor wall, W/(m2 · K) k0 ) frequency factor, mol/m3/s. k0a ) frequency factor for absorption kc ) mass transfer coefficient, m/s L ) reactor length, m m ) number of carbon atom m ˙ ) mass flow rate, kg/s M ) molecular mass, kg/mol j ) average molecular mass, kg/mol M n ) number of hydrogen atom r ) reaction rate, mol/m3/s R ) specific gas constant, J/(mol · K) t ) time, s T ) temperature, K

V ) volume, m3 w ) mass fraction x ) time-dependent coefficient of eigenfunction, K · m0.5 y ) mole fraction z ) axial distance, m

∆H ) reaction heat, J/mol δ ) washcoat thickness, m ε ) porosity ω ) frequency, 1/m λ ) thermal conductivity, W/(m · K) νij ) stoichiometric coefficient σcs ) surface area of catalyst per volume of support, m-1 σws ) surface area of wall per volume of support, m-1 σww ) surface area of wall per volume of wall, m-1 F ) density, kg/m3 ξ ) conversion of WGS, mol/kg Subscripts and Superscripts g ) gas phase j ) the jth component s ) solid phase w ) reactor metal wall ° ) initial

Literature Cited (1) Ahmed, S.; Ahluwalia, R.; Lee, S. H. D.; Lottes, S. A Gasoline Fuel Processor Designed to Study Quick-Start Performance. J. Power Sources 2006, 9, 214. (2) Hu, Y.; Chmielewski, D. J.; Papadias, D. Autothermal Reforming of Gasoline for Fuel Cell Applications: Controller Design and Analysis. J. Power Sources 2008, 182, 298. (3) Qin, J. S.; Badgwell, T. A. An Overview of Nonlinear Model Predictive Control Applications. Presented at the International Symposium on Nonlinear Model Predictive Control, Ascona, Switzerland, 1998. (4) Manousiouthakis, V.; Chmielewski, D. J. On Constrained InfiniteTime Nonlinear Optimal Control. Chem. Eng. Sci. 2002, 57 (1), 105–114. (5) Tenny, M. J.; Wright, S. J.; Rawlings, J. B. Nonlinear Model Predictive Control via Feasibility-Perturbed Sequential Quadratic Programming. Comput. Optim. Appl. 2004, 28 (1), 87–121. (6) Papadias, D.; Lee, S. H. D.; Chmielewski, D. J. Autothermal Reforming of Gasoline for Fuel Cell Applications: A Transient Reactor Model. Ind. Eng. Chem. Res. 2006, 45, 5841–5858. (7) Pistikopoulos, E. N.; Dua, V.; Bozinis, N. A.; Bemporad, A.; Morari, M. On-Line Optimization via Off-Line Parametric Optimization Tools. Comput. Chem. Eng. 2002, 26, 175–185. (8) Sakizlis, V.; Kakalis, M. P.; Dua, V.; Perkins, J. D.; Pistikopoulos, E. N. Design of Robust Model-Based Controllers via Parametric Programming. Automatica 2004, 40 (2), 189–201. (9) Ray, W. H. AdVanced Process Control; McGraw-Hill: New York, 1981. (10) Fletcher, C. A. J. Computational Galerkin Methods; SpringerVerlag: New York, 1984.

ReceiVed for reView February 5, 2008 ReVised manuscript receiVed July 16, 2008 Accepted July 28, 2008 IE800209K

Autothermal Reforming of Gasoline for Fuel Cell Applications: A

Recommend Documents