Microsoft Word - 476hernandez.docx CHEMICAL ENGINEERING TRANSACTIONS VOL. 43, 2015 A publication of The Italian Association of Chemical Engineering Online at www.aidic.it/cet Chief Editors: Sauro Pierucci, Jiří J. Klemeš Copyright © 2015, AIDIC Servizi S.r.l., ISBN 978-88-95608-34-1; ISSN 2283-9216 Bifurcation Analysis of Perfectly Stirred Reactors with Large Reaction Mechanisms Luigi Acampora*a, Erasmo Mancusia, Francesco S. Marrab aUniversità degli Studi del Sannio, Dipartimento di Ingegneria, Benevento, Italy bIstituto di Ricerche sulla Combustione – CNR, Napoli, Italy luigi.acampora@unisannio.it An accurate characterization of new generation fuels gets through a definition of the ignition and extinction conditions. According to the terminology of non-linear analysis the ignition and extinction conditions (critical conditions) are defined as the turning-points on the S-curve and are uniquely identified by two saddle-node bifurcations. The exact location of the critical points that mark the ignition and extinction conditions, the onset of instabilities, as well as the dependence of the location on parameter values, cannot be obtained easily using temporal simulations tools. For a systematic and accurate analysis, the bifurcation analysis and the parametric continuation technique are the tools of choice. Even when the reactive mixture is described by a simple surrogate, but in conjunction with very complex and detailed chemical mechanism, the computation of the critical conditions become computationally very demanding. Usually, these difficulties are overcome by using reduced schemes and allowing simplifications in the model, like constant averaged thermodynamic properties. In this work we explore these issues. Particularly, we show that the adoption of a Broyden type corrector in the continuation algorithm leads to performances that made affordable the computations even with desktop computers. Consequently, we introduce a suitable algorithm to investigate ignition, extinction and linear stability of the air-fuel mixtures in a Perfectly Stirred Reactor (PSR). The algorithm is based on the well-known Keller pseudo-arclength continuation method in order to compute steady state solution curves. The solutions stability and then ignition and extinction states are identified by test functions based on the numerical eigenvalues of the Jacobian of the governing system of equations. The algorithm is implemented in the numerical computing environment Matlab coupled with the CANTERA Toolbox for managing of complex chemical kinetic mechanisms. The algorithm is thus easily applicable to chemical schemes available in the standard ChemKin format. To illustrate the capability of the resulting method, the characteristic S-Shaped curve, including non-stable branches, for several air-hydrocarbon mixtures have been computed. 1. Introduction With the introduction of new generation fuels derived from alternative sources, the need to deal with very large chemical mechanisms is quickly increased. Indeed chemical composition of new fuels is usually very complex and their composition is likely to vary depending upon the properties of initial raw material. Several efforts are being conducted to properly characterize all the relevant properties of the new generation fuels (Biagini et al., 2013). Limit phenomena, such as ignition, are rather sensitive to chemical kinetics and these properties are therefore useful to physically characterize the behaviour of these fuels. In the framework of bifurcation theory, the ignition and extinction correspond to saddle-node bifurcation points in a classical S-shaped steady-state curve (Law, 2006). Therefore, the location of ignition and extinction conditions and their dependence on the main parameters (like pressure, equivalence ratio, residence time or inlet temperature in reactors) can be reformulated as a problem of bifurcation analysis. However, this DOI: 10.3303/CET1543147 Please cite this article as: Acampora L., Mancusi E., Marra F.S., 2015, Bifurcation analysis of perfectly stirred reactors with large reaction mechanisms, Chemical Engineering Transactions, 43, 877-882 DOI: 10.3303/CET1543147 877 approach represents a computationally difficult task when reaction mechanisms with several hundreds of species and thousands of chemical reactions are considered. In the past, a full bifurcation analysis has been successfully conducted adopting detailed schemes for some elementary fuels like hydrogen (Kalamatianos and Vlachos, 1995; Olsen et al., 1999), or methane with the simpler versions of the GRI mechanism (Park and Vlachos, 1997; Park and Vlachos, 1998) by using software written specifically for modelling combustion reactions. To avoid the repetitive encoding of subroutines for calculation of thermodynamic properties and reaction rates and managing of transport and kinetic database, a specific tool able to deal with these data has to be coupled with a bifurcation analysis module. The idea of coupling bifurcation analysis modules to chemical engineering software has been already successfully adopted by Vadapalli and Seader (2001). They modified a well-known predictor-corrector continuation algorithm for using Aspen Plus software as corrector. However, with this approach the Jacobian matrix is not available to the user, therefore the stability of steady state solutions and the exactly location of the bifurcation points cannot be determined. An alternative approach based on analytical Jacobian evaluation was recently proposed by Shan and Lu (2012). These authors presented the complete bifurcation map for methane (Shan and Lu, 2012), Di-Methyl Ether (DME) (Shan and Lu, 2012, Shan and Lu, 2014) and Jet Fuels (JF) (Shan and Lu, 2013) by using computer code that automatically generate mechanism-specific subroutines for analytical Jacobian evaluation from mechanisms described in CHEMKIN format. In this work Cantera (Goodwin, 2003), a suite of object-oriented software tools is adopted for the management of kinetic mechanism (including thermodynamic and transport properties) and the computation of the properties derived from the state variables (e.g. density, heat capacity). It is coupled with predictor-corrector pseudo-arclength continuation algorithm (Keller, 1987). Some basic elements of this algorithm are chosen specifically to reduce the huge computational cost in dealing with reaction mechanisms with several hundred of species and thousands of chemical reactions. The proposed approach has been applied to compare, ignition and extinction of three JF mixtures. They are complex mixtures of hundreds of hydrocarbons spanning a wide range of carbon numbers and chemical classifications (Dooley et al., 2010), posing a significant test case for the developed method. 2. Methodologies The equations of unsteady adiabatic constant pressure PSR can be written as: , , 1, 2, ,j j f j j j s dY Y Y W r j N dt Vτ ρ − …= + = (1) , , 1 ( )sN j f j f j j j j p pj Y h h h W rdT dt c Vcτ ρ =  − = −       (2) where , , , , , , , , , ,s pY t W r V N T h cτ ρ are the mass fraction, the time, the residence time, the molecular weight, the net production rate, the density, the volume of reactor , the number of chemical species, the mass specific enthalpy and constant pressure specific heat. The subscripts j and f refer to species index and to feed conditions. The stationary solutions curve (Seydel, 2010), including ignition and extinction states, has been identified by using a one-parameter continuation technique. Let ( ), pF x denotes the augmented system formed by the r.h.s. of the Eq(1) and Eq(2) and an additional scalar equation representing the chosen parameterization. The pseudo arc-length continuation algorithm (Keller, 1987) is used to numerically solve the steady state vector system: ( ), 0p =F x (3) where x is the state vector, 1[ , , , ]sNY Y T= …x and p the chosen parameter that is, in all the computations here reported, the residence time τ in the PSR. An appropriate choice of the following sub-modules is required for a successful continuation algorithm: • predictor; • corrector; • step-size control. 878 The first two of these three items can be chosen independently of each other. Instead the step-size control depends on the choice of the parameterization strategy, the predictor and the corrector (Seydel, 2010). There is not an ideal choice for all these items but the choice is problem dependent. For the problem at hand, the main obstacle besides in the high-dimensional chemical scheme formulation, that depends on the number of species and number of reactions. The corrector step is usually performed by means of a Newton-like method thus requiring evaluation of J , i.e. the Jacobian matrix of Eq(3). The number of entries of this matrix increases by increasing the number of species while the computation of each entry involves the very expensive evaluations of the reaction rates. Also expensive, being for this case the Jacobian matrix dense, is its factorization. To make affordable the computational cost, matrix factorization should be avoided together with a significant reduction of function F evaluations. This is achieved by adopting the so- called Broyden's method (Broyden, 1965; Kelley, 1995) that constructs a sequence of states that, starting from an initial guess 0Y and 11 0 0( ) −− ≈B J Y , converge to the solution of Eq(3) by using the iterative scheme: 1 1 ( ) ( )k k k k k − + = −Y Y B Y F Y (4) which correspond to a secant update, in conjunction with the formula (Broyden, 1965): 1 1 1 1 1 1 Tk k k kT k − − − − + − Δ − Δ = + Δ Δ Δ Y B F B B Y B Y B F (5) The cost of every Broyden step is approximatively equal to the cost of a vector function evaluation. It is locally super-linearly convergent (Kelley, 1995). The algorithm has been implemented by using the Matlab© programming environment in conjunction with the CANTERA open source library (Goodwin, 2003) to provide a very flexible and general tool for the description of the governing equations and an easy handling of large chemical kinetic mechanisms, including species thermodynamic and transport properties. In particular, CANTERA is used in order to manage the kinetic scheme and to compute , , ,pc h W r as function of the state of the system. 3. Results and discussion The procedure illustrated has been adopted to study the chemical stability of three different JF, respectively: • JET-A POSF 4658 • Shell GTL (Gas-To-Liquid) • Sasol CTL (Coal-To-Liquid) These fuels are complex mixture of several hundred hydrocarbons including alkanes, aromatics, and cycloalkanes or naphtenes (Dagaut and Cathonnet, 2006), that does not permit the complete development of a detailed reaction model (Kick et al., 2012). Therefore surrogates have been defined to study their combustion behaviour as reported in Table 1. Table 1: Composition of selected surrogates. Jet Fuel Reference Formula MW Surrogate [kg/kmol] Species mol% n-decane 42.67 JET-A Dooley et al., (2010) C8.60 H17.27 120.67 iso-octane 33.02 toluene 24.31 ShellTM GTL Kick et al., (2012) C9.81 H21.62 139.3 n-decane 90.60 iso-octane 9.40 n-decane 72.00 SasolTM CTL Kick et al., (2012) C9.59 H19.98 135.1 iso-octane 13.00 n-propylbenzene 15.00 Obviously, in order to describe the conversion of reactants to products, a detailed kinetic scheme is needed. To allow a consistent comparison of the different fuels, a single chemical mechanism should be adopted for all the three fuels. Such a scheme must be able at least to describe the combustion of any hydrocarbon that constitutes the surrogates, if not a specific validation for the adopted surrogate is available. A detailed scheme suitable for this purpose has been proposed by the Chemical Reaction Engineering and Chemical Kinetics 879 (CRECK) group, ver. 1212 (Frassoldati et al., 2010; Ranzi et al., 2012). It consists of 466 species and 14631 reactions. The one-parameter continuation of the stationary solution from cold conditions to ignited ones in a PSR with the residence time τ as parameter and inlet temperature of 900 K has been performed at different pressure conditions (10 atm and 50 atm), and with different air/fuel equivalence ratios AΦ (0.6 to 1.4 by 0.1). Figure 1a shows the equilibrium state temperature in the reactor as function of residence time for the three selected JF under pressure of 50 atm and inlet temperature 900KinT = . The curves are typical folded S- curves with multiple solutions and distinct ignition and extinction states (Law, 2006). Figure 1a confirms that the three fuels have similar behaviour. This result was expected, because the alternative fuels have been designed to replace the traditional jet fuels without any change of engine technology. The most relevant difference between the three jet fuels is at ignition and extinction states. In particular JET-A has ignition and extinction residence times greater than that of the alternative jet fuels. The stability of the three branches can be established by studying the sign of all eigenvalues a ibλ = + of the Jacobian matrix, according the Lyapunov's indirect method. This analysis is performed as a post processing of the S-curve of Figure 1a, allowing for a further reduction of the computational effort by properly choosing the points where the eigenvalues are required to not miss the detection of bifurcations. The eigenvalues can be sorted in descending order of the real part without loss of generality (Shan and Lu, 2012). They have been computed initially every 10 points. Then, in range where a change of sign of the first eigenvalue occurs, computation has been performed adopting information from all the computed solutions. Then, starting from the points enclosing the bifurcation, the continuation is repeated with a step size small enough to detect the bifurcation with the desired accuracy. The resulting picture for the JET-A fuel is shown in Figure 1b by reporting the two largest real part eigenvalues as function of residence time. An investigation of Figure 1b shows that the stability depends only on the first eigenvalue 1λ : it is the real part of this eigenvalue that changes sign at the bifurcation points. The behaviour for the other two fuels is very similar and it has not been reported in Figure 1b. The effect of mixture composition and pressure on extinction and ignition conditions can be estimated by collecting results from one-parameter continuations. Specifically, only the ignition and extinction as function of mixture composition and for two values of the pressure are reported in Figure 2a. We can observe that the ignition reactor temperature is almost insensitive to mixture composition while the extinction temperature is very sensitive to this parameter, with a maximum at 0.8AΦ = (fuel rich mixtures). Both decrease with decreasing pressure. While the difference in ignition temperature keeps almost the same over the whole range of mixture composition investigated (slightly decrease), the extinction temperatures clearly diverge by increasing AΦ . a) b) Figure 1: a) Reactor Temperature vs. Residence Time for a non-isothermal CSTR, at 50 atm, stoichiometric conditions, and 900 K inlet temperature. The circles represent the location of the fold bifurcations, the solid lines indicating stable branches, and dotted lines indicating unstable branches solutions. Symbols distinguishes among fuel types. Circle cross indicates bifurcation points. b) Real part of two largest real part eigenvalues as function of residence time (JET-A fuel). 880 The different behaviour at ignition of the three fuels can be highlighted by reporting the extinction and ignition residence time versus the air/fuel equivalence ratio. This is shown in Figure 2b. The ignition residence time slowly (almost linearly) grows with AΦ . Therefore a longer residence time is required to allow the activation of the reaction for leaner mixtures, and ignition happens when a critical temperature is reached, which is independent from mixture composition. The extinction residence time and reactor temperature behave almost mirror, i.e. when extinction residence time grows the extinction reactor temperature decrease and vice-versa. Further, the minimum of extinction residence time does not correspond to the maximum of extinction temperature. Indeed, the extinction residence time is minimum at air/fuel ratio 0.9 (slightly fuel rich mixtures). a) b) Figure 2: Extinction and ignition reactor temperature, a), and residence time, b), as function of air/fuel equivalence ratio at 10 atm and 50 atm with mixture inlet temperature 900 K. Consequently, we can deduce that the combustion of these fuels in the PSR is better sustained in conditions different from stoichiometry, in particular for slightly fuel rich conditions. Limiting residence times are more sensitive to pressure variations. Both largely decrease by increasing pressure. The differences are slightly larger at the lowest values of AΦ for the ignition residence time, while they are maximum at 0.9AΦ = for the extinction residence time. The substantial chemical equivalence of the traditional fuels and the new synthetic fuels is confirmed: appreciable differences arise especially at higher values of the pressure. 4. Conclusions A predictor-corrector continuation method suitable for the parametric study of fuel-air mixtures combustion modelled with large detailed kinetic mechanisms has been presented. In particular, we have shown the potential of the Broyden's method as corrector in the continuation of system of ODE arising from 0- dimensional models adopting complex chemical schemes. The numerical code has been developed in Matlab and coupled with the CANTERA libraries for the prompt inclusion of the chemical and thermodynamic data of the mixtures. The overall code results very compact and highly customizable, offering the opportunity of an easy extension to different reactor configurations or to the investigation of different mixtures compositions and chemical mechanisms. Actually, it has proven to be effective also in application to parametric continuation of problems governed by PDE equations (Mancusi et al., 2015). Finally, we have applied this code to the study of ignition, extinction of conventional and alternative JF-air mixtures in a non-isothermal adiabatic PSR. Several one-parameter continuations have been performed at different pressure conditions and with different air/fuel equivalence ratios by using the residence time as bifurcation parameter. Ignition and extinction points have been collected over a single map to show that, even if the general trends of conventional fuels are equally reproduced by the alternative fuels, some differences exist especially close to ignition and extinction conditions and that these differences are sensitive to pressure and, close to extinction, to air/fuel equivalence ratio. The analysis can be extended to verify how the choice of the surrogate or of the kinetic mechanisms affects the ignition and extinction conditions. 881 Acknowledgements The support of POR FSE Campania 2007-2013 is acknowledged. The authors are grateful to Prof. Antonio Viviani for his helpful comments. References Biagini E., Simone M., Barontini F., Tognotti L., 2013, A Comprehensive Approach to the Characterization of second Generation Biofuels, Chemical Engineering Transactions, 32, 853-858. DOI: 10.3303/CET1332143 Broyden C. G., 1965, A class of methods for solving nonlinear simultaneous equations, Mathematics of computation, 19(92), 577–593. Dagaut P., Cathonnet M., 2006, The ignition, oxidation, and combustion of kerosene: A review of experimental and kinetic modeling, 32, 48–92, DOI: 10.1016/j.pecs.2005.10.003 Dooley S., Won S. H., Chaos M., Heyne J., Ju Y., Dryer F. L., Kumar K., Sung C.-J., Wang H., Oehlschlaeger M. A., Santoro R. J., Litzinger T. A., 2010, A jet fuel surrogate formulated by real fuel properties, Combustion and Flame, 157(12), 2333–2339, DOI: 10.1016/j.combustflame.2010.07.001 Frassoldati A., Cuoci A., Faravelli T., Niemann U., Ranzi E., Seiser R., Seshadri K., 2010, An experimental and kinetic modeling study of n-propanol and iso-propanol combustion, Combustion and Flame, 157(1), 2– 16, DOI: 10.1016/j.combustflame.2009.09.002 Goodwin D. G., 2003, An open-source, extensible software suite for CVD process simulation, Chemical Vapor Deposition XVI and EUROCVD, 14(40), 2003-08. Kalamatianos S., Vlachos D. G., 1995, Bifurcation Behavior of Premixed Hydrogen/Air Mixtures in a Continuous Stirred Tank Reactor, Combustion Science and Technology, 109(1-6), 347–371. DOI:10.1080/00102209508951909 Keller H. B., 1987, Lectures on numerical methods in bifurcation problems. Springer Verlag, Heidelberg, West Germany Kelley C. T., 1995, Iterative Methods for Linear and Nonlinear Equations, Frontiers in Applied Mathematics. Society for Industrial and Applied Mathematics, Philadelphia, Pennsylvania, United States. Kick T., Herbst J., Kathrotia T., Marquetand J., Braun-Unkhoff M., Naumann C., Riedel U., 2012, An experimental and modeling study of burning velocities of possible future synthetic jet fuels, Energy, 43(1), 111–123, DOI: 10.1016/j.energy.2012.01.035 Law C. K., 2006, Combustion Physics. Cambridge University Press, Cambridge, United Kingdom Mancusi E., Acampora L., Marra F. S., Altimari P., 2015, Hysteresis in autothermal methane reforming over Rh catalysts: Bifurcation analysis, Chemical Engineering Journal, 262, 1052–1064. doi:10.1016/j.cej.2014.10.061 Olsen R. J., Vlachos D. G., 1999, A Complete Pressure−Temperature Diagram for Air Oxidation of Hydrogen in a Continuous-Flow Stirred Tank Reactor, The Journal of Physical Chemistry A, 103, 7990–7999, DOI: 10.1021/jp991148b Park Y. K., Vlachos D. G., 1997, Kinetically driven instabilities and selectivities in methane oxidation, AIChE Journal, 43(8), 2083–2095 DOI: 10.1002/aic.690430816 Park Y. K., Vlachos D. G., 1998, Isothermal Chain-Branching, Reaction Exothermicity, and Transport Interactions in the Stability of Methane/Air Mixtures, Combustion and Flame, 114, 214–230, DOI: 10.1016/S0010-2180(97)00285-X Ranzi E., Frassoldati A., Grana R., Cuoci A., Faravelli T., Kelley A., Law C., 2012, Hierarchical and comparative kinetic modeling of laminar flame speeds of hydrocarbon and oxygenated fuels, Progress in Energy and Combustion Science, 38(4), 468–501, DOI: 10.1016/j.pecs.2012.03.004 Seydel R., 2010, Practical Bifurcation and Stability Analysis, vol. 5: Interdisciplinary Applied Mathematics, Eds. Antman S., Greengard L., Holmes P., Springer, New York, NY, United States Shan R., Lu T., 2012, Ignition and extinction in perfectly stirred reactors with detailed chemistry, Combustion and Flame, 159(6), 2069–2076, DOI: 10.1016/j.combustflame.2012.01.023 Shan R., Lu T., 2013, Effects of Surrogate Jet-Fuel Composition on Ignition and Extinction in High Temperature Applications, in 8th U. S. National Combustion Meeting. Shan R., Lu T., 2014. A bifurcation analysis for limit flame phenomena of DME/air in perfectly stirred reactors, Combustion and Flame, 161(7), 1716–1723. DOI:10.1016/j.combustflame.2013.12.025 Vadapalli A., Seader, J. D., 2001, A generalized framework for computing bifurcation diagrams using process simulation programs, Computers & Chemical Engineering, 25(2-3), 445–464. DOI:10.1016/S0098- 1354(01)00624-X 882