Prediction of Multiple Steady States in Distillation through Simple

Nov 18, 2010 - Separation processes play a critical role in industry, accounting for 40−70% of ..... steady states through simple mass and heat bala...
2 downloads 0 Views 846KB Size
Download PDF
1346

Ind. Eng. Chem. Res. 2011, 50, 1346–1351

Prediction of Multiple Steady States in Distillation through Simple Mass and Heat Balance Analysis Manabu Kano,* Hiroshi Makita, and Shinji Hasebe Department of Chemical Engineering, Kyoto UniVersity, Kyoto, Japan

A novel method is proposed to judge the possibility of multiple steady states in distillation. Interestingly, the method is inspired by the analysis of multiple steady states in a nonisothermal continuous stirred tank reactor, in which the balance between heat generation and heat removal plays a critical role. In fact, the overall heat balance is visualized to predict multiple steady states in the developed method; any detailed static or dynamic simulation is not necessary. Several case studies demonstrate the usefulness of the proposed method. Introduction Separation processes play a critical role in industry, accounting for 40-70% of both capital and operating costs. In addition, distillation, which handles more than 90% of all separations,1 is the most common unit operation in the chemical process industry. Therefore, efficient design and operation of distillation processes are critical. Multiple steady states are phenomena such that several sets of dependent variables, e.g., product compositions, exist for a given set of independent variables, e.g., reflux and boilup flow rates. A well-known example is the multiplicity observed in exothermic chemical reactors, in which two stable and one unstable steady states exist.2 It is difficult but sometimes required in practice to operate a reactor at its unstable steady state. In such a situation, catastrophic jumps might be observed in operating conditions.3 Here, it is important that multiple steady states appear in distillation.4 Research concerning multiple steady states in distillation has a long history. In the last two decades, multiple steady states of simple (binary) distillation,5-7 azeotropic distillation,8-19 and reactive distillation20-27 have been actively investigated theoretically and also experimentally. In addition, there are several research works focusing on the operational issue28 and the startup of the distillation process with multiple steady states.16,17,29,30 More recently, multiple steady states in a heat integrated distillation column (HIDiC), which is a highly energy-efficient distillation process with a compressor, are analyzed.31 In the theoretical analysis of multiplicity, infinity/infinity analysis (infinite reflux and an infinite number of trays)8 is the most popular and has been widely used.9-11,13,22,23,14,15,18 Other approaches include the arclength continuation method,32 the interval-Newton method,33 bifurcation diagram/analysis,6,8,9,12,13,18,22,23,34 nonlinear wave propagation theory,6 and so on. However, these methods require a process model or simulator of a distillation column to judge the presence of multiple steady states. In other words, in conventional approaches, process simulations are repeated by changing manipulated variables until multiple steady states are found. In addition, all of the calculations have to be made again if design parameters such as the number of trays are changed. Clearly, such approaches are demanding. In the present work, a novel method is proposed to predict multiple steady states in distillation without detailed models or simulations. In fact, the proposed method can predict * To whom correspondence should be addressed. E-mail: manabu@ cheme.kyoto-u.ac.jp.

multiple steady states without information on the number of trays. The method is inspired by the analysis of multiple steady states in a nonisothermal continuous stirred tank reactor (CSTR) with cooling jacket, in which the balance between heat generation and heat removal plays a critical role.35 The usefulness of the proposed method is demonstrated through case studies. Modeling and Visualization of Heat Balance The thermal stability analysis of CSTR, which is based on material and energy balances, is first explained. Then, the concept is extended so that multiple steady states in distillation can be analyzed. Multiple Steady States in CSTR. The existence of multiple steady states in nonisothermal CSTR can be explained by using Figure 1. The straight energy balance line and the S-shaped material balance line are drawn in the T-xA plane.2 There are three intersections, that is, three solutions to the material and energy balance equations, A, B, and C. Each solution corresponds to a steady state. The steady states A and C are stable. However, B is an unstable steady state, because with a small rise in temperature the heat generation (material balance curve) is greater than the heat removal (energy balance curve). The excess heat generation increases the temperature until A is reached. An opposite situation happens when the temperature drops slightly below B. In this case, the temperature decreases until C is reached. In summary, any steady state appears when and only when the heat generation and the heat removal are balanced. Overall Heat Balance of Distillation. In the present work, a novel method is proposed to judge the possibility that multiple steady states appear in a targeted distillation process without detailed models. The proposed method does not need informa-

Figure 1. Multiple steady states in CSTR.

10.1021/ie1007823  2011 American Chemical Society Published on Web 11/18/2010

Ind. Eng. Chem. Res., Vol. 50, No. 3, 2011

QF )

{

L FH(x F,TF,PF) V FH(yF,TF,PF)

(1 -

1347

liquid

(5)

vapor

L rv)FH(T F,PF)

+

V rVFH(T F,PF)

mixed

L QD ) DH(x D,PD)

(6)

L QB ) BH(x B,PB)

(7)

Here, rv denotes vapor fraction. It is assumed that distillate and bottoms are liquid at the boiling point. To formulate reboiler heat duty, the energy balance equation around the reboiler is derived on the basis of Figure 3.

Figure 2. Binary distillation column.

tion about the number of trays and the feed tray; nevertheless, it can judge the potential possibility of multiplicity. In the proposed analysis, heat input Qin and output Qout around the entire distillation column are formulated as functions of light key component composition in distillate and bottoms xD and xB. They form curved surfaces and their intersection curve corresponds to steady states in the targeted distillation because they must be balanced at any steady state. The heat input Qin is the sum of heat of feed and reboiler duty; the heat output Qout is the sum of heat of distillate and bottoms and condenser duty. Although multiple steady states never appear in binary distillation with constant molar flow,4 simple distillation columns with ideal vapor-liquid equilibrium may display multiple steady-states from two different sources.5 The first type is found for columns with mass or volume inputs such as mass reflux. The second type depends on the presence of an energy balance in the model. A molar basis is popular in academia, but in industry most streams are not given on a molar basis but rather on a mass or volume basis. Therefore, in this work, two types of control structures are investigated. (1) LV configuration: reflux flow rate L and boilup flow rate V on a molar basis are manipulated variables. (2) LWQreb configuration: reflux flow rate LW on a mass basis and reboiler heat duty Qreb are manipulated variables. The binary distillation column investigated here is shown in Figure 2. LV Configuration Case. In the proposed method, the overall heat input Qin and output Qout are calculated as functions of L, V, xD, xB, column pressure P, and feed condition F, xF, TF, and PF. No heat loss from the column is assumed. For binary distillation, the material balance equations around the entire column are given as follows:

L V L 0 ) LN-1H(x + Qreb - VH(y - BH(x N-1,PN-1) N,PB) B,PB)

(8)

Unknown variables are yN, LN-1, and xN-1. yN is calculated from physical property estimation formula as a function of xB and PB. LN-1 and xN-1 can be calculated from known variables through the material balance equations around the reboiler. 0 ) LN-1 - V - B

(9)

0 ) LN-1xN-1 - VyN - BxB

(10)

LN-1 ) V + B

(11)

Consequently,

xN-1 )

VyN + BxB V+B

(12)

Similarly, to formulate condenser heat duty, the energy balance equation around the condenser is derived on the basis of Figure 3. V L 0 ) V2H(y + Qcnd - (L + D)H(x 2,P2) D,PD)

(13)

Unknown variables y2 and V2 can be calculated from known variables through the material balance equations around the condenser. 0 ) V2 - L - D

(14)

0 ) V2y2 - (L + D)xD

(15)

V2 ) L + D

(16) (17)

Consequently,

0 ) FxF - DxD - BxB

(1)

y2 ) xD

0)F-D-B

(2)

Finally, Qin and Qout are formulated as follows when feed is liquid.

From these two equations, D and B can be expressed as:

L H(x ) N-1,PN-1)

x F - xB D) F xD - xB

(3)

xD - xF F xD - xB

(4)

B)

L V L L Qin ) FH(x + V(H(y - H(x ) + B(H(x F,TF,PF) N,PB) N-1,PN-1) B,PB)

The heat of feed, distillate, and bottoms are expressed by using liquid and vapor molar enthalpy, HL and HV.

(18)

L L V L Qout ) DH(x + BH(x + (L + D)(H(x - H(x ) D,PD) B,PB) D,P2) D,PD) (19)

where D and B are expressed by using eqs 3 and 4. LWQreb Configuration Case. In industry, most streams are not given on a molar basis but on a mass or volume basis. Therefore, the formulation of Qin and Qout is extended for LWQreb configuration.

1348

Ind. Eng. Chem. Res., Vol. 50, No. 3, 2011

Figure 3. Binary distillation column.

In LWQreb configuration, the reboiler heat duty Qreb can be directly specified. As a result, L Qin ) FH(x + Qreb F,TF,PF)

(20)

On the other hand, with molar weight MW(xD) of distillate, Qout is formulated as follows: L L + BH(x + Qout ) DH(x D,PD) B,PB)

(

)

LW V + D (H(x D,P2) MW(xD) L ) (21) H(x D,PD)

Visualization. Given product compositions xD and xB, the overall heat input Qin and output Qout can be calculated. To visualize the heat balance of the distillation column, curved surfaces representing Qin and Qout are drawn in xD-xB-Q space. The visualization procedure is as follows. (1) Set feed conditions and assume column pressure P. (2) Set independent variables: (L, V) or (LW, Qreb). (3) Calculate Qin and Qout for various sets of xD and xB. (4) Draw curved surfaces representing Qin and Qout in xDxB-Q space. An example of this visualization is shown in Figure 4 (top), where methanol-propanol separation is targeted. Any steady state is on the intersection curve, which is the boundary in Figure 4 (bottom), because Qin and Qout must be balanced at any steady state. The boundary corresponds to the trajectory of product compositions when the height of packed column is continuously changed. On the other hand, the product compositions realized by the tray column appear discretely along the boundary. Method for Analyzing Multiplicity In this section, the necessary conditions for the existence of multiple steady states are derived on the basis of the visualization technique proposed in the previous section. In addition, the usefulness of the proposed method is demonstrated through case studies. Necessary Conditions. Different values of manipulated variables result in different product compositions. A steady state moves on the xB-xD plane when one of the manipulated variables, L, V, LW, or Qreb, is changed. The trajectory of steady states has the following characteristics. First, xD increases as xB increases: the trajectory increases monotonically on the xB-xD plane regardless of the number of trays. Second, the slope decreases monotonically as xB increases. Third, the trajectory moves toward the point (0,1) when the number of trays increases and separation becomes easier accordingly. These characteristics are easily confirmed through rigorous simulations.

Figure 4. Visualization of heat balance. (Top) 3D plot. (Bottom) 2D plot.

Figure 5. Necessary conditions that multiple steady states exist.

Figure 5 illustrates three steady-state trajectories (dashed lines) drawn together with the boundary of Qin ) Qout. Each trajectory is drawn by changing one of the manipulated variables; different trajectories correspond to different numbers of trays. Steady states must exist at the intersection points (black points) of the trajectory and the boundary of Qin ) Qout. Therefore, the necessary condition for the existence of multiple steady states is that there is a region where the boundary has positive slope, because there is only one intersection point when the boundary decreases monotonically as in Figure 4. This necessary condition can be checked without any process simulation; the calculation and visualization of Qin and Qout are sufficient. Figure 5 indicates that multiple steady states appear only when the number of trays is not too large or too small. In the past, the condition for the existence of multiple steady states in binary distillation was derived from the different viewpoint.5 For distillation columns without multiple steadystates, the mole fraction of the most volatile component in distillate xD is expected to increase when the reflux flow rate L increases and the boilup vapor flow rate from a reboiler V is fixed. Therefore, a negative slope, ∂xD ∂L

|

0

(23)

V

where D denotes a distillate flow rate. This condition means that an unstable operating point appears when the increase of L causes the increase of D. Later, the similar condition for the existence of multiple steady states was derived for a heat integrated distillation column (HIDiC).31 These researches have clarified that the multiplicity appears when the relationship between key variables has opposite sign. Of course, the condition depends on design and operation as well as physical properties of mixtures. Case Study 1: Single Steady State. To test the validity of the proposed method, it is applied to binary distillation of methanol-propanol with LV configuration. The operating conditions are summarized in Table 1. Here ∆P is the estimated pressure drop at each tray. Qin and Qout are calculated by setting V ) 100 kmol/h and L ) 65, 70, and 75 kmol/h. The visualization results shown in Figure 6 (top) indicate that multiple steady states will not appear in this distillation because the boundaries at various reflux flow rates have negative slopes.

Figure 7. Rigorous steady-state simulation results: Case study 2.

To confirm this prediction, rigorous steady-state simulations are executed by using the RadFrac model on the platform of Aspen Dynamics (Aspen Technology, Inc.). The number of trays is fixed at six and the feed tray is the fourth. A large number of simulations with various reflux flow rate L are conducted. The results have clarified that multiple steady states do not exist. The trajectory is drawn in Figure 6 (top). The black points are the steady states at L ) 65, 70, and 75 kmol/h. Case Study 2: Multiple Steady States. In the next case study, Qin and Qout are calculated by setting V ) 300 kmol/h and L ) 298.5 kmol/h. The other conditions are the same as those in case study 1. The visualization result shown in Figure 6 (bottom) indicates that multiple steady states will appear in this distillation because the boundary has positive slope when xB is between 0.05 and 0.25. In fact, this prediction was verified by rigorous steady-state simulations; multiple steady states appear when the number of trays is six while they do not appear when the number is five as shown in Figure 7. These case study results clarify that the existence of multiple steady states depends on both design and operation and also that the proposed method is useful to judge the possibility of multiple steady states through simple mass and heat balance analysis. Changes in Steady States

Figure 6. Visualization and rigorous simulation results: Case studies 1 and 2.

In the previous section, the possibility of multiple steady states was judged by visualizing curved surfaces of both Qin and Qout, in particular, the boundary of Qin ) Qout. The dynamics of a targeted distillation column as well as the possibility of multiple steady states would be predicted without drawing boundaries many times if the influence of manipulated variables on the boundary was known in advance. In this section, the influence is investigated and the validity of prediction is confirmed through rigorous dynamic simulations. Influence of Manipulated Variables on Boundary. The influence of manipulated variables on curved surfaces of Qin and Qout can be calculated by partially differentiating Qin and

1350

Ind. Eng. Chem. Res., Vol. 50, No. 3, 2011

Figure 9. Dynamic simulation results: Case study 3. Figure 8. Visualization and rigorous simulation results: Case study 3.

Qout with respect to the manipulated variables. The LV configuration is investigated here; eqs 18 and 19 are partially differentiated with respect to L or V. For simplicity, P is assumed constant. First, the influence of L is studied. ∂Qin )0 ∂L

(24)

∂Qout V L ) H(x - H(x >0 D,P) D,P) ∂L

(25)

These partial differential coefficients mean that Qin is unchanged and Qout increases as L increases. As shown in Figure 6 (top), the area of Qin < Qout becomes larger when L increases. Next, the influence of V is studied. L

∂H(xN-1,P) ∂Qin V L ) H(y H (V + B) ,P) (x ,P) N N-1 ∂V ∂V

V L ) H(y - H(x - (V + B) N,P) N-1,P)

V L - H(x ) H(y N,P) N-1,P)

(26)

L ∂H(x ∂xN-1 N-1,P)

∂xN-1

∂V

L xB) ∂H(xN-1,P)

B(yN V+B

∂xN-1

∂Qout )0 ∂V

(27) In most cases, the third term of 26 is much smaller than the first and the second terms, and vapor molar enthalpy is greater than liquid molar enthalpy. Thus, the sign will be positive. As a result, Qin increases and Qout is unchanged when V increases. In other words, the area of Qin < Qout becomes smaller when V increases. In LWQreb configuration case, similar results can be obtained more easily.

Case Study 3: Influence of L. The binary distillation of methanol-propanol with LV configuration is investigated again. The settings are the same as those in case study 2. A large number of steady-state simulations are conducted while L is changed and V is constant at 300 kmol/h. The results shown in Figure 8 (top) have clarified that multiple steady states exist. In addition, Qin and Qout are calculated at L ) 298.5 and 299.0 kmol/h. The visualized boundaries of Qin ) Qout in Figure 8 (bottom) also indicate the existence of multiple steady states. Three steady states, i, ii, and iii, exist at L ) 298.5 kmol/h; similarly, three steady states, I, II, and III, exist at L ) 299.0 kmol/h. First, consider the case where L increases stepwise from 298.5 to 299.0. Since large L makes xD higher, it is predicted that steady states move such that i f I, ii f I, and iii f III. In addition, as expected, the area of Qin < Qout becomes larger when L increases as shown in Figure 8 (bottom). Next, consider the case where L decreases stepwise from 299.0 to 298.5. Contrary to the previous case, it is predicted that steady states move such that I f i, II f iii, and III f iii. It is obvious that the steady states ii and II are unstable. To confirm these predictions, rigorous dynamic simulations are executed by using Aspen Dynamics. The step responses are shown in Figure 9. In both cases, the simulation results are corresponding to the predictions. Similar results can be obtained from the proposed analysis, steady-state simulations, and dynamic simulations in the case where V is changed. It is important that any rigorous steady-state or dynamic simulation is not required to judge the possibility of multiple steady states. The proposed approach based on simple mass and heat balances can provide useful information about the multiplicity. Conclusions A novel prediction method for multiple steady states in distillation was proposed, and its validity was confirmed through

Ind. Eng. Chem. Res., Vol. 50, No. 3, 2011

case studies. The proposed method focuses on overall heat balance around the column and does not require rigorous distillation models or design parameters (the number of trays). Nevertheless, the proposed method can judge the possibility of multiple steady states accurately. It should be pointed out here that the prediction results depend on the physical properties of mixtures. If the physical properties are different, then the predicted and real results would become different. The proposed approach, which is inspired by the thermal stability analysis of a nonisothermal continuous stirred tank reactor (CSTR), is quite unique and different from conventional techniques. It is interesting to notice the analogy between exothermic reaction and distillation. However, the target of the present work is limited to binary distillation; the further extension is necessary to cope with more complex distillation processes. Literature Cited (1) Humphrey, J. L.; Keller, G. E. Separation Process Technology; McGraw-Hill: New York, 1997. (2) Levenspiel, O. Chemical Reaction Engineering, 3rd ed.; Wiley: New York, 1998. (3) Poston, T.; Stewart, I. Catastrophe Theory and Its Applications; Pitman: London, 1978. (4) Doherty, M. F.; Perkins, J. D. On the Dynamics of Distillation ProcessessIV Uniqueness and Stability of the Steady State in Homogeneous Continuous Distillations. Chem. Eng. Sci. 1982, 37, 381–392. (5) Jacobsen, E. W.; Skogestad, S. Multiple Steady States in Ideal TwoProduct Distillation. AIChE J. 1991, 37, 499–511. (6) Kienle, A.; Groebel, M.; Gilles, E. D. Multiple Steady States in Binary Distillation-Theoretical And Experimental Results. Chem. Eng. Sci. 1995, 50, 2691–2703. (7) Koggersbol, A.; Andersen, T. R.; Bagterp, J.; Jorgensen, S. An Output Multiplicity in Binary Distillation: Experimental Verification. Comput. Chem. Eng. 1996, 20, S835–S840. (8) Bekiaris, N.; Meski, G. A.; Radu, C. M.; Morari, M. Multiple Steady States in Homogeneous Azeotropic Distillation. Ind. Eng. Chem. Res. 1993, 32, 2023–2038. (9) Bekiaris, N.; Meski, G. A.; Morari, M. Multiple Steady States in Heterogeneous Azeotropic Distillation. Ind. Eng. Chem. Res. 1996, 35, 207– 227. (10) Bekiaris, N.; Morari, M. Multiple Steady States in Distillation: Infinity/Infinity Predictions, Extensions, and Implications for Design, Synthesis, and Simulation. Ind. Eng. Chem. Res. 1996, 35, 4264–4280. (11) Guttinger, T. E.; Dorn, C.; Morari, M. Experimental Study of Multiple Steady States in Homo-Geneous Azeotropic Distillation. Ind. Eng. Chem. Res. 1997, 36, 794–802. (12) Muller, D.; Marquardt, W. Experimental Verification of Multiple Steady States in Heterogeneous Azeotropic Distillation. Ind. Eng. Chem. Res. 1997, 36, 5410–5418. (13) Esbjerg, K.; Andersen, T. R.; Muller, D.; Marquardt, W.; Jorgensen, S. B. Multiple Steady States in Heterogeneous Azeotropic Distillation Sequences. Ind. Eng. Chem. Res. 1998, 37, 4434–4452. (14) Bekiaris, N.; Guttinger, T. E.; Morari, M. Multiple Steady States in Distillation: Effect of VL(L)E Inaccuracies. AIChE J. 2000, 46, 955– 979. (15) Gaubert, M. A.; Gerbaud, V.; Joulia, X.; Peyrigain, P. S.; Pons, M. Analysis and Multiple Steady States of an Industrial Heterogeneous Azeotropic Distillation. Ind. Eng. Chem. Res. 2001, 40, 2914–2924. (16) Scenna, N. J.; Benz, S. J.; Franceseoni, J. A.; Rodriguez, N. H. Start-up of Homogeneous Azeotropic Distillation Columns with Multiple Steady States. Ind. Eng. Chem. Res. 2004, 43, 553–565.

1351

(17) Scenna, N. J.; Benz, S. J.; Rodriguez, N. H.; Klaric, J. I. Startup of Homogeneous Azeotropic Distillation Sequences with Multiple Steady States. Ind. Eng. Chem. Res. 2005, 44, 9208–9220. (18) Kannan, A.; Joshi, M. R.; Reddy, G. R.; Shah, D. M. MultipleSteady-States Identification in Homogeneous Azeotropic Distillation Using a Process Simulator. Ind. Eng. Chem. Res. 2005, 44, 4386–4399. (19) Lee, H. Y.; Huang, H. P.; Chien, I. L. Design and Control of an Acetic Acid Dehydration Column with p-Xylene or m-Xylene Feed Impurity. 2. Bifurcation Analysis and Control. Ind. Eng. Chem. Res. 2008, 47, 3046–3059. (20) Nijhuis, S. A.; Kerkhof, F. P. J. M.; Mak, A. N. S. Multiple SteadyStates during Reactive Distillation of Methyl tert-Butyl Ether. Ind. Eng. Chem. Res. 1993, 32, 2767–2774. (21) Eldarsi, H. S.; Douglas, P. L. Methyl-tert-butyl-ether catalytic distillation columnsPart I: Multiple Steady States. Chem. Eng. Res. Des. 1998, 76, 509–516. (22) Guttinger, T. E.; Morari, M. Predicting Multiple Steady States in Equilibrium Reactive Distillation. 1. Analysis of Nonhybrid Systems. Ind. Eng. Chem. Res. 1999, 38, 1633–1648. (23) Guttinger, T. E.; Morari, M. Predicting Multiple Steady States in Equilibrium Reactive Distillation. 2. Analysis of Hybrid Systems. Ind. Eng. Chem. Res. 1999, 38, 1649–1665. (24) Higler, A. P.; Taylor, R.; Krishna, R. Nonequilibrium Modelling of Reactive Distillation: Multiple Steady States in MTBE Synthesis. Chem. Eng. Sci. 1999, 54, 1389–1395. (25) Chen, F. R.; Huss, R. S.; Doherty, M. F.; Malone, M. F. Multiple Steady States in Reactive Distillation: Kinetic Effects. Comput. Chem. Eng. 2002, 26, 81–93. (26) Huang, K.; Nakaiwa, M.; Tsutsumi, A. Towards Further Internal Heat Integration in Design of Reactive Distillation ColumnssPart II. The Process Dynamics and Operation. Chem. Eng. Sci. 2006, 61, 5377– 5392. (27) Lee, H. Y.; Tang, Y. T.; Huang, H. P.; Chien, I. L. Bifurcation in the Reactive Distillation for Ethyl Acetate at Lower Murphree Plate Efficiency. J. Chem. Eng. Jpn. 2006, 39, 642–651. (28) Jacobsen, E. W.; Skogestad, S. Multiple Steady-States and Instability in DistillationsImplications for Operation and Control. Ind. Eng. Chem. Res. 1995, 34, 4395–4405. (29) Benz, S. J.; Scenna, N. J. An Extensive Analysis on the Start-up of a Simple Distillation Column with Multiple Steady States. Can. J. Chem. Eng. 2002, 80, 865–881. (30) Benz, S. J.; Francesconi, J. A.; Scenna, N. J. Strategies for Starting up Azeotropic Distillation Columns with Multiple Steady States. Can. J. Chem. Eng. 2004, 82, 920–929. (31) Kano, M.; Fukushima, T.; Makita, H.; Hasebe, S. Multiple SteadyStates in a Heat Integrated Distillation Column (HIDiC). J. Chem. Eng. Jpn. 2007, 40, 824–831. (32) Morgan, A.; Sommese, A. Computing All Solutions to Polynomial Systems Using Homotopy Continuation. Appl. Math. Comput. 1987, 24, 115–138. (33) Kearfott, R. B., III. M. N. Algorithm 681 INTBIS, a Portable Interval Newton/Bisection Pack- age. ACM Trans. Math. Software 1990, 16, 152–157. (34) Doedel, E.; Champneys, A. R.; Fairgrieve, T. F.; Kuznetsov, Y. A.; Sandstede, B.; Wang, X. AUTO97: Continuation And Bifurcation Software for Ordinary Differential Equations; Con-cordia University: Montreal, Canada, 2007. (35) van Heeden, C. Autothermic Processes. Ind. Eng. Chem. 1953, 45, 1242–1247.

ReceiVed for reView April 1, 2010 ReVised manuscript receiVed September 21, 2010 Accepted October 26, 2010 IE1007823