The development of computational models for the numerical simulation of chemically reacting ows operating in the turbulent regime requires the solution of partial di erential equations that represent the bal- ance of mass, linear momentum, chemical species, and energy. The chemical reactions of the model may in- volve detailed reaction mechanisms for the description of the physicochemical phenomena. One of the biggest challenges is the sti ness of the numerical simulation of these models and the nonlinear nature of species rate of reaction. This work presents a study of in situ adaptive tabulation (ISAT) technique, focusing on the accuracy, eciency, and memory usage in the simulation of homo- geneous stirred reactor models using simple and com- plex reaction mechanisms. The combustion of carbon monoxide with oxygen and methane with air mixtures are considered, using detailed reaction mechanisms with 4 and 53 species, 3 and 325 reactions respectively. The results of these simulations indicate that the developed implementation of ISAT technique has a absolute global error smaller than 1%. Moreover, ISAT technique pro- vides gains, in terms of computational time, of up to 80% when compared to the direct integration of the full chemical kinetics. However, in terms of memory usage the present implementation of ISAT technique is found to be excessively demanding. Keywords thermochemistry reduction in situ adaptive tabulation stirred reactor simulation  detailed reaction mechanism A. Cunha Jr (corresponding author) L. F. Figueira da Silva Department of Mechanical Engineering Pontif cia Universidade Cat olica do Rio de Janeiro Rua Marqu^ es de S~ ao Vicente, 225 - G avea Rio de Janeiro, 22451-900, Brasil E-mail: americo.cunhajr@gmail.com1 Introduction Combustion is a physicochemical phenomena charac- terized by a sequence of chemical reactions (most of then exothermic) which converts fuel/oxidizer mixtures (fresh gases) into combustion products (burned gases). Due to the large amount of energy released in these exothermic reactions, the combustion process is widely used in the operation of industrial devices such as gas turbines, process furnaces, etc. Design and optimization of devices which operate based on combustion processes is an important engi- neering task, which demands the development of com- putational models as predictive tools. Such computa- tional models may require the solution of partial dif- ferential equations that represent the balance of mass, linear momentum, chemical species, and energy. These models may include a detailed kinetic mechanism for the description of physicochemical phenomena involved [28, 6, 22, 35]. Typically, such reaction mechanisms for mixtures of hydrocarbons with air involve tenths of species, hundreds of elementary reactions, and timescales that vary up to nine orders of magnitude [36, 20, 21]. Concerning the numerical simulation of these mod- els, the challenge is related to reaction rate of the chem- ical species, which is nonlinear and presents a strong de- pendence with reaction mechanism dimension. There- fore, numerical simulation of a model with a detailed re- action mechanism is computationally demanding, which justi es the development of techniques to reduce the as- sociated computational cost. Several techniques to reduce the computational cost of these models are available in the literature. Essen- tially, they may be split into two classes. The rst class includes techniques that seek to reduce the dimension of reaction mechanisms, such as Quasi-Steady State Ap- 2 Americo Cunha Jr, Lu s Fernando Figueira da Silva proximation [8], Rate-Controlled Constrained Equilib- rium [13], Computational Singular Perturbation [15], the reduction to Skeleton Mechanisms [16], Intrinsic Low-Dimensional Manifold [19], Proper Orthogonal De- composition [11], and Invariant Constrained Equilib- rium Edge Pre-image Curve [26]. The second class in- cludes techniques that seek ecient algorithms to solve complex models, such as Look-Up Table [3], Repro- Modeling [32], Piece-Wise Reusable Implementation of Solution Mapping [31], Arti cial Neural Network [12], andIn Situ Adaptive Tabulation (ISAT) [24, 29, 17, 18]. A characterization of one type of transient homoge- neous reactor, using simple and complex thermochem- istry, is presented in this this work. The analysis ex- plores the e ects of mixing to residence time scale ratios as reactor controlling parameter. The ISAT technique is employed to reduce the computational time spent to approximate the solution of the equations that governs the evolution of reactive mixture in this transient reac- tor. This work also investigates accuracy, performance, and memory usage of the present ISAT technique imple- mentation. The accuracy analysis uses local and global metrics to quantify the error incurred by ISAT com- pared with a reference solution and investigates the sta- tistical seed in uence on the results. The performance is investigated by comparing the computational time gain obtained by ISAT with respect to a standard numeri- cal integration procedure. The memory usage analysis quanti es the amount of memory used by ISAT imple- mentation. This work is organized as follows: the transient reac- tor of interest, its governing equation as well the numer- ical scheme used are presented in the second section. The third section presents the basic theory of ISAT technique. The fourth section discusses about the re- actor con gurations studied and ISAT implementation issues. Finally, in the fth section, the main conclu- sions are summarized and a path for a future work is suggested. 2 Transient homogeneous reactor model 2.1 The geometry of reactive systems In the framework of the transported probability den- sity function (PDF) models [23], a reactive ow may be described by an ensemble of stochastic particles, which mimic the behavior of the uid system. Consider a transient spatially homogeneous reactive mixture evolving adiabatically and at constant pres- sure in a continuous ow reactor. The thermodynam- ical state of a uid particle in this reactor may be completely determined by the mass fraction Yi, wherei= 1;;ns, of thenschemical species and the speci c enthalpyh, which can be lumped in the composition vector de ned as = (h;Y 1;;Yns)T; (1) where the superscriptTdenotes the transposition oper- ation. One should note that, due to the invariance of the system number of atoms, which ensures the total conservation of the mass, the components of vector are not linearly independent. The composition vector of each particle in this ow reactor evolves according to d dt=( ;t) +S( ;t); (2) wheretdenotes the time, is the rate of change due to mixing, andSis the rate of change associated with the chemical reactions. Integrating Eq.(2) from an initial timet0to a timetone obtain the reaction mapping R( 0;t) = 0Zt t0( ;t0)dt0+Zt t0S( ;t0)dt0;(3) which is the solution of Eq.(2) starting from the initial composition (t0) = 0. The reaction mapping corre- sponds to a trajectory in composition space , which, for large values of t, tends to the equilibrium composition for the given enthalpy and mass fractions on 0. The composition space is the ( ns+1){dimensional Euclidean space where the rst direction is associated with the en- thalpy and the remaining nsare related to the chemical species. 2.2 Pairwise mixing stirred reactor The classical Partially Stirred Reactor (PaSR), used in [4, 27], describes by the interaction by exchange with the mean (IEM) micromixing model [7], but, for the purpose of testing a thermochemistry reduction tech- nique, it is desirable to employ a mixing model that leads to a composition region accessed during the so- lution process which is \wider" than that provided by the IEM model. A modi ed version of PaSR model called Pairwise Mixing Stirred Reactor (PMSR), [24], is designed to yield a much larger accessed region, and, hence, should provide a stringent test on the ability of ISAT technique to yield a reduction in computational time. In the PMSR model the reactor consists of an even numbernpof particles, initially arranged in pairs ( j1;j2) Assessment of a transient homogeneous reactor through in situ adaptive tabulation 3 such that the particles (1 ;2), (3;4),, (np1;np) are partners. Given a time step, t, for each discrete time kt, wherekis an integer, the model is characterized by three types of events: in ow ,out ow , and pairing . The in ow and out ow events consist of randomly se- lectingnin=ceil1 2np(t= r) pairs of particles, being rthe residence time within the reactor, and exchang- ing their thermodynamical properties by the proper- ties of a prescribed in ow. The pairing event consists of randomly selecting for pairing a number of pairs of particles, di erent from the in ow particles, equal to npair=ceil1 2np(t= p) , being pthe pairwise time. Then the chosen particles (in ow/out ow and paring) are randomly shued. Between these discrete times, the pairs of particles ( j1;j2) evolve according to d dt(j1) = (j1) (j2) m+S( (j1);t); (4) d dt(j2) = (j2) (j1) m+S( (j2);t): (5) where mis the mixing time. 2.3 Numerical Integration An operator splitting technique [37] is employed to nu- merically integrate Eq.(2). The overall process of inte- gration via operator splitting technique can be repre- sented as (t)mixing! mix(t+t)reaction! (t+t); (6) where given an initial composition 0and a time step t, the rst fractional step integrates the pure mixing system, d dt=( ;t); (7) to obtain (j) mix(t+t). Then, the pure chemical reac- tion system, d dt=S( ;t): (8) is integrated from an initial composition mix(t+t) over a time step tand gives (t+t). The operator splitting technique allows to solve each term in the evolution equation, Eq.(2), separately, using speci c ecient numerical methods to treat the par- ticular features inherent to the physical phenomenon modeled by each term [7].3 In situ adaptive tabulation 3.1 Linearized reaction mapping Consider a composition and an initial composition 0, so that Taylor expansion of the reaction mapping of around 0is R( ;t) =R( 0;t) +A( 0;t) +O(k k2); (9) where =  0,A( 0;t) is an n matrix, called mapping gradient matrix , with components given by Aij( 0;t) =@Ri @ 0j( 0;t); (10) kkdenotes the Euclidean norm of a vector, and O(k k2) denotes terms that have order k k2. The linear approximation Rl( ;t), obtained by ne- glecting the high order terms of Eq.(9), is second-order accurate at a connected region of composition space centered at 0. The shape of this region is unknown before the calculations, but ISAT algorithm approxi- mates this region by a hyper-ellipsoid, as will be shown in section 3.2. The local error of this linear approxima- tion is de ned as the Euclidean norm of the di erence between the reaction mapping at and the linear ap- proximation for it around 0, "= R( ;t)Rl( ;t) : (11) 3.2 Ellipsoid of accuracy The accuracy of the linear approximation at 0is con- trolled only if the local error is smaller than a positive error tolerance "tol, which is heuristically chosen. The region of accuracy is de ned as the connected region of the composition space centered at 0where local error is not greater than "tol. As shown in [24], this region is approximated by a hyper-ellipsoid centered at 0which is dubbed ellipsoid of accuracy (EOA), and is mathe- matically represented by the following inequation  TLLT "2 tol; (12) where the EOA Cholesky matrix, denoted by L, is lower triangular [10]. The adaptive step of ISAT algorithm involves the solution of the following geometric problem: given a hyper-ellipsoid centered at 0and a query composi- tion, q, outside it, determine a new hyper-ellipsoid of 4 Americo Cunha Jr, Lu s Fernando Figueira da Silva minimum hyper-volume, centered at 0, which encloses both the original hyper-ellipsoid and the point q. The solution of this problem is presented by [25] and is not shown here for sake of brevity. 3.3 Tabulation procedure Initially ISAT algorithm receives the time step tand the tolerance "tol. Then, in every time step, ISAT algo- rithm receives a query composition qand returns an approximation for the corresponding reaction mapping R( q;t). This approximation is obtained via numeri- cal integration of Eq.(8) or by the linear approximation Rl( q;t). During the reactive ow calculation, the computed values are sequentially stored in a table for future use. This process is known as in situ tabulation. The ISAT table, which is created by the tabulation process, in- cludes the initial composition 0, the reaction mapping R( 0;t), and the mapping gradient matrix A( 0;t). Using these elements it is possible to construct the lin- ear approximation. As the calculation proceeds, a new query composition, q, is received by ISAT and the ta- ble is transversed until a 0is found that is close to q. Depending on the accuracy, the linear approxima- tion around 0is returned or the reaction mapping of qis obtained by direct integration of Eq.(8). The ISAT table is a binary search tree, since this data structure allows for searching an information in O(log2ntab) operations, where ntabis the total num- ber entries in the tree, if the tree is balanced [14]. The binary search tree is basically formed by two types of elements, nodes, and leaves. Each leaf of the tree stores the following data: { 0: initial composition; {R( 0;t): reaction mapping at 0; {A( 0;t): mapping gradient matrix at 0; {L: EOA Cholesky matrix. Each node of the binary search tree has an associated cutting plane . This plane is de ned by a normal vector v= q 0; (13) and a scalar a=vT q+ 0 2 ; (14) such that all compositions withvT >aare located to the \right" of the cutting plane, as sketched in Fig- ure 1. The cutting plane construction de nes a search criterion in the binary search tree. φ0 φqv vTφ<a vTφ=avTφ>aFigure 1 Cutting plane in relation to EOA position. If during the calculation a query point qis encoun- tered within the region of accuracy, i.e., ""tol, but outside the estimate of EOA, then the EOA growth pro- ceeds as detailed in [25]. The rst three items stored in the binary search tree leaf 0;R( 0;t) andA( 0;t) are computed once, whereas Lchanges whenever the EOA is grown. Once a query composition qis received by ISAT table, the binary search tree is initialized as a single leaf and the exact value of the reaction mapping is returned. The subsequent steps are: 1. Given a query composition the tree is transversed until a leaf ( 0) is found. 2. Equation (12) is used to determine if qis inside EOA or not. 3. If qis inside EOA, the reaction mapping is given by the linear approximation. This is the rst of four outcomes, called retrieve . 4. If qis outside EOA, direct integration is used to compute the reaction mapping, and the local error is measured by Eq.(11). 5. If the local error is smaller than tolerance, "tol, the EOA is grown according to the procedure presented in [25] and the reaction mapping is returned. This outcome is called growth . 6. If local error is greater than the tolerance "toland the maximum number of entries in the binary search tree is not reached, a new record is stored in the bi- nary search tree based on qand the reaction map- ping is returned. The original leaf is replaced by a node with the left leaf representing the old composi- tion 0and the right leaf the new one qas shown in Figure 2. This outcome is an addition . 7. If the local error is greater than the tolerance "tol and the maximum number of entries in the binary search tree is reacted, the reaction mapping is re- turned. This outcome is called direct evaluation . Assessment of a transient homogeneous reactor through in situ adaptive tabulation 5 φ0v, a φ0 φq Figure 2 Binary search tree before and after the addition of a new node. 4 Results and discussion This section presents the results of numerical simu- lations performed in order to assess characteristics of PMSR reactor. This study investigates three con gu- rations of a PMSR reactor, two mixtures of CO=O 2, and one mixture of CH 4/air, using di erent time scales to de ne their behavior. Accuracy, performance, and memory usage of ISAT technique implementation are also investigated. For this purpose, the approximation for solution of the model equation, Eq.(2), obtained by ISAT technique are compared with a reference solution from direct integration (DI) procedure, described in sec- tion 2.3. 4.1CO=O 2mixture The rst two cases studied consider a 1024 particles PMSR, initially lled with a fuel-lean mixture of CO=O 2 at 2948:5 K and 1 atm. At every time step, a fuel-lean mixture of CO=O 2enters the reactor at Tin= 300 K andpin= 1 atm. Both fuel-leam mixtures have equiv- alence ratio equal to 0 :7. The constant pressure and enthalpy equilibrium state associated to the in ow mix- ture is reached at Teq= 2948:5 K. The reaction mech- anism used to describe COwithO2kinetics involves 4 species and 3 reactions, [9]. For the rst case a time scales con guration with m= r= p= r= 1=2 is used, so that the pairwise/mixing time scales are of the same order of magnitude as the residence time. In this situa- tion, partially stirred reactor (PaSR) conditions result within the reactor. For the second case, time scales con- guration adopted is m= r= p= r= 1=10, so that the pairwise/mixing time scales are small when compared to the residence time. Thus, the reactor should behave almost as a perfect stirred reactor (PSR), where the processes of mixing and pairing occur instantaneously. These cases use a binary search tree with a maximum of 50;000 entries, t= 10s, and"tol= 103. 4.2CH 4/air mixture The third case studied consists of a PMSR, with 1024 particles, initially lled with the combustion productsof a stoichiometric mixture of CH 4/air at 2100 K and 1 atm. At every time step, a stoichiometric mixture ofCH 4/air enters the reactor at Tin= 300 K and pin= 1 atm. The constant pressure and enthalpy equi- librium state associated to this mixture is reached at Teq= 2225:5 K. The reaction of CH 4with air kinetics is described by GRI mechanism version 3.0 [30], with 53 species and 325 reactions. For this PMSR one have m= r= p= r= 1=4, so that it is expected to behave like a PaSR. This simulation uses a binary search tree with a maximum of 60,000 entries, t= 0:1 ms and "tol= 103. 4.3 Analysis of reactor behavior Figures 3 and 4 show the comparison between DI and ISAT calculations, as function of the dimensionless time =t= r, for the ensemble average hTiand the en- semble variance T02 of the reduced temperature in the cases 1 and 2, where the reduced temperature is de ned as T=TTin TeqTin; (15) and the ensemble average and the ensemble variance operators are respectively de ned as h i=1 npnpX j=1 (j)andD 02E =D 2E h i2; (16) being a generic property of the reactive system. For case 1 results, which span over a range of 500 residence times, one can observe an excellent qualita- tive agreement for hTiand T02 . The ensemble av- erage value rapidly drops from the initial value to reach the statistically steady-state regime around hTi= 0:3. The ensemble variance rapidly grows, then decreases until it reaches the statistically steady-state value around T02 = 0:13. The analysis of the two gures shows that the statistically steady-state regime is reached af- ter = 10. In case 2, where a range of 250 residence times have been computed, one can also observe an excellent quali- tative agreement for hTiand T02 . Again, the overall history of the PMSR is the same for DI and ISAT. Simi- lar results, not shown here for sake of brevity, have been obtained for the other thermochemical properties of the reactors. The di erence among both cases is due to the be- havior of each reactor at the statistically steady state 6 Americo Cunha Jr, Lu s Fernando Figueira da Silva 0 100 200 300 400 50000.20.40.60.811.2 τ∗/angbracketleftT/angbracketright∗ DIISAT (a) Evolution of hTifor case 1. 0 50 100 150 200 25000.20.40.60.811.2 τ∗/angbracketleftT/angbracketright∗ DIISAT (b) Evolution of hTifor case 2. Figure 3 Comparison between DI (|) and ISAT (- - -) calculations of the ensemble average of the reduced temperature. 0 100 200 300 400 50000.050.10.150.2 τ∗/angbracketleftT/prime2/angbracketright∗ DIISAT (a) Evolution of T02 for case 1. 0 50 100 150 200 25000.050.10.150.2 τ∗/angbracketleftT/prime2/angbracketright∗ DIISAT (b) Evolution of T02 for case 2. Figure 4 Comparison between DI (|) and ISAT (- - -) calculations of the ensemble variance of the reduced temperature. regime. Indeed, the case 2 reactor behavior is governed by a competition between the chemical and residence times mostly. Therefore, the thermodynamical proper- ties steady-state PDF is concentrated over a smaller range than in case 1, where the mixing and pairing time scales are large when compared to the residence time. This behavior is illustrated in Figures 5 and 6, which present the comparison between DI and ISAT computations of mean histograms, averaged over the last 50 residence times, of TandYOfor cases 1 and 2. These gures underscore the in uence of PMSR con- trolling parameters, i.e., the time scales ratios m= r and p= r, on the thermochemical conditions prevailing within each reactor. Indeed, the temperature within the case 2 reactor is such that almost only burned gases are found. On the other hand, case 1 reactor is character- ized by a bimodal temperature distribution with a large probability of nding T= 0:04 and a broader tempera- ture distribution leaning toward the burned gases. Such a distribution, as it could be expected, is re ected on YO histogram, which also exhibits a bimodal distribution. The comparison between DI and ISAT results for the rst two statistical moments of TandYOH, for case 3, are presented in Figures 7 and 8. One can notereasonable and good agreements for the Tand theYOH statistics, respectively. The results of case 3 show a large discrepancy for the Tthan that obtained in case 1. The reaction mechanism of methane is much more complex than the mechanism used to model the carbon monox- ide chemical kinetics in case 1, which would in princi- ple lead to the PMSR with methane to assume a wider range of possible thermodynamic states. Thus, it could be expected that the present binary search tree, with 60k entries (almost similar to that used in case 1), would yield comparatively lower accuracy. Finally, Figure 9 presents the comparison between DI and ISAT computations of the mean histograms, av- eraged over the last 50 residence times, of the Tand YHCO for case 3. The temperature histogram presents a bimodal distribution with large probability of nd- ingT0:9 and a broader temperature distribution leaning to the fresh gases. On the other hand, YHCO histogram shows a distribution essentially concentrated in the fresh gases region, and nearly homogeneous else- where. This behavior illustrates the fact that HCO , an intermediate species, appears in low concentration in the burned gases. Assessment of a transient homogeneous reactor through in situ adaptive tabulation 7 0 0.2 0.4 0.6 0.8 1 1.200.20.40.60.8 T∗Relative Frequency DIISAT (a) Histogram of Tfor case 1. 0 0.2 0.4 0.6 0.8 1 1.200.20.40.60.8 T∗Relative Frequency DIISAT (b) Histogram of Tfor case 2. Figure 5 Comparison between DI and ISAT calculations of the mean histograms (averaged over the last 50 residence times) of the reduced temperature. 0 0.01 0.02 0.03 0.04 0.05 0.0600.20.40.60.8 YORelative Frequency DIISAT (a) Histogram of YOfor case 1. 0 0.01 0.02 0.03 0.04 0.05 0.0600.20.40.60.8 YORelative Frequency DIISAT (b) Histogram of YOfor case 2. Figure 6 Comparison between DI and ISAT calculations of the mean histograms (averaged over the last 50 residence times) of theOmass fraction. 0 50 100 15000.20.40.60.811.2 τ∗/angbracketleftT/angbracketright∗ DIISAT (a) Evolution of hTifor case 3. 0 50 100 15000.0050.010.0150.02 τ∗/angbracketleftYOH/angbracketright DIISAT (b) Evolution of hYOHifor case 3. Figure 7 Comparison between DI (|) and ISAT (- - -) calculations of the ensemble average of the reduced temperature and OHmass fraction. 4.4 Analysis of ISAT accuracy The relative local error of , which is a time dependent function, is de ned as "r; (t) = DI(t) ISAT (t) DI(t) ; (17)where subscripts DIandISAT denote DI and ISAT cal- culations of , respectively. This metric represents a local measure of error incurred by ISAT. The evolution of the relative local errors of Tand YOstatistical moments are presented in Figures 10 and 11. Concerning case 1 errors, one can observe a large statis- tical variation due to the stochastic nature of the PMSR 8 Americo Cunha Jr, Lu s Fernando Figueira da Silva 0 50 100 15000.050.10.150.2 τ∗/angbracketleftT/prime2/angbracketright∗ DIISAT (a) Evolution of T02 for case 3. 0 50 100 15001234x 10−5 τ∗/angbracketleftYOH/prime2/angbracketright DIISAT (b) Evolution of Y02 OH for case 3. Figure 8 Comparison between DI (|) and ISAT (- - -) calculations of the ensemble variance of the reduced temperature and OHmass fraction. 0 0.2 0.4 0.6 0.8 1 1.200.20.40.60.8 T∗Probability Density Function DIISAT (a) Histogram of Tfor case 3. 0 1 2 3 x 10−400.20.40.60.8 YHCOProbability Density Function DIISAT (b) Histogram of YHCO for case 3. Figure 9 Comparison between DI and ISAT computations of the mean histograms (averaged over the last 50 residence times) of the reduced temperature and HCO mass fraction. model, with amplitudes reaching 13%. In case 2, rela- tive errors of the order of 1% only can be observed. All realizations of the stochastic processes in DI and ISAT calculations shown above use the same statisti- cal seed for pseudorandom number generator. However, since the mixing model is statistical in nature, it could be expected that the seed value may in uence the ISAT behavior. Therefore, if the DI calculation seed is kept xed and ISAT seed is changed, the results for hTiand "r;hTiare expected to be modi ed. This can be seen in Figure 12, where one may observe, by comparison with Figures 3 and 10, an increase in the "r;hTifor both cases. The analysis of the graphs in Figure 12 indicates that 1024 particles are not sucient to guarantee the statistical independence of the computed results. This hypothesis is also con rmed if one considers Figure 13, where it is possible to see discrepancies in DI and ISAT mean histograms of Tfor both cases, even if, quali- tatively these histograms are not much di erent than those shown in Figure 5. Thus, one can conclude thatalthough the results not exhibit the statistical indepen- dence with a sample of 1024 particles, they do not vary much with a sample of this size. Possibly a sample of 4096 particles is sucient to ensure the independence of the results. The mean relative error of over an interval  , de ned as "r; M=1  Zt+ t"r; (t0)dt0; (18) is a metric that gives a global measure of the error in- curred by ISAT. The mean relative errors associated to the statistical moments of the Tand theYOare presented in Table 1. One can observe that, when the same seed is used for ISAT and DI, the mean error of the average properties is smaller than 1.2% in case 1, and of 0.1% in case 2 only. ForTandYOensemble variances, the mean error is large in case 1. Concerning the simulations where the used seeds are di erent, the mean error of the proper- Assessment of a transient homogeneous reactor through in situ adaptive tabulation 9 0 100 200 300 400 50000.050.10.150.2 τ∗/epsilon1r,/angbracketleftT/angbracketright∗ (a) Evolution of "r;hTifor case 1 0 50 100 150 200 25000.050.10.150.2 τ∗/epsilon1r,/angbracketleftT/angbracketright∗ (b) Evolution of "r;hTifor case 2 0 100 200 300 400 50000.050.10.150.2 τ∗/epsilon1r,/angbracketleftYO/angbracketright (c) Evolution of "r;hYOifor case 1 0 50 100 150 200 25000.050.10.150.2 τ∗/epsilon1r,/angbracketleftYO/angbracketright (d) Evolution of "r;hYOifor case 2 Figure 10 Evolution of the relative local error for the ensemble average of the reduced temperature and of the Omass fraction, using the same statistical seeds in both cases. ties average/variance are an order of magnitude greater than the previous one. Table 1 Mean relative error of for cases 1 and 2. Equal Seeds Di erent Seeds Case 1 Case 2 Case 1 Case 2 hTi1:0 % 0:0 % 8.1 % 5.3 % hYOi 1:2 % 0:1 % 11.7 % 1.3 % T02 1:0 % 0:5 % 6.6 % 10.1 % Y02 O 0:9 % 0:5 % 7.0 % 6.1 % In the early development of ISAT technique [24] it was noted that the choice of the tolerance could a ect the accuracy of the problem solution. In order to inves- tigate the e ect of the tolerance on the present results, Figure 14 presents the absolute global error ,"g, as a function of ISAT error tolerance for cases 1 and 2. The absolute global error is de ned, over a time interval  , as "g=1  Zt+ t h i(t0)DIh i(t0)ISAT dt0;(19)whereh idenotes a vector which the components are the ensemble averages of components. As may be seen in Figure 14, when the same statisti- cal seed is used, "gdecreases linearly as "tolis reduced, for both studied cases. However, if di erent statistical seeds are used, one can note a limit behavior where "g does not decrease if "tolis reduced. This saturation in "gvalue (102) indicates that a decrease in "tolvalue below 103is not e ective in cases where the statistical seeds are di erent. The evolution of the relative local error of the statis- tical moments of TandYOH, for case 3, are presented in Figure 15. Compared with the results of case 1, the errors of case 3 present more oscillations. Furthermore, for this case, "g= 2:4103is an order of magnitude larger than the values obtained in cases 1 and 2, also using"tol= 103. Since di erent seeds are used for DI and ISAT calculations in case 3, large errors in the results are to be expected. 4.5 Analysis of ISAT performance The comparison of time evolution of ISAT outputs (num- ber of additions, growths, retrieves, and direct evalua- tions) and of binary search tree height, as well as the 10 Americo Cunha Jr, Lu s Fernando Figueira da Silva 0 100 200 300 400 50000.050.10.150.2 τ∗/epsilon1r,/angbracketleftT/prime2/angbracketright∗ (a) Evolution of "r;hT02ifor case 1 0 50 100 150 200 25000.050.10.150.2 τ∗/epsilon1r,/angbracketleftT/prime2/angbracketright∗ (b) Evolution of "r;hT02ifor case 2 0 100 200 300 400 50000.050.10.150.2 τ∗/epsilon1r,/angbracketleftYO/prime2/angbracketright (c) Evolution of " r;D Y02 OEfor case 1 0 50 100 150 200 25000.050.10.150.2 τ∗/epsilon1r,/angbracketleftYO/prime2/angbracketright (d) Evolution of " r;D Y02 OEfor case 2 Figure 11 Evolution of the relative local error for the ensemble variance of the reduced temperature and of the Omass fraction, using the same statistical seeds in both cases. corresponding rates of change, for cases 1 and 2, are presented in Figures 16 and 17. A rst important ob- servation is that the number of additions in both cases reaches the maximum prescribed value in the binary search tree of 50k. As a consequence of binary search tree saturation, the additions curve reaches a steady state after 5 :4 and 1:2 residence times ( r), in the rst and second cases, respectively. This indicates that ISAT table is saturated earlier when mixing is faster ( m= r small). Figures 16 and 17 also show the evolution of binary tree height, which reaches steady state 1 :8 rin the rst case and 0:4 rin the second case. It is also noteworthy that, in both cases, the tree height is an order of magni- tude smaller than the total number of entries in the tree (17k in case 1 and 7k in case 2). This di erence between height and total entries in data table ensures the eciency of search process for a new query, which could be performed up to three and seven times faster in cases 1 and 2, respectively, than a vector search. In the rst case, Figures 16 and 17 show that the number of growths presents a sharp rate of change around r, whereas, in the second case, this occurs around 5 r. In both cases, growth steady state occurs after 10 r.During both simulations the number of growths is al- ways smaller than additions. This indicates that the desirable massive increase of the ellipsoids of accuracy, to improve the estimate of the region of accuracy, is not observed. This behavior might circumstantial to the re- action system of the carbon monoxide considered, since, due to its simplicity (only three reactions), a small part of the realizable region should be assessed during the course of the calculation. The Figure 16 shows that, after tree saturation oc- curs, the number of retrieves, and direct evaluations exceed the number of additions in both cases. In case 1 there is a higher occurrence of retrieve events, whereas in case 2 direct evaluations prevail. The number of re- trieves exhibits a linear limit behavior in both cases. The ISAT behavior for the second case re ects the fact that the binary tree of this case is poor, i.e., contains too few compositions in the region accessed by the calcula- tion. As a consequence, the number of direct evaluations vastly outnumbers ISAT operations. This explains the error behavior observed in section 4.4. When ISAT technique is more ecient than DI pro- cedure, the computational time spent by ISAT is smaller than computational time spent by DI. The computa- Assessment of a transient homogeneous reactor through in situ adaptive tabulation 11 0 100 200 300 400 50000.20.40.60.811.2 τ∗/angbracketleftT/angbracketright∗ DIISAT (a) Evolution of hTifor case 1. 0 50 100 150 200 25000.20.40.60.811.2 τ∗/angbracketleftT/angbracketright∗ DIISAT (b) Evolution of hTifor case 2. 0 100 200 300 400 50000.10.20.30.4 τ∗/epsilon1r,/angbracketleftT/angbracketright∗ (c) Evolution of "r;hTifor case 1 0 50 100 150 200 25000.10.20.30.4 τ∗/epsilon1r,/angbracketleftT/angbracketright∗ (d) Evolution of "r;hTifor case 2 Figure 12 Comparison between DI (|) and ISAT (- - -) calculations, using di erent statistical seeds, of the ensemble average of the reduced temperature and the corresponding relative local errors. 0 0.2 0.4 0.6 0.8 1 1.200.20.40.60.8 T∗Relative Frequency DIISAT (a) Histogram of Tfor case 1. 0 0.2 0.4 0.6 0.8 1 1.200.20.40.60.8 T∗Relative Frequency DIISAT (b) Histogram of Tfor case 2. Figure 13 Comparison between DI and ISAT calculations, using di erent statistical seeds, of the mean histograms (averaged over the last 50 residence times) of the reduced temperature. tional time spent by ISAT is the sum of the computa- tional time spent at each of its possible outputs. There- fore, the eciency condition can be stated as nA A+nG G+nR R+nDE DE<nDI DI; (20) wherenA,nG,nR,nDE,nDIare the number of addi- tions, growths, retrieves, direct evaluations, and direct integrations, respectively, whereas A, G, R, DE, and DIare the corresponding average duration at each op- eration. Assuming that G DI1; DE DI1; and R DI 1: (21) it is possible to show [5] that a necessary (but not suf- cient) condition for a calculation using ISAT to be faster than the same calculation using DI is 12 Americo Cunha Jr, Lu s Fernando Figueira da Silva 10−810−610−410−210010−810−610−410−2100 /epsilon1tol/epsilon1g dif. seedsame seed (a)"gvs"tolfor case 1. 10−810−610−410−210010−810−610−410−2100 /epsilon1tol/epsilon1g dif. seedsame seed (b)"gvs"tolfor case 2. Figure 14 Absolute global errors as function of the error tolerance, using a binary search tree with 50k entries. 0 50 100 15000.10.20.30.4 τ∗/epsilon1r,/angbracketleftT/angbracketright∗ (a) Evolution of "r;hTifor case 3. 0 50 100 15000.10.20.30.4 τ∗/epsilon1r,/angbracketleftYOH/angbracketright (b) Evolution of "r;hYOHifor case 3. 0 50 100 15000.10.20.30.4 τ∗/epsilon1r,/angbracketleftT/prime2/angbracketright∗ (c) Evolution of "r;hT02ifor case 3. 0 50 100 15000.10.20.30.4 τ∗/epsilon1r,/angbracketleftYOH/prime2/angbracketright (d) Evolution of " r;D Y02 OHEfor case 3. Figure 15 Evolution of the relative local error of the ensemble average mean and the ensemble average variance of the reduced temperature and OHmass fraction. nR nA> A DI1: (22) ForCO=O 2mixtures estimates lead to A= DI10, [5]. It is possible to see in Figure 16 that, for case 1, nR approximately exceeds nAby a factor of 130, whereas in case 2 this factor is only 6. Therefore, as A= DI<10 in case 2, ISAT calculations are not expected to be faster than DI procedure.As can be seen in Table 2, where a comparison of computational time is shown, cases 1 and 2, for "tol= 103are computed using DI in 4 :0 ks and 2:0 ks respec- tively, whereas with the use of ISAT, the same cases spent 2:3 ks and 2:1 ks, respectively. Speed-up factors of 1:7 (case 1) and 1 :0 (case 2) are obtained, where thespeed-up factor is de ned as the ratio between the computational time spent by DI and the computational time spent by ISAT. This table also allows to com- pare the computational time spent by DI and ISAT Assessment of a transient homogeneous reactor through in situ adaptive tabulation 13 10−2100102104100102104106108 τ∗ISAT outputs AdditionGrowthRetrieveD.Eval.Height (a) Case 1. 10−2100102104100102104106108 τ∗ISAT outputs AdditionGrowthRetrieveD.Eval.Height (b) Case 2. Figure 16 Evolution of ISAT algorithm outputs and of the height of ISAT binary search tree. 10−210010210400.511.52x 104 τ∗ISAT outputs rate of change AdditionGrowthRetrieveD.Eval.Height (a) Case 1. 10−210010210400.511.522.5x 104 τ∗ISAT outputs rate of change AdditionGrowthRetrieveD.Eval.Height (b) Case 2. Figure 17 Evolution of the rates of change of each ISAT algorithm outputs and of the height of ISAT binary search tree. for di erent values of error tolerance. An increase in processing time is obtained as ISAT error tolerance is reduced, which is to be expected, given the fact that lower values of "tolcorrespond to a smaller region of accuracy. Indeed, as "tolis decreased, it is less likely that ISAT returns a retrieve, which is ISAT output with lower computational cost. In case 1, for all values of "tol, ISAT technique o ers an advantage, in terms of com- putational time, when compared to the process of di- rect integration, reducing on average the computational time in 42%. On the other hand, in case 2, no reduc- tion in computational time is seen. As discussed above, this behavior is natural, once the necessary condition for eciency of the algorithm is not reached. Aiming to analyze the asymptotic behavior of ISAT speed-up, case 1 has also been simulated for 50,000 time steps. This computation with ISAT spends 204 :7 ks, whereas if DI were to be used the same calculation would spend1024:0 ks (speed-up factor of 5.0). There- fore, in this asymptotic case, ISAT spends approximately 80% less time than DI. The pioneer work of Pope (1997) [24] reports an asymptotic speed-up factor of 1000, butthis impressive factor has not been observed in the study developed here. The evolution of ISAT algorithm outputs, ISAT bi- nary search tree height, and the corresponding rates of change, for case 3, are presented in Figure 18. As in cases 1 and 2 the number of additions reaches the max- imum allowed value in the binary search tree, 60k for this case. The additions curve reaches a steady state after = 3:2. The binary search tree maximum height for case 3 is7k as in case 2. Here the steady state of the tree height is observed at = 1:28, as can be seen in Figure 18. The behavior of the number of growths is quite similar to that of case 1, where the greater rate of change occurs near = 2 but, now, the steady state is reached before = 10. For case 3, nevertheless, the number of growths is not always smaller than the number of additions. Ini- tially the growths exceed additions by nearly an order of magnitude. After = 0:2 the additions exceeded the number of growths, only to be overcome again once steady state is reached. This behavior indicates that the ellipsoids of accuracy show a considerable period 14 Americo Cunha Jr, Lu s Fernando Figueira da Silva Table 2 Comparison between the computational time spent by DI and ISAT and the corresponding speed-up factors. Case 1 Case 2 "tol time spent (ks) speed-up time spent (ks) speed-up DI4:0 2:0 1021:92:12:11:0 1032:31:72:11:0 1042:41:72:20:9 1052:51:62:20:9 1062:61:52:20:9 10−2100102104100102104106108 τ∗ISAT outputs AdditionGrowthRetrieveD.Eval.Height (a) Outputs for case 3. 10−21001021040123456x 104 τ∗ISAT outputs rate of change AdditionGrowthRetrieveD.Eval.Height (b) Outputs rates of change for case 3. Figure 18 Evolution of ISAT algorithm outputs, height of ISAT binary search tree, and the corresponding rates of change. of adaptation, which allows to access, with precision, a larger portion of the realizable region. As in the previous cases, after tree saturation, the number of retrieves and direct evaluations exceed the number of additions. The number of direct evaluations overcomes total number of retrieves in approximately 60%, which is not negligible, but less than that observed in case 2. This indicates that the binary search tree covers a signi cant portion of the realizable region in the composition space. From Figure 18 one can estimate A= DI50 in order that ISAT technique be more ecient than DI. For case 3, the computational time spent by DI and ISAT are 689 ks and 455 ks, respectively. Therefore, a speed-up of 1.5 is observed. In this case, the e ect of other values of ISAT error tolerance has not been inves- tigated. This PMSR has 1024 particles and its evolution is computed during = 150. Thus, the system of gov- erning equations de ned by Eq.(8) is solved 7 ;680;000 times. In this case ISAT technique allows to save 34%, in terms of computational time. For problems that require solving Eq.(8) several times, one could speculate that ISAT technique would provide an even better perfor- mance improvement, since more retrieves are expected when compared with the more costly direct integration.4.6 Analysis of ISAT memory usage The studied cases 1 and 2 are both modeled by a reac- tion mechanism with 4 species and use a binary search tree with 50,000 entries for ISAT simulations. These parameters lead to a memory consumption by ISAT algorithm of approximately 40 Mbytes, which is very small when compared with the available memory in the workstation used. On the other hand, case 3 is modeled by a reac- tion mechanism with 53 species and use a binary search tree with 60,000 entries. This case uses approximately 3:3 Gbytes of memory. It is noteworthy that this value of memory usage is not negligible, when compared to the total memory available on the workstation. Thus, for practical purposes, this tree is the largest one that may be used to simulate a PMSR with this methane re- action mechanism and 1024 particles. This underscores what is perhaps the greatest shortcoming of ISAT tech- nique, its huge expense of memory. 5 Concluding remarks This work assessed characteristics of a PMSR reactor through the numerical simulation of reactor con gu- rations with simple ( CO=O 2) and complex ( CH 4/air) reaction mechanisms and also investigated ISAT tech- Assessment of a transient homogeneous reactor through in situ adaptive tabulation 15 nique as an option to evaluate a computational model with detailed thermochemistry. For this purpose, the accuracy, performance, and memory usage of the cor- responding technique were analyzed through the com- parison of ISAT results with a reference result, obtained from direct integration of PMSR model equation. The studied cases analysis showed that ISAT technique has a good accuracy from a global point of view. Also, it was possible to identify the statistical seed e ect on the control of absolute global error, which decreases mono- tonically as ISAT error tolerance was reduced, when the same seed is used for ISAT and DI calculations. On the other hand, when di erent seeds were used, a limit value for ISAT error tolerance was observed. In terms of performance, ISAT technique allows to reduce the com- putational time of the simulations in all cases studied. ForCO=O 2cases, a speed-up of up to 5.0 was achieved, whereas for CH 4/air, the algorithm allowed to save 34% in terms of computational time. Moreover, ISAT tech- nique presented the desired feature of speed-up factor increase with the complexity of analyzed system. Re- garding the memory usage, ISAT technique showed to be very demanding. This work will be extended by cou- pling ISAT technique to the hybrid LES/PDF model by [2, 34, 1, 33] for description of turbulent combus- tion, when a detailed reaction mechanism would allow a better description of combustion. In this context, ISAT could be a viable option that may be able to decrease, to an acceptable level, the simulation time. Acknowledgements The authors acknowledge the support awarded to this research from CNPq, FAPERJ, ANP and Brazilian Combustion Network. The authors are indebted to Professor Guenther Carlos Krieger Filho from Universidade de S~ ao Paulo (USP), for providing the code that served as example to the code developed in this work, and for his hos- pitality during the authors' visit to USP. This work was per- formed while the second author was on leave from Institut Pprime, CNRS, France.