Autoignition and burning rates of fuel droplets under microgravity

43 (2 ing icro ris, Dipartimento di Chimica, Materiali e Ingegneria Chimica, Politecnico di Milano, 20133 Milano, Italy Received 10 May 2004; received in revised form 2 June 2005; accepted 3 June 2005 Available online 20 July 2005 Abstract This paper presents a mathematical model for the unsteady evaporation, ignition, and combustion of isolated fuel droplets under microgravity. The model consists of a large structured system of differential algebraic equa- tions, where the numerical complexity is due both to the stiff nature of the kinetic mechanism and to the flame structure around the droplet. A very general, detailed kinetic scheme, consisting of ∼200 species and over 5000 reactions, is used to describe the gas-phase combustion of different fuels. Several comparisons with experimental measurements, carried out under various operating conditions, confirm that the proposed model is a useful tool for characterizing low-temperature and high-temperature ignition delay times. The predicted explosion diagrams of n-alkanes, as well as their ignition delay times, agree with the experimental measurements; the various oxida- tion regions are closely reproduced too. In addition to this, recent experimental results, relating to the influence of the initial diameter on droplet burning rates in cold and hot environments, are also presented and discussed. Lastly, an analysis of the extinction diameters for the combustion of n-heptane droplets allows a discussion of the role of radiative heat transfer, as well as further emphasizing the importance of the low-temperature oxidation mechanisms.  2005 The Combustion Institute. Published by Elsevier Inc. All rights reserved. Keywords: Fuel droplets; Microgravity; Burning rates; Detailed kinetics; Autoignition 1. Introduction Because of their high-energy density per unit vol- ume, the combustion of liquid fuels is of great interest in many practical applications from industrial burners to diesel engines. The characterization of this process is complicated by the strong interactions between sev- eral physicochemical processes dominating the igni- tion of fuel droplets. In fact, it is necessary not only to * Corresponding author. Fax: +39 02 7063 8173. E-mail address: [email protected] (E. Ranzi). account for the heating of the liquid droplet, its evap- oration, and finally the ignition delay time of the fuel vapor, but also for the interactions between the differ- ent droplets of a spray. It is thus useful to take apart the overall system and to study simpler and possibly ideal conditions. The chemistry of the ignition and autoignition of the pure components and their mix- tures is only a first step toward a proper understanding and the optimal design of various combustion devices and chambers. Higher combustion efficiencies offer a way forward for sustainable development and eco- nomic use of fuels, while reducing the production of greenhouse gases, particularly CO . The formation Combustion and Flame 1 Autoignition and burn under m A. Cuoci, M. Mehl, G. Buzzi-Ferra 0010-2180/$ – see front matter  2005 The Combustion Institute. doi:10.1016/j.combustflame.2005.06.003 005) 211–226 www.elsevier.com/locate/combustflame rates of fuel droplets gravity T. Faravelli, D. Manca, E. Ranzi ∗ 2 Published by Elsevier Inc. All rights reserved. and It is quite evident that the burning rates are con- trolled mainly by coupled diffusion and heat trans- fer, especially radiative heat transfer [2,3]. Therefore combustion chemistry plays only a marginal role. Nevertheless, the use of a detailed kinetic scheme, which also includes low-temperature oxidation mech- anisms, becomes necessary not only when describ- ing the typical ignition delay times and the explo- sion diagrams, but also when searching for a better understanding of burning rates under specific condi- tions [4]. At low temperatures, typically when the droplets are in an atmosphere at a controlled temperature, it is possible to observe and verify the characteristic fea- tures of the homogeneous oxidation of a hydrocarbon. Thus, the presence of cool flames, a negative temper- ature coefficient (NTC), or a zero temperature coef- ficient (ZTC) region, as well as two-stage ignitions around a fuel droplet, are a clear indication of the im- portance of the detailed chemistry of peroxide com- The general mathematical model used to describe the transient evaporation, ignition, and combustion of isolated pure fuel droplets in a large environment at constant pressure is briefly reported below. Micro- gravity simplifies the problem, because of the absence of buoyancy and the resulting 1-D nature of the evap- oration and combustion of a droplet. Time-dependent, spherically symmetrical numerical computations are well suited to such a problem. The creation and loca- tion of a droplet can be a source of asymmetry in ex- perimental measurements. Nevertheless, the predic- tions indicate that these uncertainties do not greatly affect the model’s assumptions. Detailed chemical ki- netics in the gas phase with multicomponent molec- ular transport are considered. Only the ordinary dif- fusion term of the mass flux of species is taken into account in the gaseous phase. Estimations using the model, including the Soret effect, did not show any 212 A. Cuoci et al. / Combustion Nomenclature Cp, C heat capacity d droplet diameter D mass diffusivity f fugacity H enthalpy k thermal conductivity Nc number of components P pressure Pv vapor pressure r radial coordinate Rd droplet radius R∞ reference volume radius T temperature t time vg radial gas velocity of pollutants (NOx, SOx, unburned and polyaromatic hydrocarbons and particles) is also influenced by the relative role of these different combustion processes. For this reason, the study of a single fuel droplet is significant in providing a basis for refining the exist- ing phenomenological understanding of burning and extinction processes, while allowing a quantitative as- sessment of the theoretical ability to predict such phe- nomena [1]. Moreover, microgravity conditions fur- ther simplify the overall process, because of the spher- ical symmetry of the system. The 1-D description of this system makes the adoption of detailed chemistry more viable. ponents [5,6]. Finally, detailed kinetic information is also required to characterize pollutant formation dur- ing the whole combustion process. Flame 143 (2005) 211–226 Vg gas diffusion velocity Vl droplet volume W mass flow rate Y mass fraction �Hv heat of vaporization ρ density ω reaction rate qR radiative flux Superscripts and subscripts 0 initial conditions a ambient conditions g gas phase l liquid phase i ith species Several efforts have been already devoted to ana- lyzing the combustion of fuel droplets [7,8], but only a few of them also fully describe the low-temperature mechanisms. Consequently, certain specific aspects of this subject still require further investigation [5,9,10]. All of these reasons justify the development of this mathematical model, dealing with the complexity of the numerical problem, arising mainly from the large dimensions of the adopted kinetic scheme. 2. Equations for the model and numerical methods significant modification of the results and for this rea- son are neither reported nor discussed. Given these assumptions, the liquid- and gas-phase conservation and ilibri 1 − T T 1C4t2 t = ( ) Vapor pressure PL [Pa] [14] exp C1 + C2T + C3 ln(T ) + C4T 2 n-Heptane C1 = 87.829, C2 = −6996.4, C3 = −9.8802, C4 = 7.2099 × 106 n-Decane C1 = 112.73, C2 = −9749.6, C3 = −13.245, C4 = 7.1266 × 106 Heat of vaporization �Hv [J/kmol] [14] C1 ( 1 − T Tc )C2 n-Heptane C1 = 5.0014 × 10−7, C2 = 0.38795 n-Decane C1 = 6.6126 × 10−7, C2 = 0.39797 equations of mass, species, and energy describe the system [2,11]. 2.1. Liquid phase Mass conservation: (1)ρ¯l dRd dt + Rd 3 dρ¯l dt = −ρg ( vg − dRddt ) . Energy conservation: (2)ρ¯lCl ∂Tl ∂t = 1 r2 ∂ ∂r ( r2kl ∂Tl ∂r ) . dRd/dt describes the movement of the droplet’s sur- face. This term can be positive (i.e., the droplet en- larges) if the expansion, due to heating, is greater than the mass loss due to evaporation. The average liquid density, defined as a volumetric average ρ¯l =∫ Vl ρl dV/Vl, is updated at each integration step. The properties of the liquid (ρl, kl,Cl,Pv,�Hv) are as- Specific attention should be devoted to the pres- ence of internal liquid motion due to heterogene- ity in the surface temperature [8]. Such internal re- circulation inside a pure fuel droplet enhances heat transfer. The heating up transient becomes slower, while the burning rate or extinction phenomena are not affected. To estimate the effect of this convective two-dimensional contribution is quite complex and would require further computational efforts, making the whole problem quite difficult for current com- puters to handle. On the basis of the proposals of Marchese et al. [2], it is possible to use an “effec- tive” liquid-phase conductivity, but this would in- troduce a new parameter. Generally good agreement with experimental measurements is observed, even when the presence of internal recirculation in the liquid phase is ignored. All the reported results as- sume unmodified transport properties. Only in the case of the analysis of sensitivity to transport proper- A. Cuoci et al. / Combustion Table 1 Thermophysical properties of liquid-phase and gas–liquid equ Properties Ref. Correlation Liquid density ρL [kmol/m3] [14] C1 C2 [ 1 + (1 − T C3 )C4 ] Liquid thermal conductivity kL [W/mK] [13] 3.5 × 10−3 T 1.2 ebN√ PMT 0.167c ( Liquid heat capacity CpL [J/mol K] [14] C21 t + C2 − 2C1C3t − C − C3C42 t4 − C0.24 t5, C1 + C2T + C3T 2 sumed from Daubert and Danner’s data base [12], from Reid et al. [13], and Perry and Green [14]; they are summarized in Table 1. Flame 143 (2005) 211–226 213 um n-Heptane C1 = 0.61259, C2 = 0.26211, C3 = 540.2, C4 = 0.28141 n-Decane C1 = 0.42831, C2 = 0.25745, C3 = 617.7, C4 = 0.28912 c )0.38( T Tc )−1/6 n-Heptane TebN = 371.58 K, Tc = 540.26 K n-Decane TebN = 447.30 K, Tc = 618.45 K n-Heptane − C2/33 t3 + · · · C1 = 61.260, C2 = 3.1441 × 105, 1 − T Tc C3 = 1824.6, C4 = −2.5479 × 103 n-Decane C1 = 2.7862 × 105, C2 = −197.91, C3 = 1.0937 ties, discussed in the coming paragraphs, is effective thermal conductivity introduced and explicitly indi- cated. and 214 A. Cuoci et al. / Combustion 2.2. Gas phase Mass conservation: (3)∂ρg ∂t + 1 r2 ∂ ∂r (r2ρgvg) = 0. Species conservation: ρg ( ∂Yg,i ∂t + vg ∂Yg,i ∂r ) = − 1 r2 ∂ ∂r (r2ρgYg,iVg,i ) (4)+ ωg,i , i = 1, . . . ,Nc. Energy conservation: ρgCp,g ( ∂Tg ∂t + vg ∂Tg ∂r ) = 1 r2 ∂ ∂r ( r2 ( kg ∂Tg ∂r )) − ρg Nc∑ i=1 (Yg,iVg,iCp,g,i ) ∂Tg ∂r (5)− Nc∑ i=1 (ωg,iHg,i ) − ∇qR. The transport properties were either taken from the CHEMKIN transport database [15] or estimated following the procedure proposed by Wang and Frenklach [16]. The thermochemical information on the gas phase was obtained primarily from the CHEMKIN thermodynamic database. Unavailable thermodynamic data were estimated by group addi- tivity methods [17]. Closure of the system and equilibrium between gas and liquid demand an equation of state, gener- ally expressed as f (V,T ,P,Yl,Yg) = 0, for the gas phase. The results reported below refer to the ideal gas equation. However, modified Redlich–Kwong– Soave and Peng–Robinson equations were also tested, without any significant modifications in the simula- tions. Radiative heat fluxes were estimated by adopt- ing the approach proposed by Kazakov et al. [18]. The Planck mean absorption coefficients were derived from polynomial expressions [19] for CO, CO2, and H2O, which are assumed to be the only molecules in- fluencing radiative transfer [18]. 2.3. Boundary conditions r = 0: ∂Tl ∂r = 0, (6)∂Yl,i ∂r = 0, i = 1, . . . ,Nc. r = Rd(t): kl ∂Tl ∂r ∣∣∣∣ R−d + Nc∑ i=1 (WtotYg,i + ρgYg,iVg,i )�Hv,i Flame 143 (2005) 211–226 (7)= kg ∂Tg ∂r ∣∣∣∣ R+d + qR, (8)Tl|R−d = Tg|R+d , (9)fl,i = fg,i , i = 1, . . . ,Nc, WtotYl,i − ρlDl,i ∂Yl,i ∂r ∣∣∣∣ R−d = WtotYg,i + ρgYg,iVg,i , (10)i = 1, . . . ,Nc, (11)Wtot = ρg ( vg − dRddt ) . r = R∞: (12)Tg = T 0g , (13)Yg,i = Y 0g,i , i = 1, . . . ,Nc. 2.4. Initial conditions (t = 0) (14)Rd = R0d . (15)r �R0d: Tl = T 0l , (16)Yl,i = Y 0l,i , i = 1, . . . ,Nc. (17)r > R0d: Tg = T 0g , (18)vg = 0, (19)Yg,i = Y 0g,i , i = 1, . . . ,Nc. The initial temperature in the gas phase differs ac- cording to the kind of experiment to be modeled. In the case of autoignition, a uniform temperature pro- file was assumed. Ignition experiments, on the other hand, require a more complex description. In fact, the deployment of the droplet into a cold environment is usually followed by the local application of a transient external ignition source (sparks or hot wires). These procedures can modify the spherical symmetric struc- ture of the problem. In line with the approach pro- posed by Marchese et al. [2], ignition was simulated by assuming a short period (∼1 s) of pure evaporation at ambient temperature, followed by a non-uniform temperature radial profile peaking at ∼2000 K. Fig. 1 shows this initial temperature profile, where the high- temperature conditions are between 1 and 5 radii from the droplet’s surface. 3. Numerical methods The overall model consists of a differential system of the conservation equations (1)–(5) with boundary and initial conditions (6)–(19). As previously noted, the attention to low-temperature phenomena leads to the use of a detailed kinetic mechanism, which in- cludes a large number of components. The resulting system of algebraic and partial differential equations and A. Cuoci et al. / Combustion Fig. 1. Initial conditions of the gas-phase temperature for the simulation of spark-ignited experiments. requires particular attention to the numerical aspects of the problem. The partial differential equations are discretized by means of a moving spatial grid. The nu- merical problem is then characterized by a large sys- tem of differential-algebraic equations (DAE), with the velocity in the gas phase evaluated using the mass conservation equation. Furthermore, nonlinear alge- braic equations are used to describe the gas–liquid interphase conditions. The complexity of this problem, coupled to the in- trinsic stiffness [20] of the DAE system, means that specific attention must be paid to the numerical meth- ods and solver routines. In fact, this problem is quite challenging, not only with regard to the precision re- quired, but also in terms of: • robustness (sudden ignition, steep profiles, high gradients); • efficiency (simulations can take more than 2 days of CPU time). We adopted the BzzMath library [21–23] compris- ing the BzzDae C++ classes (BzzMath is freeware for noncommercial use and can be directly down- loaded at http://www.chem.polimi.it/homes/gbuzzi). Appendix A addresses in some detail the main numer- ical features of this solver. 3.1. The spatial grid The droplet first enlarges when heated and then shrinks during vaporization. Therefore, a moving spatial grid was introduced to describe the moving boundaries between the droplet and the surrounding gas. Thirty equispaced grid points describe the liquid phase, while 100–300 discrete shells, with variable step sizes, characterize the gaseous phase. By apply- Flame 143 (2005) 211–226 215 ing the simple rule �ri+1 = α�ri , with α ∼= 1.05, the final step size becomes ∼500 times the initial one, de- pending, of course, on the total number of grid points. Hence, a large gas volume around the droplet is de- scribed by maintaining a very fine grid size at the gas–liquid interface, where the critical gradients are located. 3.2. The structure of the DAE system The structure of the DAE system is mainly a tridiagonal block, due to spatial discretization, with the exception of the column corresponding to the droplet’s radius, where the nonzero elements refer to the moving grid. Radiation terms in the energy balance too are spread over the whole matrix, even though diagonal terms dominate. Most of the equa- tions are devoted to the chemical species involved in the detailed kinetic scheme, with ∼200 equations for each discretization point. Consequently, the over- all number of DAEs ranges between ∼25,000 and 65,000 according to the adopted grid. The possibil- ity of exploiting the quasi-tridiagonal block structure is of crucial importance in drastically reducing CPU time. For this purpose, a specifically developed C++ DAE class of the BzzMath library was adopted (see Appendix A [21,23]). 4. The kinetic model A very general, detailed kinetic scheme, consist- ing of ∼200 species and more than 5000 reactions, was used to describe combustion in the gas phase. The details of the scheme are discussed elsewhere [24]. Only a few steps are reported here to highlight the important role of peroxy radicals. Because of its large dimensions, the whole mechanism is not re- ported. However, it is available in CHEMKIN for- mat, with the corresponding thermodynamic data and transport properties, upon request. The kinetic model simultaneously covers all the different typical operat- ing ranges of hydrocarbon oxidation, particularly the high- and low-temperature regions. At high tempera- tures (>1000 K), the decomposition of alkyl radicals (R) and the successive molecular decomposition of alkanes largely prevail on the different reaction paths. The transition between low- and high-temperature mechanisms is determined by the equilibrium: R + O2 ↔ ROO. (A) Below 850–900 K, the direct addition of O2 to the alkyl radical is favored, and the resulting peroxy and 216 A. Cuoci et al. / Combustion radical undergoes the branching path through the fol- lowing reaction channels: ROO ↔ QOOH (internal isomerization), QOOH + O2 ↔ OOQOOH (O2 addition), OOQOOH ↔ OH + OQOOH (isomerization–decomposition), OQOOH ↔ OH + products (chain branching). The propagation reactions consuming fuel oc- cur via H-atom abstraction primarily by OH and HO2, and to a lesser extent by H, CH3, and O2. Specifically, OH is the prevalent abstracting agent be- low 700 K, while at higher temperatures HO2 also plays a significant role. Peroxy radicals (ROO) iso- merize to give hydroperoxy-alkyl radicals (QOOH), which can add O2 and form peroxy alkyl hydroper- oxy radicals (OOQOOH). The successive decompo- sition of OOQOOH generates the ketohydroperoxy- des (OQOOH) and OH radical. Lastly, OQOOH is a degenerate chain-branching agent, since it is a very unstable intermediate, which decomposes to two rad- icals in a unimolecular decomposition. The shift in equilibrium (A) is the primary reason for the existence of a NTC region (negative tempera- ture coefficient, i.e., the range in which the conversion decreases as the temperature is increased) in the over- all reaction rate at temperatures above 650 K. There is a reduction of reactivity not only because the predom- inant branching reaction via the sequence is displaced in favor of the essentially nonbranching route, but more importantly because OH is considerably more reactive than HO2. Enhanced formation of HO2 can also be attributed to the decomposition of QOOH to- ward the conjugate alkene. In fact, at intermediate temperatures (between 700 and 900 K), the alkyl hy- droperoxy radical may decompose in three different reaction paths: QOOH ↔ HO2 + conjugate alkene, QOOH ↔ OH + O-heterocycle, QOOH ↔ OH + β-scission products (aldehyde + smaller alkanes). Above ∼900 K, alkyl radicals are favored and this leads to two important reaction routes. Either the di- rect abstraction of an H atom occurs via R + O2 ↔ HO2 + conjugate alkene, which gives HO2 as the primary propagating species, or there are β-scission reactions of the alkyl radical (R) with the formation of alkenes and smaller alkyl radicals. Flame 143 (2005) 211–226 5. Results and discussion 5.1. Modeling predictions Two examples of predictions are reported in Figs. 2 and 3. They refer to two typical experimen- tal conditions, relating to spark ignition at ambient temperature and to an autoignition at an intermedi- ate temperature (633 K). Both simulations refer to experimental measurements on n-decane droplets of ∼1.5 mm diameter performed by Xu et al. [4]. Ra- dial profiles of the temperature and concentration of relevant species in the gas are reported in Figs. 2a and 2b. The temperature profile in the spark-ignited experiments shows that [O2] becomes zero where the maximum temperature is reached. A typical diffusive burning regime is then observed. The concentration of normal-decane at the interface is dictated by vapor– liquid equilibrium and rapidly decreases to zero in the pyrolysis region, well before the maximum tempera- ture. The main products are H2O, CO, and CO2 under these conditions. Ethylene and acetylene are impor- tant intermediates, together with methane, benzene, and naphthalene, as shown in Fig. 2b. Moreover, the maximum temperature, initially located very close to the droplet, moves toward larger standoff ratios. CO2 and H2O are the end products in the complete com- bustion of the fuel droplet, but naphthalene is also predicted as a minor pollutant in the gas phase. The kinetic scheme used for these simulations does not account for soot formation, but naphthalene is the heaviest component and a clear indicator of the pos- sibility of soot formation. A very different situation is observed in Fig. 3, where the autoignition of n-decane with ambient air at 633 K results in a cool flame with a maximum tem- perature of ∼780 K. Normal-decane does not fully react under these conditions and ketohydroperoxide (OQOOH) is the main component with H2O, H2O2, and CO. Normal-decane and oxygen concentration profiles are smoother than in the previous case, as a consequence of the lower reactivity of the system. This lower reactivity means that complete combustion of the high-temperature oxidation mechanism cannot be activated, with the result that partially premixed combustion occurs. This system displays oscillating cool flames, as is clearly shown in Fig. 4, where the time evolution of the maximum temperature in the gas phase is re- ported. After an initial induction time of ∼0.5 s, an initial cool flame with a maximum temperature of 825 K is observed. Successive dumped cool flames then appear between 0.6 and 0.8 s, after which the system reaches a practically steady temperature of ∼780 K. This oscillating behavior can be easily ex- plained on the basis of the low-temperature oxidation and A. Cuoci et al. / Combustion Fig. 2. Calculated gas-phase species and temperature pro- file surrounding n-decane droplets at different times (a, b); time evolution of the droplet diameter squared and flame diameter (c) (spark ignition simulation, d0 = 1.46 mm, TG = 298 K, time = 1.20 s). mechanism [25], where a limited increase in tempera- ture reduces the importance of the branching reactions of ketohydroperoxides with a decrease in the overall reactivity of the system (NTC), as will be more clearly illustrated in the next paragraph. Flame 143 (2005) 211–226 217 Fig. 3. Calculated gas-phase species and temperature profile surrounding n-decane droplets at different times (a, b); time evolution of the droplet diameter squared and flame diame- ter (c) (autoignition simulation, d0 = 1.57 mm, TG = 633 K, time = 2.30 s). In order to verify the convergence of the numeri- cal method, different solutions are reported by vary- ing the grid points in the gas phase from 100 to 200 and 250. Small deviations are observed when the number of grid points is doubled. However, the behavior of the system remains unchanged and the and 218 A. Cuoci et al. / Combustion Fig. 4. Time evolution of the predicted maximum gas-phase temperature for different grid number of points in the gas phase (a); time and radial evolution of the predicted max- imum gas-phase temperature (b) (autoignition simulation, n-decane droplet, d0 = 0.91 mm, P = 0.1 MPa, TG = 633 K). solutions converge for 200 and 250 points. Further- more, the vaporization rates practically coincide in the three simulations. CPU times for these simula- tions are more than 2 days on a Pentium IV processor at 2.8 GHz, when 250 points are considered. In order to further confirm the absence of numerical instabil- ities in these results, Fig. 4b shows the time evolu- tion of the radial gas temperature profiles, where the cool flames are more evident. The main features of these low-temperature oxidation mechanisms will be further analyzed in the next paragraph, where com- parisons with experimental measurements will be also discussed. 5.2. Explosion regions and ignition delay times The oxidation and combustion of the fuel droplets display all the typical behavior of the homogeneous combustion of a hydrocarbon. In fact, in a pressure– Flame 143 (2005) 211–226 Fig. 5. Regions of the explosion diagram of n-decane droplets in air (d0 = 0.70 mm). Comparison between ex- perimental [26] (lines) and predicted (points) data. Fig. 6. Time evolution of the maximum gas-phase temper- ature for different initial ambient temperatures (n-decane droplets, d0 = 0.700 mm, P = 0.3 MPa). Acronyms are the same as in Fig. 5. temperature diagram, it is possible to observe a slow reaction region at low temperatures, and also the ex- plosion region with the typical one-step or hot igni- tion at high temperatures and pressures. Moreover, the combustion of fuel droplets also allows the singling out of the regions related to cool flames and to two- stage ignition, which are the result of the rapid transi- tion between low- and high-temperature mechanisms. According to the experimental measurements of Moriue et al. [26], Fig. 5 shows reasonable agreement between the experimental and the predicted autoigni- tion regions for n-decane droplets in air under micro- gravity. Examples of different profiles of the maxi- mum temperature in the gas phase around the droplet are reported in Fig. 6 to demonstrate the behavior and A. Cuoci et al. / Combustion of the system under the different operating condi- tions. All such complex behavior can be explained on the basis of the same critical reaction steps, i.e., the competition between the low- and high-temperature mechanisms. As already noted, cool flames are caused by ketohydroperoxide decomposing at relatively low initial temperatures. In fact, these branching reactions increase the reactivity of the system with a conse- quent sudden, but limited, temperature increase. The higher temperature favors peroxy radical decompo- sition with a drop in overall reactivity, due to the reduced importance of the low-temperature mecha- nism. If the system does not effectively exchange heat with the surrounding environment, the temper- ature remains constant via a second induction time, which leads to high-temperature ignition. This two- stage ignition (SI) is then observed when the tem- perature and/or pressure of the reacting system are high enough to make oxidation reactions competitive with heat transfer. If, on the other hand, the chem- istry is slower than the cooling, the heat flow lowers the system’s temperature to such an extent that suc- cessive dumped cool flames (CF) once again become possible. Finally, at ambient temperatures higher than ∼900 K, the high-temperature mechanism drives the system directly to a hot ignition. At higher pressures, it is also possible to observe a further zone of high-temperature ignition (HI) at low temperatures, as clearly shown in Fig. 5. This behav- ior can be very simply explained by the kinetic model on the basis of a two-stage ignition, where the sec- ond induction time becomes so limited that it is no longer possible for it to be observed experimentally. This situation is due to the relatively higher intermedi- ate temperatures reached by the system, after the first ignition. Fig. 7 shows a comparison between modeling pre- dictions and experimental measurements of ignition times for n-decane droplets in air with an initial di- ameter of 0.7 mm at 298 K and 0.3 MPa. The ignition time of cool flames, as well as the total ignition time, is properly predicted, even though there are slight overestimates of the first and total ignition delay times at low temperatures. It is interesting to observe that the effect of different temperatures and stoichiome- tries around the droplet is to smooth the NTC re- gion. The average of these different conditions easily leads to a more or less pronounced ZTC (zero tem- perature coefficient), where the ignition delay times remain unchanged when the temperature increases. This ZTC behavior is more evident for n-heptane droplets. A similar agreement between modeling and experimental measurements [5,6] for n-heptane and n-dodecane is discussed elsewhere [11]. It is worthwhile observing that the low-tempera- ture mechanism is also still effective at relatively high Flame 143 (2005) 211–226 219 Fig. 7. First and total ignition time for different initial ambient temperatures. Comparison between experimental [26] (points) and predicted (lines) data (n-decane droplets, d0 = 0.70 mm, P = 0.3 MPa). temperatures. In fact, the initial temperature of the droplet quenches the ambient gas around the droplet where peroxy radicals play a significant role. Thus the direct use of the high-temperature mechanism would produce higher ignition delay times. 5.3. Influence of droplet diameter on burning rates Experimental measurements of the autoignition and spark ignition of n-decane droplets with vari- ous initial diameters at atmospheric pressure have been performed [4]. The fuel droplets, which had an initial temperature of 298 K, were obtained from a syringe and suspended on a quartz fiber. Model- ing predictions are now compared with experimen- tal measurements in order to investigate the overall burning rate under different conditions. A first se- ries of data refers to the autoignition of n-decane droplets in air at 633 K. As already discussed and presented in Fig. 3, the relatively low ambient temper- ature prevents high-temperature ignition occurring. Under these conditions, the burning rates are con- ditioned by the low-temperature mechanism. Fig. 8 shows that the model properly agrees with the exper- imental measurements of the fuel’s burning rates for the three different diameters. The typical two-phase evolution of the droplet’s size is reported for conve- nience in terms of the diameter squared. For an initial period of about ∼1 s, the heating and expansion of the droplet prevail on the evaporation phase. Only after this first step can the linear dependence of d2 on time be observed. The slight increase in the vaporization rate when the droplet’s diameter rises is also notewor- thy. These vaporization rates are indeed the result of the low-temperature oxidation mechanism. This fact and 220 A. Cuoci et al. / Combustion Fig. 8. Droplet diameter squared for different initial di- ameters during autoignition experiments. Comparison be- tween experimental [4] (points) and predicted (lines) data. Plain line, low- and high-temperature mechanisms; dashed line, high-temperature mechanism only (n-decane droplets, P = 0.1 MPa, TG = 633 K). Fig. 9. Time evolution of the predicted maximum gas phase temperature for different initial droplet diameter (n-decane droplets, P = 0.1 MPa, TG = 633 K). is clearly singled out in Fig. 8, where the broken lines indicate the predicted vaporization rates when only the high-temperature reaction mechanism is adopted; i.e., the n-heptyl radical addition on an oxygen mole- cule to form peroxy radicals is ignored and is con- sidered to decompose immediately into small radicals and alkenes. Low-temperature chemistry dictates the vaporization rates and the time evolution of the maxi- mum gas-phase temperatures, for the three different diameters, as reported in Fig. 9. The rapid oscilla- tions in the maximum temperature observed in Fig. 9 are not the result of numerical instabilities as already Flame 143 (2005) 211–226 Fig. 10. Droplet diameter squared for different initial di- ameters during autoignition experiments. Comparison be- tween experimental [4] (points) and predicted (lines) data (n-decane droplets, P = 0.1 MPa, TG = 1093 K). observed in Fig. 4b, where the 3D map of the tem- perature profile is shown for a droplet with an initial diameter of 1.22 mm. All these profiles show the initial dumped cool flames and also the attainment of a more stable tem- perature plateau. The dumping of cool flames, clearly shown also in Fig. 4, is the combined result of the progressive consumption of fuel and gas-phase heat- ing. Once the first flame ignites, fuel is consumed and is only partially replaced by relatively slow vapor- ization. The successive cool flames exhibit a higher frequency, produce less heat, and gradually move the system from the cool flame to the slow combustion re- gion. In other words, this dumping corresponds to the transition from the CF to the SR zone of the explosion diagram of Fig. 5, at atmospheric pressure. The predicted peak temperature of the cool flame is ∼830 K for the 0.91-mm droplet with a succes- sive maximum temperature of ∼770 K, while the largest droplet reaches only 810 K during the first cool flame with a successive, more stable temperature of ∼780 K. The higher temperatures reached for the larger droplets explain the higher vaporization rates under these conditions. Experimental and predicted burning rates relat- ing to the autoignition of n-decane droplets in air at 1093 K are compared in Fig. 10. Reasonable agree- ment is also observed under these conditions, in which a hot ignition (HI) with flame temperatures higher than 2000 K is obtained. Therefore, these re- sults remain almost unchanged if a reduced kinetic scheme, containing only the high-temperature oxi- dation mechanism, is adopted. A systematic initial overprediction of burning rates can be observed in Fig. 10 for the initial time period, while the oppo- and A. Cuoci et al. / Combustion Fig. 11. Sensitivity analysis of droplet diameter squared to different thermal conductivities in the liquid phase. Comparison with experimental data [4] (n-decane droplets, P = 0.1 MPa, TG = 1093 K). site happens at longer times. In order to explain these deviations, which are more evident for the larger droplets, the model’s sensitivity to the thermal con- ductivity of the liquid was tested. As already noted, temperature differences on its surface induce internal motion inside the fuel droplet. A rigorous evaluation of the effect of these convective flows is quite com- plex and would make the whole problem difficult to handle. For this reason, a simple sensitivity analysis is presented in Fig. 11, where the burning rates are evaluated with different values of the liquid’s thermal conductivity. The increase in thermal conductivity re- sults in a more elevated initial heating of the droplet with a successive slight increase of the burning rates, in accordance with experimental measurements. This analysis indicates that an “effective” conductivity can help in reaching better agreement. Both experimental measurements and modeling predictions show higher burning rates for larger droplets, under these condi- tions too. Burning rates under spark-ignition conditions are more difficult to simulate, not only because the spher- ical symmetry is slightly lost but also due to the com- plications in the proper modeling of igniters and spark conditions. As previously noted when analyzing the boundary and initial conditions of the mathematical model, a specific approach is applied based on a short period of pure evaporation at ambient temperature, followed by a non-uniform radial profile of the tem- perature. Fig. 12 shows predicted and experimental burning rates of ignited n-decane droplets at 298 K in air. Once again reasonable agreement is observed between experimental and predicted burning rates, al- ways with a systematic deviation at least partially Flame 143 (2005) 211–226 221 Fig. 12. Droplet diameter squared for different initial droplet diameters during spark ignition experiments. Comparison between experimental [4] (points) and predicted (lines) data (n-decane droplets, P = 0.1 MPa, TG = 298 K, td is the time between starting and ignition). caused by the increased thermal conductivity of the liquid. The time evolutions of the diameter and tem- perature of the three different flames are reported in Fig. 13, together with experimental measurements of flame diameters. The model overpredicts these flame diameters by ∼20%. A sensitivity analysis shows that significant reduction in flame diameters could be ob- tained by decreasing the gas’s thermal conductivity, while radiative heat losses play only a very limited role. Flame temperatures, however, are significantly affected by the radiative heat loss and lower tempera- tures are reached by the larger droplets. For this rea- son, burning rates decrease when the initial droplet diameter is increased, as is clearly demonstrated in Fig. 14, which also shows that the model really is able to correctly predict the inverse influence of the initial diameter on burning rates in cold and hot ambi- ences, in line with experimental measurements [4]. In fact, the burning rate increases with the initial diam- eter in autoignition tests, whereas it decreases in the spark-ignited ones in cold ambiences. As already dis- cussed [4], the difference between the heat gain by the droplet and the heat loss by the flame is responsible for these effects. Finally, the internal circulation in a droplet affects the burning velocity, especially during an autoignition at high temperatures. Predicted burn- ing velocities with fourfold increases in liquid thermal conductivity are also shown in Fig. 14; they agree closely with the experiments, consistently maintain- ing the correct trend. The importance of radiative heat exchange and of radiative extinction diameter is ana- lyzed in the next paragraph. and 222 A. Cuoci et al. / Combustion Fig. 13. Comparison between experimental [4] (points) and predicted (lines) data of the flame diameter transients for different initial droplet diameters (a); predicted maximum gas-phase temperature for different initial droplet diame- ters (b) (n-decane droplets, P = 0.1 MPa, TG = 633 K, td is the time between starting and ignition). 5.4. Burning rates and extinction diameter The experimental measurements of Nayagam et al. [1] constitute a further interesting test of the model and help in the evaluation of burning rates at high temperatures. These measurements refer to n-heptane droplets in different atmospheres of oxygen + he- lium, nominally at 300 K and a total pressure of 0.1 MPa. Initial droplet diameters ranged from ∼2 to 4.1 mm. As already observed [1], the total rate of radiative energy loss from the flame is propor- tional to the third power of the flame’s diameter, since the gas is nearly transparent. For this reason, larger droplets show a larger extinction diameter, while smaller droplets cannot reach flame diameters large enough for extinction to occur. These data can Flame 143 (2005) 211–226 Fig. 14. Burning rates at different ambience tempera- tures. Comparison between experimental [4] (points) and predicted (lines) data. Dashed line, apparent liquid ther- mal conductivity four times higher (n-decane droplets, d0 = 1.54 mm, P = 0.1 MPa). be used to confirm different extinction diameters by varying the initial droplet diameter. The emissivities of H2O, CO, and CO2 contribute greatly to radiative heat transfer and consequently these species affect ra- diative loss rates and flame extinction. The maximum flame temperature is reduced by this radiant loss with a partial reduction in soot formation. Consequently, radiative energy loss from soot should be less impor- tant in flame extinction [1]. As already discussed, in line with Marchese et al.’s instructions [2], simula- tions were performed assuming preliminary evapora- tion of the n-heptane into the ambient atmosphere for 1 s, followed by the imposition of an initial temper- ature profile with a maximum of 2000 K in the gas phase, for about 4 droplet diameters. Predicted flame characteristics, burning rate, and extinction phenom- ena are only marginally affected by these assumptions regarding initial ignition conditions. As confirmed by the authors, the initial large expansion of a fuel droplet is due to gas from the chamber’s atmosphere initially present in the liquid fuel. For this reason, comparisons with experimental measurements are al- ways shifted by ∼0.25 s. Fig. 15 compares the ex- perimental and predicted burning rates for the spark ignition of a n-heptane droplet with an initial diameter of 2.86 mm in a mixture of 25 vol% oxygen in helium. Fig. 15b compares flame diameters and again shows the model overpredicting by ∼20%. The model seems able to predict the extinction diameter. The artificial reduction of this heat loss, of course, increases the flame temperature and reduces the extinction droplet diameter without significant changes in flame diame- ter. and A. Cuoci et al. / Combustion Fig. 15. Sensitivity analysis of droplet diameter squared experiments (a) and flame diameter (b) to different radia- tive losses during autoignition in a 25% O2–75% He mix- ture. Comparison between experimental [1] (points) and predicted (lines) data (n-heptane droplet, d0 = 2.86 mm, P = 0.1 MPa, TG = 293 K). The initial diameter of a droplet influences the flame’s diameter and its maximum value, but it has only a small effect on the time where the maximum temperature drops and approaches the liquid inter- face. The flame’s diameter is significantly affected by the thermal conductivity of the gas. For this reason, the same calculations are reported in Fig. 16 by re- ducing the thermal conductivity of the gas to 90 and 80% of its original value. Better agreement is ob- tained for the flame’s diameter, without significant variation in the initial burning rate. The extinction diameter decreases with a reduction in the gas’s con- ductivity. Fig. 16 also shows the predicted maximum temperature profiles in the gas phase. These clearly indicate the transition between the high- and the low- temperature mechanisms, with the attainment of an Flame 143 (2005) 211–226 223 Fig. 16. Sensitivity analysis of droplet diameter squared (a), flame diameter (b), and maximum gas phase temperature (c) to different gas thermal conductivity during autoignition ex- periments in a 25% O2–75% He mixture. Comparison be- tween experimental [1] (points) and predicted (lines) data (n-heptane droplet, d0 = 2.86 mm, P = 0.1 MPa, TG = 293 K). oxidation–evaporation phase at ∼700 K. The flame’s diameter, which is where the gas phase temperature is a maximum, moves rapidly toward the droplet dur- and 224 A. Cuoci et al. / Combustion Fig. 17. Droplet diameter squared for different initial droplet diameters during experiments with induced ignition. Com- parison between experimental (points) [1] and predicted (lines) data (n-heptane droplets, P = 0.1 MPa, Ta = 300 K, ambience: 25 vol% oxygen in helium). ing this second phase. Experimental confirmation of this behavior does not seem to have been reported in the literature. Nevertheless, the important role of low temperature mechanisms is strongly supported by the experiments of Xu et al. [4], as previously discussed. Similar results in terms of flame diameters and extinc- tion diameters have also been obtained for droplets with different initial diameters (see Fig. 17); they con- firm that the model is able to catch some features, but further work is necessary for a clearer understanding of the complex behavior of these systems. 6. Conclusions and future developments A time-dependent, spherically symmetrical model of an isolated droplet in microgravity, coupled with a detailed hydrocarbon oxidation mechanism, has been developed. The results compare reasonably with ex- periments across a broad range of combustion condi- tions. The overall kinetic scheme, including the low- temperature oxidation mechanism, allows the correct estimation of both the initial and the total ignition times, as well as the identification of the typical phe- nomena of aliphatic hydrocarbon oxidation. Multi- stage ignitions and cool flames are the results of com- petition between the low- and high-temperature ox- idation mechanisms. Modeling predictions correctly reproduce the transition between these two chemical regimes and are able to properly predict the differ- ent ignition regions. Moreover, the low-temperature mechanism plays a relevant role in governing the va- porization rate at low temperatures, whereas in the high-temperature experiments, radiative heat transfer Flame 143 (2005) 211–226 mainly controls the burning rate and the extinction phenomenon. The model properly predicts the influ- ence of droplet diameters both in the low-temperature and in the hot ignition experiments. The simulations confirm the opposite result of the effect of the droplet size on the burning rate, as summarized in Fig. 14. When the flame temperature is relatively low, the heat losses to the surrounding are not relevant. At high flame temperatures, on the other hand, the radiative losses are more relevant for larger diameters and the burning rate decreases with an increase in the initial diameter of the droplet. To the best of our knowledge, these comparisons are the first modeling confirmation of the opposite effect of droplet diameter on burning rates in cold and hot ambiences as experimentally ob- served by Xu et al. [4]. As far as the radiative terms are concerned, it should be noted that the extinction diameters are very sensitive to the estimates of modeling parameters. In this light, the assumptions of ignoring the contribu- tion of soot and/or other species should be recon- sidered. In fact, a similar agreement with the exper- imental measurements of extinction diameter was not obtained; the system is very sensitive to the accuracy with which the radiative flux is evaluated and further research is required to completely validate this model across the whole range of experimental conditions. This detailed kinetic model of the evaporation and combustion of fuel droplets has the intrinsic poten- tial of offering a characterization of the combustion environment, in terms of NOx, aldehydes, butadiene, PAH, and soot formation. Preliminary results on the effect of pressure and oxygen concentration on the formation of soot precursors, i.e., PAH, from ethanol droplets are encouraging. The extension to real and complex fuels, such as diesel or jet fuels, whose ki- netic schemes are already available, makes this model a powerful tool for the investigation and optimization of the combustion of liquids. Acknowledgments This work was partially supported by ASI (Agen- zia Spaziale Italiana). The work of Andrea Castellano, Simone Corno, Tito Paratico, and Michele Perego is also gratefully acknowledged. The authors thank the referees for their useful and detailed comments and suggestions. Appendix A. Numerical methods for DAE systems and the BzzMath library Before addressing differential-algebraic equations (DAE), it is worth spending some time on the numeri- cal algorithms for the solution of ordinary differential and A. Cuoci et al. / Combustion equation (ODE) systems. As the ancestors of DAE packages, ODE solvers share several peculiarities and a common theoretical background with them. Start- ing from DIFSUB [27] and passing to GEAR [28], LSODE [29], VODE [30,31], and BzzOde [32], the improvement in terms of the features and capabilities of ODE solvers has been continuous, mainly in terms of robustness and efficiency. As far as the DAE sys- tems are concerned, besides the well-known LSODI routine [29], the most recent decades were character- ized by the evolution of the DASSL routine [33] into DASPK [34] and by the introduction of BzzDae [35]. BzzMath library [21–23] comprises the aforemen- tioned BzzOde and BzzDae C++ classes. BzzMath is freeware for noncommercial use and can be down- loaded directly at http://www.chem.polimi.it/homes/ gbuzzi. While a more complete discussion is reported elsewhere, a rather concise description of the main features of BzzDae solver is reported in the follow- ing. BzzDae integrates DAE systems in the form{ y′ = f1(y, t), f2(y, t) = 0. On the robustness side, BzzDae is characterized by the following features: 1. According to Brenan et al. [36] and Brown et al. [34], BzzDae normalizes the algebraic portion of the Jacobian matrix through a simple a pri- ori division by h. With reference to matrix G, the BzzDae formulation exploits the following structure, G = [ I − hr0 ∂f1∂y1 hr0 ∂f1∂y2 ∂f2 ∂y1 ∂f2 ∂y2 ] , where y1 and y2 are the differential and algebraic variables with h being the stepsize and r0 the first coefficient of the BDF method [27]. This simple normalization significantly improves the robust- ness and precision of the solver. 2. If the model variables are physically bounded, the solver deals with constraints. DASPK, a com- mon Fortran routine, allows the user to assign a nonnegative scalar constraint to the solution vec- tor y throughout the integration path. This op- tion has been extended in BzzDae, where the user may simply assign maximum and/or min- imum constraint vectors. The solver automati- cally handles the constraints, taking care not to violate the assigned bounds. The control is per- formed before passing any illegal values to the DAE system routine. The correction vector, b, is accepted only when the nonlinear system, result- ing from the DAE problem, is accurately solved Flame 143 (2005) 211–226 225 to the assigned precision and when y simultane- ously complies with the constraints. As a result, the DAE function is always computed with safe y values and math errors are avoided a priori. 3. As suggested by Brenan et al. [36] the order is re- duced when the elements of the Nordsieck vector are not decreasing. 4. Both the order and the step size are reduced when there are convergence problems. 5. The integration order is automatically reduced to one and BzzDae restarts from the last suc- cessful convergence point when repeated conver- gence failures occur [32,35]. The following features contribute to the efficiency of BzzDae: 1. Droplet combustion, similar to most of the DAE chemical problems, is characterized by a relevant number of equations and by the need to numeri- cally evaluate the Jacobian matrix, J. Therefore, function evaluations have the greatest impact on CPU time [37]. The Jacobian evaluation becomes more and more exacting when the equation num- ber increases. BzzDae uses a distinct memory allocation for the Jacobian matrix and its factor- ization G = (I − hr0J). A direct consequence is an overall improvement in efficiency since when either a different step size or a method order is chosen, there is no need to reevaluate the Jaco- bian matrix and then superimpose its factoriza- tion, which is required by the nonlinear system solver based on the Newton method. 2. BzzDae uses the following criterion for updat- ing the Jacobian matrix. It checks whether J should be updated through the following equa- tion, y′ n+1 ∼= y′n + J · (yn+1 − yn)+ ft · (tn+1 − tn), where y′n+1, y′n, yn+1, yn are the variables at the n and n+1 iterations and ft is the time deriv- ative of the DAE system. With a large number of equations, the Jacobian numerical evaluation is rather time consuming, so it is advisable to de- lay the update of J as far as possible. Conversely, if the system has few equations, it is convenient to evaluate J more frequently in order to increase the efficiency of the Newton method. On this ba- sis, BzzDae updates the Jacobian matrix with a proper frequency dependent on J dimension. 3. DAE systems characterized by sparse and not necessarily structured Jacobian matrices can be easily solved by exploiting the automatic mem- ory allocation and matrix rearrangement of the C++ classes. Namely, structured systems with a tridiagonal block structure are efficiently solved by BzzDaeBlockTridiagonal. and 226 A. Cuoci et al. / Combustion The aforementioned features improve the overall performance of the BzzDae solver not only when dealing with constrained integration variables, but also under the following conditions: • Highly oscillating problems, that is, Jacobian matrix of the DAE system with complex λ eigen- values having a negative real part and a large imaginary part (Re(λ) < 0 and |Im(λ)| � 1). Ex- amples are reported in Figs. 4, 6, and 9. • DAE systems with large discontinuities in the derivatives, i.e., discontinuous or piecewise ini- tial conditions or physical properties, IF . . . THEN structures, code branching (see, for in- stance, Fig. 1). References [1] V. Nayagam, J.B. Haggard Jr., R.O. Colantonio, A.J. Marchese, F.L. Dryer, F.A. Williams, AIAA J. 36 (1998) 1369–1378. [2] A.J. Marchese, F.L. Dryer, V. Nayagam, Combust. Flame 116 (1999) 432–459. [3] A.N. Hayhurst, R.M. Nedderman, Chem. Eng. Ed. 21 (1987) 126–152. [4] G. Xu, M. Ikegami, S. Honma, K. Ikeda, X. Ma, H. Nagaishi, D.L. Dietrich, P.M. Struk, Int. J. Heat Mass Transfer 46 (2003) 1155–1169. [5] M. Tanabe, M. Kono, J. Sato, J. Konig, C. Eigenbrod, F. Dinkelacker, H.J. Rath, Combust. Sci. Technol. 108 (1995) 103–119. [6] M. Tanabe, T. Bolik, C. Eigenbrod, H.J. Rath, Proc. Combust. Inst. 26 (1996) 1637–1643. [7] S.K. Aggarwal, Prog. Energy Combust. Sci. 24 (1998) 565–600. [8] A.J. Marchese, F.L. Dryer, Combust. Flame 105 (1996) 104–122. [9] S. Schnaubelt, O. Moriue, T. Coordes, C. Eigenbrod, H.J. Rath, Proc. Combust. Inst. 28 (2000). [10] J.R. Yang, S.C. Wong, Combust. Flame 132 (2003) 475–491. [11] M. Perego, S. Corno, G. Buzzi-Ferraris, T. Faravelli, D. Manca, E. Ranzi, T. Parra, in: Proceedings of 6th Inter- national Conference on Engines for Automobile, Capri, Italy, 2003, SAE_NA Technical Paper Series 2003-01- 25. [12] T.E. Daubert, R.P. Danner, Data Compilation Tables of Properties of Pure Compounds, Design Institute for Physical Property Data, American Institute of Chemi- cal Engineering, 1985. [13] R.C. Reid, J.M. Prausnitz, B.E. Poling, The Properties of Gases and Liquids, McGraw–Hill, New York, 1987. [14] R.H. Perry, D.W. Green, Perry’s Chemical Engineers’ Handbook, seventh ed., McGraw–Hill, New York, 1998. Flame 143 (2005) 211–226 [15] R.J. Kee, P.M. Rupley, J.A. Miller, The Chemkin Ther- modynamic Data Base, Report SAND86-8215B, San- dia National Laboratories, 1987. [16] H. Wang, M. Frenklach, Combust. Flame 96 (1994) 163–170. [17] S.W. Benson, Thermochemical Kinetics, second ed., Wiley, New York, 1976. [18] A. Kazakov, J. Conley, F.L. Dryer, Combust. Flame 134 (2003) 301–314. [19] H. Zhang, M.F. Modest, J. Quantitative Spectrosc. Ra- diat. Transfer 73 (6) (2002) 649–653. [20] L.F. Shampine, in: R.C. Aiken (Ed.), What is Stiffness? Stiff Computation, Oxford Univ. Press, New York, 1985. [21] D. Manca, G. Buzzi-Ferraris, in: Proceedings of ICheaP 7, Taormina, Italy, 15–18 May 2005. [22] G. Buzzi-Ferraris, BzzMath: Numerical Libraries in C++, Politecnico di Milano; available at http:// www.chem.polimi.it/homes/gbuzzi. [23] G. Buzzi-Ferraris, Scientific C++: Building Numeri- cal Libraries the Object-Oriented Way, Addison Wes- ley/Longman, New York, 1993. [24] E. Ranzi, P. Gaffuri, T. Faravelli, P. Dagaut, Combust. Flame 103 (1995) 91–106. [25] P. Gaffuri, T. Faravelli, E. Ranzi, N. Cernansky, D. Miller, A. D’Anna, A. Ciajolo, AIChE J. 43 (1997) 1278–1286. [26] O. Moriue, C. Eigenbrod, H.J. Rath, J. Sato, K. Okai, M. Tsue, M. Kono, Proc. Combust. Inst. 28 (2000) 969–975. [27] C.W. Gear, Commun. Am. Chem. Soc. 14 (3) (1971) 185–190. [28] A.C. Hindmarsh, GEAR: Ordinary Differential Equa- tion System Solver, Report UCID 30001, Rev. 3, Lawrence Livermore Laboratory, Livermore, CA, 1974. [29] A.C. Hindmarsh, Am. Chem. Soc. SIGNUM News- lett. 15 (1980) 10–11. [30] P.N. Brown, G.D. Byrne, A.C. Hindmarsh, SIAM J. Sci. Stat. Comput. 10 (1989) 1038–1051. [31] G.D. Byrne, A.M. Dean, Comput. Chem. 17 (1993) 297–302. [32] G. Buzzi-Ferraris, D. Manca, Comput. Chem. Eng. 22 (11) (1998) 1595–1621. [33] L.R. Petzold, in: R.S. Stepleman, et al. (Eds.), Scientific Computing, North-Holland, Amsterdam, 1983. [34] P.N. Brown, A.C. Hindmarsh, L.R. Petzold, A Descrip- tion of DASPK: A Solver for Large-Scale Differential- Algebraic Systems, Lawrence Livermore National Re- port UCRL, 1992. [35] D. Manca, G. Buzzi-Ferraris, T. Faravelli, E. Ranzi, Combust. Theory Model. 5 (2001) 185–199. [36] K.E. Brenan, S.L. Campbell, L.R. Petzold, Numeri- cal Solution of Initial-Value Problems in Differential- Algebraic Equations, North-Holland, New York, 1989. [37] D. Manca, T. Faravelli, G. Pennati, G. Buzzi-Ferraris, E. Ranzi, Numerical Integration of Large Kinetic Sys- tems, ICheaP-2 Conf. Ser., ERIS, 1995, pp. 115–121. Autoignition and burning rates of fuel droplets under microgravity Introduction Equations for the model and numerical methods Liquid phase Gas phase Boundary conditions Initial conditions (t = 0) Numerical methods The spatial grid The structure of the DAE system The kinetic model Results and discussion Modeling predictions Explosion regions and ignition delay times Influence of droplet diameter on burning rates Burning rates and extinction diameter Conclusions and future developments Acknowledgments Numerical methods for DAE systems and the BzzMath library References