Novel Reactor Concepts for Thermally Efficient Methane Steam

4796

Ind. Eng. Chem. Res. 2004, 43, 4796-4808

Novel Reactor Concepts for Thermally Efficient Methane Steam Reforming: Modeling and Simulation G. Kolios,* A. Gritsch, B. Glo1 ckler, and G. Sorescu Institute for Chemical Process Engineering (ICVT), Bo¨ blinger Strasse 72, 70199 Stuttgart, Germany

J. Frauhammer Division of Gasoline Systems, Robert-Bosch GmbH, Postfach 300240, 70442 Stuttgart, Germany

Simulation of chemical reactors is often considered as a routine task, tacklable with standard models and standard simulators. However, the analysis, design, and optimization of novel integrated reactor concepts require specialized approaches. A countercurrent wall reactor concept for the autothermal coupling of high-temperature endothermic and exothermic reactions is discussed as a case study in order to demonstrate the importance of individual solutions. The custom of insisting on a given model formulation regardless of its numerical implications is addressed in particular. It is demonstrated how simple modifications and heuristics improve the solvability of a problem. The role of model reduction in capturing the multiplicity of solutions is discussed. The comprehensive picture obtained in this way provides a safe basis for deciding on the practical feasibility of the investigated process. It finally generates guidelines for a specific process optimization. 1. Introduction Multifunctional autothermal reactor concepts have been proposed for the thermal coupling of methane steam reforming with an auxiliary combustion reaction to an efficient heat-integrated process. These concepts aim at minimizing fuel consumption and waste heat production in stand-alone hydrogen production, which gains increasing interest for innovative power generation systems, e.g., fuel cells. They feature countercurrent heat exchange of two process streams as a prerequisite for efficient heat recovery.1 Figure 1 shows a sketch of the countercurrent reactor together with the desired operating conditions. The major characteristics of the targeted solution can be summarized in the following points: (i) equal heat capacity fluxes in the heat-exchanger zones and a minimal heat surplus as a driving force for integrated heat recovery and (ii) axial overlapping of reforming and combustion and efficient heat transfer between the two compartments as a prerequisite for a moderate maximum temperature close to the temperature required for the desired equilibrium conversion of reforming. A priori considerations led to a countercurrent wall reactor as the obvious implementation of this concept. Depositing the reforming and the combustion catalyst on opposite surfaces of the separating wall promotes the direct thermal coupling of heat sources and heat sinks. Narrow channels and the resulting high specific heattransfer area facilitate efficient heat recovery from the product to the feed streams.2 Details on the properties of the countercurrent reactor can be found in refs 3-5. Here we are focusing on the software developed in order to gain a better process understanding and to optimize the process design. The major goal is to clarify * To whom correspondence should be addressed. Tel.: +49 (0)711 6412214. Fax: +49 (0)711 6412242. E-mail: [email protected].

Figure 1. Operation principle of a countercurrent reactor for coupling endothermic and exothermic reactions. The diagram shows schematical temperature and conversion profiles as expected from an optimized design.

whether operating conditions with the prespecified properties lie within the attainable region of the reactor concept. Mathematical modeling of the countercurrent reactor leads to a set of parabolic partial differential equations in time t and in the axial reactor coordinate z, which is quite common in chemical engineering.5 Two major trends prevail nowadays in the simulation of these processes. On the one hand, there are integrated packages for computer-aided process engineering.6-8 These packages support the user during modeling, process analysis, design, and optimization. They provide special simulation languages and standard interfaces that enable users to formulate the problem in a convenient way. On the other hand, CFD software is widely applied for the simulation of chemical reactors.9,10 The idea behind these concepts is to release the user from routine tasks and to open process simulation to engineers, who are not familiar with the details of numerical analysis. The standard procedure for numerical treatment of the model equations is based on the approximation of the spatial derivatives by finite differences or finite elements. In this way the model equations are converted into large sets of algebraic equations for stationary problems or sets of differential algebraic equations in the case of dynamic modeling. These tools passed

10.1021/ie034158e CCC: $27.50 © 2004 American Chemical Society Published on Web 07/14/2004

Ind. Eng. Chem. Res., Vol. 43, No. 16, 2004 4797

through a tremendous progress during the past decade regarding their numerical performance as well as their handling. However, the implemented numerical algorithms are in many cases proprietary and optimized primarily with respect to robustness in order to cover a broad spectrum of applications. Users have generally a restricted access to the numerical kernel of the packages: Only selection among different preinstalled solvers and manipulation of a few optional control parameters of the numerical algorithms are permitted. Therefore, these tools are very useful for applications where modeling and the numerical procedure do not interfere with each other. Theoretical analysis of the countercurrent wall reactor required a specific numerical treatment. The major reason was a lack of empirical knowledge about the process behavior: It was impossible to estimate numerical hurdles and to provide adequate remedies in advance because even qualitative features of the expected solution were unclear. It was indispensable to develop the numerical procedure and the model in parallel. Furthermore, when different aspects of the system are studied, different models have been introduced that are not compatible with standard numerical concepts. These reasons led to the development of a library of in-house algorithms. This approach is elucidated in the present paper, which is organized as follows: Numerical aspects of the countercurrent reactor model and implementation details are discussed in section 2. Selected examples highlight the importance of model tuning for the sake of numerical efficiency. In section 3, we discuss a simplified model derived for systematic parametric analysis that yields the complete set of possible steady-state solutions. Finally, in section 4, we introduce a modified concept with optimized properties. Fully transparent numerical algorithms and a package for interactive simulation and online visualization are essential components of the developed software. The most sophisticated tool developed in this context is PDEXPACK.11,12 PDEXPACK has been designed for solving systems of partial differential equations of parabolic type in chemical engineering applications. It is a derivate of the algorithm introduced by Eigenberger and Butt.13 Many of the original ideas, e.g., adaptivity in time and space, error estimation, and the user interface design, can be found in PDEXPACK in an advanced and refined form. 2. Modification and Simplification of the Model The reference model for simulation of the countercurrent reactor is described elsewhere5 together with the main assumptions. In the following, we discuss the modifications and simplifications that have been introduced in order to improve its solvability. The manipulations apply mainly to nonlinear terms, i.e., the reaction rate expressions, the correlations for caloric constants, and heat-transfer coefficients. 2.1. Reaction Rate Expressions. The numerical complexity of the problem is mainly determined by the nonlinearity of the reaction rate expressions. This is due to the Arrhenius term in the rate constant and due to the nonlinear order of the reaction. This issue is discussed in detail because it appears quite frequently in chemical reaction modeling: Many reaction rate expressions published in the literature are given in an

inappropriate form for numerical simulation. The kinetics of the reforming reaction over a Ni catalyst adopted from ref 14 illustrate the difficulty. The reaction rate of the main reaction CH4 + H2O h CO + 3H2 is given by

r1 )

(

kCO(T) pCH4pH2O -

)

pCOpH23 K1(T)

1 ads 2 ads 2 ads ads [KH p + (1 + K p + K p )p + K p ] p H2 CO CO CH4 CH4 H2 H 2O H 2O x H 2 2

This expression obviously tends to infinity as pH2 approaches zero, which is the case at zero conversion. The problem has been overcome by linerarization. A threshmin old pH is specified, below which the linearized form of 2 the rate expression holds: min min pH2 < pH ri ) ri(pH )+ 2 2

∂ri ∂pH2

|

min (pH2 - pH ) 2

min pH 2

An approximation through numerical differentiation is preferred to analytical evaluation because ∂ri/∂pH2|pHmin2 leads to a very complex expression. However, it is sufficient to linearize not the entire rate expression but only the factors that generate the singularity. Applying this approach to the above rate expression, we obtain min w r1 ) pH2 < pH 2

(

kCO(T) pCH4pH2O -

)

pCOpH23 K1(T)

ads ads 2 ads [KH pH22 + (1 + Kads CO pCO + KCH4pCH4)pH2 + KH2OpH2O] 2

(

1-

×

)

min pH2 - p H 2 min 2 pH 2

xp

min H2

min A reasonable threshold value of pH2 is pH ) 2 norm 0.5rtolpH2 . In this way, the inaccuracy of the manipulated kinetic expression is restricted within the tolerance limit of numerical integration while avoiding round-off errors. A similar pretreatment is necessary for the kinetics of the homogeneous combustion of methane. The rate expression is adopted from ref 15

r ) k(T) cCH4-0.3 cO21.3

(1)

This expression is numerically unfavorable in several aspects: (i) It is not defined for cCH4 e 0 or cO2 < 0. However, these conditions may occur during the numerical procedure, particularly during ignition of the homogeneous combustion. (ii) The reaction rate tends to infinity, which is physically infeasible, as cCH4 approaches zero. These drawbacks can be eliminated through a modified expression that holds in the shaded area in the cCH4-cO2 plane (Figure 2). The modified expression must continue smoothly to the original expression at the boundary of the two half-planes. In addition, the reaction rate must be zero if the concentration of one of the

4798 Ind. Eng. Chem. Res., Vol. 43, No. 16, 2004

Figure 2. cCH4-cO2 plane divided into two sections where different rate expressions hold: (1) original rate expression according to eq 1; (2) modified rate expression according to eq 2.

reactants is zero. The following polynomial expressions satisfy these conditions:

r ) k(T) fsign f1 f2

{ {

(2)

fsign ) sign[sign(cCH4) + sign(cO2) + 1]

f1 )

f2 )

eCH4 cCH 4

cCH4 g cCH4,min

eCH4-2 2 -(1 - eCH4)cCH c + 4,min CH4 eCH4-1 (2 - eCH4)cCH4,mincCH4 eCH4-1 (2 - eCH4)cCH c 4,min CH4

0 < cCH4 < cCH4,min cCH4 e 0

(eCH4 ) -0.3)

eO2 cO 2

cO2 g cCO2,min

eO2-2 2 -(1 - eO2)cO c + 2,min O2 eO2-1 (2 - eO2)cO2,mincO2 eO2-1 (2 - eO2)cO c 2,min O2

0 < cO2 < cO2,min

it would be attractive to use constant values for a coarse analysis of the basic countercurrent reactor properties. In this way the integration procedure becomes very robust and the required CPU time is reduced about one order of magnitude as compared with the detailed correlation. However, proper selection of approximate values for cp is not a trivial task. An obvious choice without further knowledge about the process would be the values at inlet conditions (right). However, this leads to significant deviations from the reference solution. A better approximation is achieved by exploiting some empirical knowledge: The values of cp are computed for a hypothetical mixture at an intermediate conversion and temperature level (middle). Nevertheless, it is important to check the effect of simplifications compared to the detailed model in order to avoid artifacts and misinterpretation of the simulation results. 2.3. Heat Radiation. The reactor model includes radiative heat transfer in the energy balance of the solid phase. This is motivated by the high temperature level and the strong axial gradients in the temperature profiles. Radiation may affect the heat transfer, particularly in the straight, empty channels of an ideal wall reactor. Radiation acts in a qualitatively different way than other heat-transfer mechanisms in the sense that it is a long-range effect and not coupled to local temperature differences. Heat received at a certain spot depends on the temperature of all surfaces visible from it and vice versa. The effective heat flux per unit area at the axial position z is given by the following expression:

cO2 e 0

(eO2 ) -1.3)

The modified reaction rate expression is crucial for achieving convergence of the numerical procedure. The described extension is generally applicable to power law expressions: It ensures continuity and smoothness of the rate expression, zero passage at zero concentration, and a meaningful extension in the range of negative concentrations in order to avoid spurious solutions. 2.2. State-Dependent Coefficients. Temperatureand composition-dependent thermodynamic, caloric, and transport parameters contribute to the nonlinearity of the model and potentially obstruct the numerical solution. It is therefore important to examine whether this degree of complexity is crucial for the reactor behavior. The effect of variable heat capacities is discussed as an example because the heat capacity fluxes of the process streams are known to have a significant impact on countercurrent heat exchange.1 A given configuration has been simulated with constant and state-dependent specific heat capacities of the reaction mixture. Figure 3 shows the effect of various correlations for the specific heat capacity cp on the steady-state solution. Using the detailed correlation (left), the required CPU time is unacceptable in the context of extended parametric studies. This is partly due to the numerical effort for evaluating the correlations but also due to very short time steps taken by the DAE solver for the integration in time. According to our experience, the performance of PDEX can be significantly affected by state-dependent coefficients of the differential terms in space. Therefore,

In the discretized system, the solid-phase energy balance of each node of the spatial grid depends on state variables of all other nodes. Hence, the system matrix no longer has a banded structure as required in the context of the PDEX algorithm. All entries outside a band that covers the coupling of three adjacent nodes of the grid are omitted. The incomplete system matrix generally affects the convergence rate of the implicit discretization scheme but not the accuracy of the solution. In the considered case, it has a negligible impact on the robustness of the algorithm because the radiation terms appear in the heat balance of the solid phase, which features the slowest dynamics. Nevertheless, the evaluation of the heat fluxes caused by radiation is extremely expensive. Computing the effective radiation term at each node requires numerical integration of Wien’s law along the entire length of the channel. Using the detailed radiation model leads to unacceptable CPU times for the dynamic simulation. Two alternatives have been considered in order to reduce the numerical effort. Truncation of the Range of Radiation. Instead of evaluation of the integral in eq 3 over the complete reactor length, the boundaries are restricted within the eff eff interval [z - seff rad, z + srad]. srad is determined heuristically in order to capture more than 90% of the actual


Figure 3. Effect of simplified expressions for the specific heat capacities (cp) of the process streams on the axial steady-state temperature profile (upper row) and the cp values (lower row) in the countercurrent reactor: reference solution with state-dependent cp values (left); simplified solution with constant cp values evaluated at an average temperature and gas composition (middle); simplified solution with constant cp values evaluated at inlet conditions (right).

radiative heat flux. This simplification leads to a significant acceleration of the evaluation of the radiation terms because L/seff rad ) O(100). Time-Delayed Update. The solution of the previous time step is used for computing the radiation term that is effective in the current time step. From the numerical point of view, the time-delayed update is equivalent to explicit time discretization of the radiation term. The radiation terms must be evaluated only once per time step, which significantly reduces the numerical effort. In this way only a small error is introduced in the transient behavior, which disappears in the steady-state solution. When these two simplifications are combined, the effect of radiative heat exchange can be considered without a noticeable additional numerical effort. 2.4. Implementation Details. The reactor model has been integrated in the in-house simulation tool IGS. IGS stands for interactive graphic system and has been developed especially for the purpose of interactive simulation. It provides the following features: (i) Online graphical output of simulation results with interactive selection of the output variables and the layout of the graphical output. (ii) Online user interaction for interrupting the numerical process in order to modify model and integration parameters and continuing or restarting the simulation run. (iii) IGS is linked during run time to the simulation model using interprocess communication in UNIX. The user has the choice to run the simulation stand-alone

Figure 4. Example for the use of interactive parameter manipulation in the IGS environment. The diagram on the left shows the variation of the critical model parameter, the inlet fuel concentration y°CH4 over the simulation time. The diagram on the right shows the temperature profile on the combustion side of the reformer immediately before varying y°CH4.

or with graphical output using the same executable. A comfortable user interface assists the integration of the simulation process in the IGS environment. (iv) IGS is compatible with arbitrary simulation codes and does not require formatted output data. Conceptual and implementation details of IGS are given in Appendix A. Online visualization and the optional interactive control of the simulation run have proven essential for developing and debugging the simulation program. Figure 4 shows a typical example, which illustrates the procedure for computing a meaningful reference solution


with the first, crude version of the code. When the simulation is started with the desired set of parameter values and with trivial initial conditions, e.g., uniform spatial profiles, the simulation most likely fails or becomes extremely time-consuming. This is due to the strong nonlinearity and the stiffness of the system. The only feasible way to find a solution is an embedding strategy regarding the critical model parameter (in the present case, the fuel feed concentration y°CH4). We start the simulation with a convenient value y°CH4 ) 6% instead of the targeted value y°CH4 ) 6.75%. The critical parameter is modified interactively toward the desired value while inspecting its effect on the solution. Simultaneous adaption of the spatial grid to the solution profile in PDEXPACK avoids failure of the numerical procedure. 2.5. Reactor Optimization Based on Dynamic Simulation. The criteria for rating the performance of the countercurrent reactor can be mapped by the following characteristic parameters:

fuel ratio: available heat of combustion maximum required heat of reforming

Figure 5. Temperature and conversion profiles of the reforming and combustion reactions of the optimal solution found by dynamic simulation of the reactor behavior.5

effectiveness factor: heat actually consumed for reforming η) heat actually released for combustion

two process streams significantly deviates from the optimal value.

Ψ)

Optimal conditions imply Ψ f 1+ and η f 1-. Conversion of the reforming and the combustion reaction higher than 99% are required as constraints, and a low maximum temperature is required as a secondary condition. The extension of the catalytic sections in both channels, the flow rate of the combustion gas, and the fuel concentration have been used as optimizing parameters. The used numerical model is too delicate in order to compute the required steady-state solutions for arbitrary parameter sets and therefore not suitable for an automated optimization procedure. A trial and error strategy is employed instead. Lacking a priori knowledge requires a broad search area in the parameter space in order to cover the optimum configuration. This is a grave problem regarding the efficiency of the procedure and the certainty of the result. Online visualization of the results has been a valuable source of information during this procedure in order to understand the crucial effects and to identify the key parameters for adjusting the process behavior. Figure 5 shows the result of this optimization procedure.5 The end of the catalytic zone on the reforming side is aligned with the position of the maximum temperature in order to prevent backreaction in the heat-exchanger zone. The catalytic zone and the flow rate on the combustion side are used to fix the position of the combustion zone. The inlet fuel concentration is adjusted in order to stabilize the ignited state with a minimum heat surplus. Approximately 76% of the heat generation is utilized by steam reforming. The excess heat is required as the driving force for heat recovery and exits the process with the product streams. Despite the significant advantage as compared with alternative concepts, this solution obviously does not fulfill the initially specified criteria (section 1): The two reaction zones tend to separate from each other, leading to an extreme maximum temperature. Additionally, the heat capacity ratio of the

3. Simplified Model At this point the question arises of whether a more favorable solution than the one discussed above can definitely be excluded. Therefore, further analysis aims at the computation of the complete set of solutions that fulfill the prespecified qualitative features. The detailed model is prohibitive for treating this issue: Its complexity hides the significance of crucial effects and impedes a clear interpretation of the obtained results. Furthermore, using the common dynamic simulation approach, it is not possible to implement design constraints directly but only through a superimposed trial and error procedure. A reduced model representing a best case scenario and neglecting secondary effects has been considered instead. The reduced model is based upon knowledge gained from the detailed model about the crucial features to be included. For the sake of simplification, a generic endothermic reaction A f B and an analogue exothermic reaction A f C are assumed. The equilibrium constant of the endothermic reaction is adopted from methane steam reforming. Likewise, catalytic and homogeneous paths are postulated for the exothermic reaction, and their kinetics are adapted to methane combustion. Further simplifying assumptions are as follows: (i) instantaneous endothermic reaction; (ii) catalyst applied on both sides of the separating wall; (iii) ideal heat transfer between heat sources (exothermic side) and heat sinks (endothermic side); (iv) negligible axial conductivity and diffusivity. The simplified model leads to a stationary boundary value problem in space given by the following equations. 3.1. Model Equations.

Energy balance on the endothermic side: m ˘ 1cp,1

dT1 ) -Rav(T1 - Tw) dz

(4)


Energy balance on the exothermic side: m ˘ 2cp,2

dT2 (5) ) -Rav(T2 - Tw) + 2(-∆HR2)rhom 2 dz

Energy balance of the wall: 0 ) Rav(T1 + T2 - 2Tw) + 1(-∆HR1)rcat 1 + (6) 2(-∆HR2)rcat 2 Component balance on the endothermic side: m ˘1

dw1,A ) -1MWArcat 1 dz

(7)

Component balance on the exothermic side: m ˘2

dw2,A hom ) -2MWA(rcat 2 + r2 ) dz

(8)

Chemical equilibrium condition: K(Tw) ) 1.0e(∆HR1/R)(1/T0-1/Tw) )

w1,B w1,A

(9)

Boundary conditions: T1(z)0) ) T01 T2(z)L) ) T02 w1,A(z)0) ) w01,A w2,A(z)L) ) w02,A

(10)

The model equations can be rearranged as follows in order to eliminate the unknown rcat 1 . Equations 6 and 7 lead to

Rav -∆HR2 cat m ˘ 1 dw1,A ) (T1 + T2 - 2Tw) + 2 r MW dz -∆HR1 -∆HR1 2 (11) Equation 9 leads to

Tw )

T0

(

RT1 w01,A + w01,B 1+ ln -1 -∆HR1 w1,A

)

(12)

The system consisting of eqs 4, 5, 8, and 11 forms an ODE system for the vector of unknowns u ) (T1, T2, w1,A, w2,A):

Conv

du )q dz

Figure 6. Illustration of the multiple-shooting algorithm. The length of the countercurrent reactor is divided into three subintervals. Integration direction: from the right to the left in intervals 1 and 3 and from the left to the right in interval 2. The defining equations of the problem are derived from boundary conditions at the nodes of the shooting algorithm or from constraint conditions. The unknowns represent variable model parameters and state variables at these fixed points.

(13)

However, it must be taken into account that the expression for Tw (eq 12) becomes nearly singular as (w01,A + w01,B)/w1,A f 1. This is the case at low wall temperatures and corresponding low conversion of the endothermic reaction. The z domain is divided into sections in order to overcome the problem. Equations 6 and 7 hold solely in the high-temperature sections. The endothermic reaction is completely neglected at low temperatures. Equations 11 and 12 are used there with rcat 1 explicitly set to zero. The significance of this procedure is explained in more detail in section 3.2.1.

3.2. Algorithm. The common procedure for solving stationary boundary value problems including discretization of the equations in space transforming the original ODE system into a large algebraic system is not adequate for the given system. The idealized model assumptions, e.g., negligible dispersive effects or instantaneous reactions, cause tremendous numerical problems: Because the model equations contain only first-order derivatives, the commonly used scheme of central differences tends to oscillations. Upwind discretization as an alternative would distort the solution by numerical diffusion. A shooting method16 has been applied instead, which is particularly suitable for the given task. 3.2.1. Multiple-Shooting Method. Figure 6 illustrates the method. The spatial domain is divided into subintervals defined by boundary points (z1, z4), internal boundaries, e.g., transition from catalytic to an inert zone (z3) and auxiliary fixed points (z2). The axial profiles are computed through numerical integration using a standard DAE solver. The direction of integration in each subinterval can be chosen in order to minimize the numerical effort. The initial conditions of each subinterval are formal unknowns of the system. They are specified through the boundary and continuity conditions at the boundaries of subintervals and superimposed constraint conditions. As an example, the conversion of the endothermic reaction is neglected in the first subinterval, which is considered as the preheating section in order to avoid a singularity in eq 12 (section 3.1). The extension of this interval is specified by the wall temperature at z2 (Tw ) 750 K). The switch of the model at z2 causes a negligible discontinuity in the axial profile of the mass fraction w1,A that jumps from w01,A (feed conditions) to a value corresponding to the equilibrium at Tw ) 750 K (∆w1,A/w01,A ) 0.5 × 10-3). The attractiveness of the shooting method results from the fact that the solution can be represented in various compatible forms of strongly differing dimensions: The full information about the steady-state solution is contained in the state variables at the nodes of the multiple-shooting algorithm. The underlying integration in z yields a very accurate solution for the axial profiles and a high spatial resolution. In this way it is possible to combine high accuracy with compactness, which is required by the superimposed processing


Figure 7. Temperature and conversion profiles obtained with perturbed initial conditions at the nodes of the multiple-shooting algorithm.

of the solution, i.e., parametric or stability analysis. Furthermore, design specifications can be implemented in an easy way. For each constraint introduced, a model parameter must be converted into a variable. In this way it is possible to restrict the range of possible solutions to those that are meaningful or practically relevant. However, applying the shooting method to countercurrent processes is numerically not a straightforward approach. The major difficulty is the sensitivity of the underlying initial value problem because one of its eigenvalues typically has a positive real part. Therefore, each deviation in the initial conditions increases exponentially along the integration path. Figure 7 illustrates this problem. Starting from a converged solution, the temperature values at the nodes of the multiple-shooting algorithm have been perturbed by ∆T ) 1 K. These small perturbations cause deviations from the solution profile up to 100 K at the end of the integration intervals. Larger perturbations may even lead to numerical overflow of the initial value problem. Adaptive adjustment of the extension of the subinterval and a good initial guess for the state variables of the boundary value problem are proper remedies for this problem. The former has been accomplished by a heuristic method of distributing the nodes for multiple shooting and the latter by embedding the boundary value problem in the continuation algorithm for parametric analysis. The continuation direction is determined by the tangent of the solution branch at the current solution.17 A heuristic step size control is applied based on the angle between the predictor and corrector steps. Details on the control criterion are given in Appendix B. 3.2.2. Stability Analysis. The algorithm has been extended in order to analyze the dynamic stability of the solution. This includes computation of the eigenvalues of the transient model equations in order to identify the dominating time constants of the system. The dynamic behavior in the present case is determined by the heat capacity of the wall. This is expressed by the following equation for the wall temperature:

w(Fc)s

∂Tw ) Rav(T1 + T2 - 2Tw) + 1(-∆HR1)rcat 1 + ∂t (14) 2(-∆HR2)rcat 2

Figure 8. Flowsheet of the procedure for parametric analysis based on a multiple-shooting algorithm.

The other variables are assumed to be quasi stationary. Hence, the stationary balance equations (4), (5), and (7)-(9) hold. The steady-state solution, which is already known from the previous procedure, can be used for the stability analysis. The spatial derivatives are approximated through central differences using a reasonable number of spatial nodes (≈30), leading to a DAE system in time:

B

du˜ ) -C‚Dz‚u˜ + q˜ dt

(15)

Equation 15 must be converted into normal form in order to determine the eigenvalues of the dynamic system. Details of the respective procedure are given in Appendix C. 3.2.3. Implementation Details. A multitasking approach is used for implementing the algorithm (Figure 8). The master process, written in Octave,18 includes the steady-state parametric and stability analysis. Octave enables a transparent and compact implementation of these algorithms through built-in functions for complex linear operations, e.g., factorization or computation of the eigenvalues or the null space of matrices, and provides according to our experience excellent numerical efficiency and robustness. Furthermore, Octave supports a modular implementation of complex algorithms by invoking external Fortran routines. The boundary value problem describing the reactor model is solved using the Gauss-Newton algorithm NLSCON.19 The residuals of the system equations are evaluated in an underlying process linking the solver of the boundary value problem with the reactor model. The interface routine translates the variables of the shooting algorithm in the context of the reactor model and triggers the integration of the stationary model equations over the specified subintervals of the spatial domain. LIMEX20,21 is used for integration of the resulting initial value problem. The relative integration accuracy is typically set to 10-9. 3.3. Reactor Analysis. To discuss the potential of the presented method, we pick up the question posed


Figure 9. Maximum wall temperature (upper diagram) and final reforming conversion (lower diagram) against the inlet fuel concentration. The circles indicate the effectively computed points. The miniaturized diagrams show the wall temperature profiles (solid line) and the reforming and combustion conversion profiles (dashed lines) at selected points of the bifurcation diagram. They illustrate the qualitative features of the steady-state solution on different branches of the bifurcation diagram.

at the beginning of section 3: Are the ideal conditions specified in section 1 attainable, in principle, with the countercurrent wall reactor? The bifurcation diagram has been constructed using the inlet fuel concentration as the distinguished continuation parameter. Figure 9 shows the maximum wall temperature and the final conversion of the endothermic reaction against the inlet fuel concentration in the combustion channel. Circles indicate the points that were actually computed. The step-size control described in Appendix B reduces automatically the step size close to the turning points and increases it in the linear parts. However, the continuation procedure stops if the shooting algorithm fails to converge. In most cases convergence is achieved after adapting the number or the position of the nodes. Reinitialization of the continuation procedure is manifested in small continuation steps in the sections a-c of the solution curve. The algorithm fails completely in section d. This is due to the extremely steep conversion profiles belonging to this part of the diagram. Stability analysis indicates that the steady-state solution in the interval AB is unstable. The leading eigenvalues are complex, with real parts alternating around the imaginary axis obviously because of numerical inaccuracies. The solutions found by dynamic simulation lie on the lowest (extinguished) and topmost (fully ignited) branches of the bifurcation diagram. Using the simplified model in conjunction with the continuation algorithm, a variety of qualitatively new solutions have been found. One additional stable branch exists in the interval CD, with solutions featuring overlapping of the endothermic and exothermic reaction zones. Of particular interest are the states close to the right turning point D, which combine high performance and thermal efficiency at moderate conditions: The effectiveness factor exceeds 90%, while the maximum temperature is identical with the equilibrium temperature determining the final conversion of the endothermic reaction. However, the respective states are very sensitive against perturbations, as indicated by their leading eigenvalue, which is very close to the imaginary axis. Therefore, their applicability is rather doubtful.

Figure 10. Integrated, autothermal methane reformer with cocurrent flow in the reaction zone and countercurrent heat recovery. Simulated temperature and conversion profiles.

4. Alternative Reactor Design: Cocurrent Mode The above result is not surprising if we recall that countercurrent cooling generally favors runaway of strongly exothermic reactions. An alternative reactor design is proposed, motivated by the fact that cocurrent cooling is much less susceptible to runaway. The new configuration is modular, consisting of two countercurrent heat-exchanger sections for heat recovery and a catalytic stage with cocurrent flow of the reforming and fuel gas4 (Figure 10). The additional fuel injector facilitates the control of the required high-temperature level. An inherent advantage of this design is that each process stream exchanges heat with itself, ensuring almost equal heat capacity fluxes for the countercurrent heat exchange. The computed temperature and conver-


sion profiles show the local overlapping of endothermic and exothermic reactions. In this way a moderate temperature profile results with cold ends and a hot section inside the reforming stage. More than 90% of the released heat of combustion is used for the production of synthesis gas. The reaction section is considered separately in order to assess the stability issue. The simplified model is described by the following set of equations:

Energy balance on the endothermic side: m ˘ 1cp,1

dT1 ) -R‚av‚(T1 - Tw) dz

(16)

Energy balance on the exothermic side: m ˘ 2cp,2

dT2 ) -R‚av‚(T2 - Tw) + 2(-∆HR2)rhom (17) 2 dz

Energy balance of the wall: 0 ) Rav(T1 + T2 - 2Tw) + w(-∆HR1)rcat 1 + (18) w(-∆HR2)rcat 2 Component balances on the endothermic side: m ˘1

dw1,j w ) avF1β(w1, j - w1, j) dz

j ) CH4, H2 (19)

w 0 ) -avF1β(w1,CH - w1,CH4) - wMWCH4 rcat (20) 1 4 w 0 ) (1 - X1eq)w1,CH w1,CH 4 4

(chemical equilibrium)

MWH2 w 0 ) w1,H w 3 Xeq w0 1,H 2 2 MWCH4 1,CH4 1,CH4

(21)

Component balances on the exothermic side: m ˘2

dw2,CH4 dz

w ) avF2β(w2,CH - w2,CH4) 4

(22)

Figure 11. Catalytic cocurrent wall reactor. Top: inlet temperature of reforming and combustion gas (upper diagram), fuel inlet concentration (middle diagram), and length of the reaction zone (lower diagram) against the mass flux of the reforming gas. Bottom: temperature and conversion profiles for three selected states of operation.

w - w2,CH4) - wMWCH4rcat (23) 0 ) -avF2β(w2,CH 2 4

X2 ) 1 -

Boundary conditions: Ti(z)0,t) ) T0i

i ) 1, 2

wi, j(z)0,t) ) wi,0 j

(24)

Initial conditions: Ti(z,t)0) ) T0i

i ) 1, 2

wi, j(z,t)0) ) wi,0 j

(25)

This model leads to a stationary DAE system with only first-order derivatives (initial value problem). Four design specifications are imposed additionally in order to restrict the set of solutions within the range of practical relevance. (i) Methane conversion of combustion and reforming are equal to 99%:

X1 ) 1 -

w1,CH4(L) 0 w1,CH 4

) 0.99

(26)

w2,CH4(L) 0 w2,CH 4

) 0.99

(27)

(ii) A temperature difference of ∆T ) 50 K is required between the inlet and outlet in both streams as a driving force inside the countercurrent heat exchanger:

T1(L) - T01 ) 50 K

(28)

T2(L) - T02 ) 50 K

(29)

For each constraint, a model parameter has to be converted into a variable. In the present case, the variable parameters are the inlet temperature of the two streams T01 and T02, the feed concentration on the 0 , and the reactor length L. The exothermic side w2,CH 4 resulting boundary value problem is treated by the procedure described in section 3.2. The upper diagrams of Figure 11 show the operating conditions and design options, which meet the introduced constraints, for variations in the mass flux of the reforming gas. The corresponding temperature and conversion profiles of three selected states of operation


are illustrated in the diagrams below. The characteristic cold spot is caused by the instantaneous reforming reaction and the assumed ideal heat transfer between the heat source and sink. The qualitative reactor analysis shows stable solutions with a favorable operation state over a wide range of operating conditions. However, depending on the kinetics, the heat-transfer characteristics, and the operating conditions, extinction or runaway of combustion may occur. 5. Conclusions The procedure for model-based analysis of a novel autothermal heat-exchanger reactor is elucidated. The close interdependence between an appropriate process model and its numerical solution is demonstrated. Crucial elements of transparent, reliable, and efficient numerical tools are state-of-the-art adaptive algorithms with integrated error control. Furthermore, online visualization of the solution profiles and interactive control of the simulation highly facilitates code development as well as process understanding. The IGS simulation environment is introduced, which supports the user in this task. Appropriate tuning of the nonlinear terms of the model is required in order to optimize its numerical solvability. The reaction rate expressions have been modified in order to eliminate singularities, which otherwise stop the numerical procedure. The numerical integration is significantly accelerated by neglecting the state dependence of heat capacities. However, this kind of simplification may distort the result and must be applied carefully. The detailed dynamic model provides a comprehensive picture of the predominant mechanisms of the considered process. For a rigorous parametric analysis, it is recommended to use a simplified and idealized version of the model that neglects secondary effects. This model requires an adequate numerical approach. The applied multipleshooting method combines a very accurate resolution of the solution profiles with a compact representation of the system for the superimposed parametric and linear stability analysis. Rigorous parametric analysis yields additional steady-state solutions besides those found by direct dynamic simulation. Some of them are very attractive but practically infeasible because of their parametric sensitivity. This allows one to conclude with a certain confidence that the countercurrent reactor concept is not suitable for coupling high-temperature steam reforming with in situ heat generation. On the other hand, the empirical knowledge gained enables us to modify the reactor design properly in order to overcome the identified shortcomings. Cocurrent flow of the reforming and fuel gas facilitates an efficient heat supply to the reforming reaction. Appendix A: The IGS Graphics IGS utilizes the multitasking capabilities of the UNIX operating system where different processes, e.g., simulation, graphics output, memory allocation, etc., run independently, synchronized by interprocess communication. Figure 12 displays schematically the communication architecture of IGS. The basic process in IGS is the so-called super manager, which triggers the simulation run and the graphics process. It further administrates the system resources for interprocess communication. The communication among different processes

Figure 12. Communication architecture of IGS: data exchange utilizing shared-memory program synchronization via messages and signals.

is accomplished via shared memory, messages, and signals. The graphics process is usually in the sleeping status in order to save CPU resources. In fact, the increase in processor load caused by IGS is usually negligible. It is activated or terminated by signals coming from the super manager or the simulation process. The simulation process receives messages from the user or from the super manager, which trigger several prespecified actions, such as starting, interrupting, or resuming simulation. On the other hand, messages from the simulation process give the super manager information on the simulation status. Shared memory segments are utilized for the exchange of large amounts of data among the simulation and the graphics process. Appendix B: Step-Size Control in the Continuation Algorithm Step size in a continuation scheme means the distance between two adjacent points of the solution branch. The selection of the step size is crucial for the efficiency and robustness of the algorithm: Small steps provide a safe convergence but are inefficient in regions of nearly linear behavior. On the other hand, large steps are much more efficient but critical regarding their convergence and may cause the failure of the procedure in regions of high parametric sensitivity. Therefore, stepsize control is useful in order to optimize the continuation algorithm with respect to efficiency and robustness. Figure 13 illustrates the strategy developed in the scope of the present paper: The control strategy is based on the angle formed by the direction ti of the initial guess ∆xjesti and the chord between the current and previous solution ∆xjsoli-1. ∆xjesti and ∆xjsoli-1 are normalized by dividing their elements by the respective scaling factors. If the system is close to linear, the directions of ∆xjesti and ∆xjsoli-1 almost coincide. The angle φ increases with growing parametric sensitivity of the system. Thus, φ can be taken as an indicator for the nonlinearity of the problem. The step size can be increased if φ is small

4806 Ind. Eng. Chem. Res., Vol. 43, No. 16, 2004 Chart 1

in order to convert the system to its normal form. Singular value decomposition of B yields

Figure 13. Illustration of the step size control of the continuation algorithm.

and must be reduced if φ is large. Accordingly, the step size is specified by the following formula:

cos φ )

〈ti, ∆xjsoli-1〉

The transformed system has the form given in Chart 1. The above system is explicitly separated into a dynamic and an algebraic part:

|ti| |∆xjsoli-1|

Sd‚x˘ ) Add‚x + Ada‚y

(33)

0 ) Aad‚x + Aaa‚y

(34)

∆xjiest ) 2 cos2 φ |∆xjsoli-1|

ti |ti|

The advantage of this strategy is its simplicity and its universal applicability. According to our experience, proper selection of the scaling factors for normalizing the elements of the vector of state variables is essential for the reliability of the step-size control. Appendix C: Transformation of a DAE into Normal Form Linearization of eq 15 yields

Equation 30 is a DAE system with a rank-deficient capacity matrix B. Further manipulation is necessary

where

( )

Sd 0 S) 0 0

(

Add Ada A) A A ad aa

)

The system given by eqs 34 and 35 can finally be transformed into its normal form:

The spectrum of eigenvalues of N determines the dynamic stability of the system. All linear operations


have been implemented in a convenient way using builtin functions of Octave.

Φ ) heat flow by radiation [W] Ψ ) fuel ratio η ) effectiveness factor

Nomenclature

Lower Indices

av ) specific surface area [m2/m3] cj ) concentration of component j [mol/m3] cp ) heat capacity of the gas phase [kJ/(mol K)] cs ) heat capacity of the solid phase [kJ/(kg K)] d ˜ ) effective channel diameter for radiative heat transfer [m] ej ) reaction order with respect to component j f i ) residual of the ith iteration in Newton’s scheme ki ) rate constant for equation i [1/s] m ˘ i ) mass flux [kg/(m2 s)] p ) distinguished continuation parameter pj ) partial pressure of component j [Pa] pmin ) threshold value of the partial pressure of compoj nent j [Pa] pnorm ) scaling value for the partial pressure of compoj nent j [Pa] q ) vector of source terms rtol ) relative integration tolerance ri ) reaction rate for equation i [mol/(m3 s)] seff rad ) effective radiation distance [m] t ) time [s] ti ) tangent to the solution curve u ) state vector u˜ l ) state vector of a spatially discretized, linearized system v ) state vector for the boundary value problem wj ) mass fraction of component j x ) normalized axial coordinate in Wien’s law xjesti ) estimation of the solution of the ith continuation step xi ) subvector of the dynamic state variables of a DAE system xji-1 sol ) solution of the (i - 1)th continuation step yi ) subvector of the algebraic state variables of a DAE system yj ) molar fraction of component j z ) axial coordinate [m] A ) surface area [m2] A, B, C ) model components A ) system matrix in a transformed system A ) linearized system matrix B ) capacity matrix Conv ) convection matrix D ) dispersion matrix Dz ) discretization scheme in spatial coordinate z ∆HR,i ) heat of reaction i [kJ/mol] Ki ) equilibrium constant of reaction i Kads ) adsorption constant for component j j L ) reactor length [m] MW ) molar mass [kg/kmol] N ) matrix of the normalized system R ) ideal gas constant [kJ/(mol K)] S ) diagonalized system matrix T ) temperature [K] U, V ) orthonormal transformation matrices X ) conversion Greek Letters R ) heat-transfer coefficient [W/(m2 K)] ) void fraction λ ) set of model parameters F ) density [kg/m3] σ ) Stefan-Boltzmann constant [W/(m2 K4)] φ ) deviation angle in continuation

a ) matrix entry index d ) matrix entry index est ) estimation i ) reaction, channel j ) component j n ) matrix entry index rad ) radiation sol ) solution w ) wall R, i ) reaction i Upper Indices ads ) adsorption cat ) catalytic eff ) effective eq ) chemical equilibrium hom ) homogeneous i ) ith iteration l ) linearized min ) minimum threshold norm ) scaling value s ) surface T ) transponed 0 ) inlet conditions Acronyms BVP ) boundary value problem CFD ) computational fluid dynamics CPU ) computational processing unit DAE ) differential algebraic equation IGS ) interactive graphic system

Literature Cited (1) Glo¨ckler, B.; Gritsch, A.; Morillo, A.; Kolios, G.; Eigenberger, G. Autothermal Reactor Concepts for Endothermic Fixed-Bed Reactions. Chem. Eng. Res. Des. 2004, 82 (A2), 148-159. (2) Zanfir, M.; Gravriilidis, A. Modeling of a catalytic plate reactor for dehydrogenation-combustion coupling. Chem. Eng. Sci. 2001, 56, 2671-2683. (3) Kolios, G.; Frauhammer, J.; Eigenberger, G. Autothermal Fixed-Bed Reactor Concepts. Chem. Eng. Sci. 2000, 55, 59455967. (4) Kolios, G.; Frauhammer, J.; Eigenberger, G. Efficient reactor concepts for coupling of endothermic and exothermic reactions. Chem. Eng. Sci. 2002, 57, 1505-1510. (5) Frauhammer, J.; Eigenberger, G.; von Hippel, L.; Arntz, D. A New Reactor Concept for Endothermic High-Temperature Reactions. Chem. Eng. Sci. 1999, 54, 3661-3670. (6) Kreul, L. U.; Fernholz, G.; Gorak, A.; Engell, S. Erfahrungen mit den dynamischen Simulatoren DIVA, gPROMS und ABACUSS. Chem. Ing. Tech. 1997, 69, 650-653. (7) Kro¨ner, A.; Holl, P.; Marquardt, W.; Gilles, E. D. DIVAs an Open Architecture for Dynamic Process Simulation. Comput. Chem. Eng. 1990, 14, 1289-1295. (8) Oh, M.; Pantelides, C. C. A Modelling and Simulation Language for Combined Lumped and Distributed Parameter Systems. Comput. Chem. Eng. 1996, 20, 611-633. (9) http://www.software.aeat.com/cfx. (10) http://www.fluent.com. (11) Frauhammer, J.; Klein, H.; Eigenberger, G.; Nowak, U. Solving moving boundary problems with an adaptive moving grid method-Rotary heat exchangers with condensation and evaporation. ZIB, Preprint SC 96-57 1996. (12) Nowak, U.; Frauhammer, J.; Nieken, U. A Fully Adaptive Algorithm for Parabolic Partial Differential Equations in One Space Dimension. Comput. Chem. Eng. 1996, 20, 547-561.

4808 Ind. Eng. Chem. Res., Vol. 43, No. 16, 2004 (13) Eigenberger, G.; Butt, J. B. A Modified Crank-Nicolson Technique with Non-Equidistant Space-Steps. Chem. Eng. Sci. 1976, 31, 681-691. (14) Xu, J.; Froment, G. F. Methane steam reforming, methanation and water-gas-shift: 1. Intrinsic kinetics. AIChE J. 1989, 35, 88-96. (15) Westbrook, C. K.; Dryer, F. L. Simplified reaction mechanisms for the oxidation of hydrocarbon fuels in flames. Comb. Sci. Technol. 1981, 27, 31-43. (16) Sto¨r, J.; Burlisch, R. Einfu¨hrung in die numerische Mathematik II. Numerical Analysis. In A First Course in Scientific Computation; Deuflhard, P., Hohmann, A., Translators; Springer: New York, 1989; p 17. (17) http//www.octave.org.

(18) Nowak, U.; Weimann, L. A Family of Newton Codes for Systems of Highly Nonlinear EquationssAlgorithm, Implementation, Application. ZIB, Technical Report TR 90-10 1990. (19) Deuflhard, P. Recent progress in extrapolation methods for ODEs. SIAM Rev. 1985, 27, 505-53. (20) Ehrig, R.; Nowak, U.; Oeverdieck, L.; Deuflhard, P. Advanced extrapolation methods for large scale differential algebraic problems. In High Performance Scientific and Engineering Computing; Springer: New York, 1999; pp 233-244.

Received for review October 1, 2003 Revised manuscript received May 27, 2004 Accepted May 31, 2004 IE034158E

Novel Reactor Concepts for Thermally Efficient Methane Steam

Recommend Documents