M Lubrano Lavadera
email: m.lubranolavadera@irc.cnr.it

Y Song

P Sabia

O Herbinet

M Pelucchi

A Stagni

T Faravelli

F Battin-Leclerc

M De Joannon

Marco Lubrano Lavadera

Oscillatory behavior in methane combustion: on the influence of the operating parameters

Keywords: Temperature oscillations, Jet Stirred Flow Reactor, chemical kinetics, diluted combustion

The influence of the main process parameters on the oscillatory behavior of methane oxidation was analyzed in conditions relevant for Low Temperature Combustion (LTC) processes. The investigation was performed by means of direct comparisons between experimental measurements realized in two Jet Stirred Flow Reactors used at atmospheric pressure. Coupling the operating conditions of the two systems, wide ranges of inlet temperature (790 -1225 K), equivalence ratio (0.5 < Φ < 1.5), methane mole fraction (XCH4 from 0.01 to 0.05), bath gases (i.e. He, N2, CO2, or H2O) as well different overall mixture dilution levels, were exploited in relation to the identification of oscillatory regimes.

Although the reference systems mainly differ in thermal conditions (i.e. heat exchange to the surroundings), temperature measurements suggested that the oscillatory phenomena occurred when systems working temperature accessed a well-identifiable temperature range.

Experimental results were simulated by means of a detailed kinetic scheme and commercial codes developed for complex chemistry processes. Simulations were also extended considering systems with different heat loss to the surroundings, thus passing from adiabatic to isothermal systems. Results highlighted the kinetic nature of the dynamic behavior.

Because of predictions were consistent with experimental tests, further numerical analyses were realized to identify the kinetics responsible for the establishment of oscillatory phenomena. Temperature oscillations were predicted for a significant reactor working temperature range, where oxidation and recombination kinetic routes, involving carbon C1-2 species as well as reactions of the H2/O2 sub-scheme, become competitive, thus boosting limit cycles behaviors. Oscillatory phenomena cease when the system working temperatures exceed characteristic threshold values with the promotion of faster oxidation routes that diminish the inhibiting effects of recombination reactions.

Introduction

Despite much progress in the application of new combustion concepts to practical systems, such as Low-Temperature Combustion (LTC) processes, there remain unresolved issues on the stabilization of the oxidation process 1 . A problem related to LTC processes is the susceptibility to thermo-kinetic oscillations, originating from non-linear dependence of relevant variables. Such instabilities may hamper the large-scale implementation of advanced combustion processes 2 . Therefore, there is a need for a fundamental understanding of critical parameters that promote instabilities.

To properly investigate oscillatory phenomena, open systems are necessary. A very useful tool for the observation of combustion regimes as a function of the operating parameters is the Well Stirred Reactor (WSR). In this device, in contrast to chemical reactions occurring in a closed vessel that proceed to a unique point of thermodynamic equilibrium, it is possible to achieve a multiplicity of permanent steady states far from equilibrium as well as sustained oscillations that may also be preserved indefinitely.

A large literature exists on the analysis of dynamic phenomena taking place in open wellstirred reactors. For example, such reactors have been widely used to investigate oscillatory cool flames and complex ignition phenomena under non-adiabatic conditions for high molecular weight organic compounds at low temperature (T < 850 K), as a result of thermal feedback coupled to complex kinetics [3][4][5][6] . Scientific interpretations of these phenomena date back several decades because modeling of the low temperature oxidation of hydrocarbons is relevant to the onset of "knock" in spark-ignition (SI) engines, or to the ignition in diesel fuels [7][8][9] . Recently, the same oscillatory behavior was also found for dimethyl ether (DME) between 530 and 650 K 10 .

Oscillatory behaviors have also been recognized for low molecular weight fuels, such as hydrogen [11][12][13][14][15][16][17] , or hydrogen and carbon monoxide mixtures [18][19][20] . Notwithstanding exothermic contributions, the oscillations of such mixtures originate from a kinetic interaction driven by the competitive elementary reactions H + O2 = OH + O and H + O2 + M = HO2 + M. For the latter, in fact, the rate of recombination waxes and wanes according to the third body efficiencies of species in the vessel.

Other experimental 2,21,22 and numerical [START_REF] Dally | 6 th Asia-Pacific Conference on Combustion[END_REF][START_REF] Wada | European Combustion Meeting[END_REF] studies have focused on methane ignition and oxidation under MILD combustion operative conditions in flow reactors at atmospheric pressure. These studies highlighted the onset of non-isothermal oscillatory regimes similar to cool-flame oscillations but at higher temperatures (1000 -1300 K). Sabia et al. [START_REF] Sabia | [END_REF] extended the studies by de Joannon et al. 2,21 adding H2 to CH4/O2/N2 mixtures. These studies showed that H2 can reduce the parameters space where oscillations were detected, with respect to methane mixtures. de Joannon et al. 2,21 and Sabia et al. 22,[START_REF] Sabia | [END_REF] identified different kinetic routes promoted by high inlet temperatures and highly diluted mixtures, with respect to conventional flames.

In particular, they pointed out that the competition between termolecular recombination reactions and branching reactions is crucial for the appearance of these regimes. Wada et al. [START_REF] Wada | European Combustion Meeting[END_REF] added that a key reaction pathway for the onset of oscillations is the route that converts CH3 to CH2O. At low temperature CH3 is oxidized to CH3O and CH3O2 that, in turns, react to form CH2O, while at intermediates one CH3 produces directly CH2O by reacting with HO2 radicals.

All these previous studies 2, [21][22][START_REF] Dally | 6 th Asia-Pacific Conference on Combustion[END_REF][START_REF] Wada | European Combustion Meeting[END_REF][START_REF] Sabia | [END_REF] claim that oscillations have a thermo-kinetic nature, thus that they originate from a balance between the heat released by chemical reactions and the heat loss to the surroundings.

On the other hand, Sabia et al. 26 have recently found that dynamic behaviors can be also numerically obtained under adiabatic conditions, so they highlighted the kinetic aspect of such phenomena.

Another process, which received considerable attention, is the thermal oxidative coupling of methane (OCM). It is generally accepted that a homogeneous coupling of methyl radical, producing ethane, becomes important when the temperature is between 1000 -1300 K 27 , which is the same temperature range in which non-isothermal oscillatory regimes have been highlighted for methane oxidation. Therefore, a detailed study of methane oscillations can be also utilized to improve OCM process.

Recently, Bagheri et al. 28 explored thermochemical oscillations of methane from a kinetic modeling perspective, also focusing on the effect of different dilution gases (N2, H2O, CO2) in promoting or inhibiting such instabilities. The authors found that the variations of heat capacity associated with the different bath gases directly impact the oscillations characteristics. On the other hand, the authors noted that in WSR isothermal studies available in the literature, oscillations were not reported; however, an attempt at modeling these cases revealed the presence of cyclic oscillatory behaviors.

Given this background, the objective of the current work is to extend the fundamentals of methane combustion stability under well-stirred flowing conditions, with particular emphasis on the physical and chemical factors that control the oscillatory phenomenon. In order to seek the means to interpret events in a wider range of external parameters, experimental tests were carried out in two Jet Stirred Flow Reactors (JSFRs), characterized by different heat release.

Oscillations onset, waveform, amplitude and frequency were characterized as a function of (1) inlet temperature, (2) equivalence ratio, (3) residence time, (4) methane mole fraction/mixture dilution level, and (5) the chemical nature of the bath gas (He, N2, CO2 or H2O), at atmospheric pressure.

The interest in the effects of the bath gas is twofold: first, for LTC systems, the combustion occurs in the presence of non-inert species (CO2 or H2O) due to the strong recirculation of exhausted gases; second, such species are commonly present in Low Calorific Value (LCV) fuels. These species may have a thermal effect (i.e. higher heat capacities) or a kinetic effect, participating directly in bimolecular elementary reactions, or as third body in termolecular reactions with high collisional Chaperon efficiencies. The comparison with mixtures diluted with inert species can reveal the role of CO2 or H2O in oxidation processes at lowintermediate temperatures.

Finally, the experimental measurements were supported by extensive numerical analyses of a chosen detailed thermo-kinetic model in order to identify the kinetics responsible for the onset and termination of oscillations with respect to operating parameters.

Experimental and numerical tools

Experimental apparatuses

In the present work, two spherical Jet Stirred Flow Reactors (JSFR) were used: the former is available at the Laboratoire Réactions et Génie des Procédés, Université de Lorraine, CNRS (Nancy) 29 , while the latter at the Institute for Research on Combustion, CNR (Napoli) 2 .

Hereafter they will be referred as the "Nancy" and the "Napoli" reactors. Photos of the two JSFRs are reported as Supporting Information. They were designed under the same construction rules 30,31 , but conceived to work under different operating conditions (i.e. residence time, mixture composition). The detailed description of the reactors is provided elsewhere 2,29 .

The opportunity to combine the experimental tests from the two facilities implies that a large field of parameters is accessible for an exhaustive investigation of methane oxidation process.

Both the JSFRs consist of a fused silica sphere where reactions take place. Quartz is commonly used to prevent catalytic reactions. The Napoli reactor consists of a sphere of 113 cm 3 while the Nancy one of 85 cm 3 . Fresh reactants enter the reactor through an injection cross (four nozzles, 1 mm diameter for the Napoli reactor; 0.3 mm for the Nancy one) located at the center of the sphere. The well-mixing condition was verified in the past by means of pulse tracer experiments [START_REF] Levenspiel | Chemical Reaction Engineering[END_REF] . The Napoli reactor behaves as an ideal WSR for residence times (τ) lower than 0.6 s, while the Nancy one lower than 10 s. Thus, the experiments were respectively carried out with residence times equal to 0.5 s and 2 s.

In both the reactors the total flow is composed by methane, oxygen and diluent. In the Napoli device, a primary flow, composed by oxygen and diluent is pre-heated through a quartz tube located within two semi-cylindrical electrically heated ceramic fiber insulated ovens supplied by Watlow Srl, regulated by PID controllers. The primary flow subsequently mixes with the secondary one, composed by fuel and diluent, in a premixing chamber before entering inside the reactor. The reactor is located within two other thermally insulated electrical fiber ovens to minimize heat loss to the surroundings. A recirculation air system provides a homogeneous temperature distribution in the oven. The homogeneity of the reactor external wall temperature is monitored by means of two movable N type thermocouples.

In the Nancy experiments, the reactor is preceded by an annular preheating zone in which the temperature of the gases is increased up to the reaction temperature before entering inside the reactor. Then the pre-mixed mixture enters the reactor. Both the spherical reactor and the annular preheating zone are heated by Thermocoax heating resistances rolled up around their walls. The external wall reactor temperature is monitored by means of a K type thermocouple.

In both experimental facilities, the preheating system allows reaching inlet temperatures (Tin) up to 1250 K. Tin is continuously monitored by a K type shielded thermocouple located in the center of the nozzles, with a precision of ± 2 K. Electronic control of the heating system gives uniform temperatures within ± 5 K over the external vessel surface and constant in time to within ± 0.5 K. This allows maintaining the reactor external wall at the same temperature as the inlet temperature.

The residence time of the fuel mixture in the preheating system is more than 100 times lower than the mean residence time in the reactor. Thus, oxidation of the fuel inside the injector is prevented. When a temperature increase is detected by the Tin thermocouple, indicating the onset of ignition/oxidation reactions during the mixing process, the experimental conditions are not considered. For diluted methane mixtures this condition does not occur in the investigated temperature range.

Gases are stored in cylinders and supplied to the system by means of calibrated digital thermal mass flow controllers supplied by BronkHorst High-Tech. They are preceded by filters to retain particulate impurities. The precision of the mass flow rates corresponds to ± 0.5% of set operating conditions.

Gases are provided with high purities: 99.999% for nitrogen, oxygen, and helium, 99.5% for carbon dioxide and methane.

The experiments carried out with steam diluted mixtures rely on a water feeding system consisting of a stainless steel cylinder, pressurized with dry nitrogen, a calibrated liquid mass flow meter and a Controlled Evaporator Mixer (CEM) provided by BronkHorst High-Tech.

The CEM consists of a control valve, a mixing device and a heat exchanger whose temperature is controlled by an integrated system controller. A carrier gas (nitrogen) is used to promote the evaporation process and to transport the vapor. The temperature of gas lines downstream the evaporator is controlled by electrical coils to prevent water condensation.

The reactor pressure is continuously monitored by means of a transducer connected to the outlet from the vessels. It is kept constant in both the apparatuses by means of a needle valve located on the gas exhaust line.

Experimental procedures and measurement methodology

In each test, the total mixture mass flow rate was set for the specified residence time τ, adjusted for the inlet temperature Tin, and the reactor pressure P. The mixture residence time is defined as τ = reactor volume/volume flow rate. For the Napoli tests, P was maintained at 1.1 atm and τ fixed at 0.5 s, whereas, for the Nancy ones, P = 1.05 atm and τ = 2 s.

For each test, vessels were flushed with inert gas before use. After fuel injection, a period of 5 minutes was allowed to reach thermal equilibrium. Initial test conditions refer to the parameters set prior to reactions taking place. Therefore, the inlet temperature has not been varied continuously, but each experimental point corresponds to a different test.

In order to accurately monitor the temperature changes during the oxidation process, homemade type R unshielded thermocouples gas welded from a very fine Pt-Pt 13% Rh wire (0.2 mm bead size) interfaced to a PC with a high-speed multi-channel module, supplied by National Instruments, were used. A photo of the thermocouple located inside the reactor is reported as Supporting Information. The global thermocouple response time is about 30 ms thus allowing for a faithful record of the reactor temperature during oscillatory behaviors with a precision of ± 2 K. A thin-wall alumina double-core tube covers thermocouples wires to prevent catalytic reactions. Herein, it has to be pointed out that thermocouples made using the same procedure were used in both the JSFRs.

The outlet concentration of the reactants was analyzed by gas chromatography. In the Napoli reactor, the exit gases were cooled down by means of a heat exchanger installed at the reactor outlet to quench the oxidation reactions. Then, gases passed through a silica gel trap to eliminate the moisture content before gas analyses. Therefore, results are reported on a dry basis. The outlet was connected to a gas chromatograph (Agilent 3000) equipped with a column Molsieve 5A and a pre-column Plot U and a Thermal Conductivity Detector (TCD).

In Nancy the quantification of the products from the reactor was realized with the aids of two gas chromatographs (GCs). The first GC equipped with a thermal conductivity detector was used to quantify O2. The second GC equipped with flame-ionization detector preceded by a methanizer and a PlotQ capillary column was used to quantify CH4, CO, CO2, C2H4, C2H6 and C2H2.

Maximum relative errors in mole fractions are estimated to be ± 10%, when the concentrations approach the detection threshold (about 10 ppm). The instrument calibration was performed daily using standard gas mixtures stored in cylinders.

Experiments were performed at least twice under all test conditions to ensure measurement reproducibility. The test repeatability was excellent under all conditions with differences among the various tests well within the reported experimental uncertainties.

Computational approach

Simulations were carried out using the OpenSMOKE++ 33 PSR code, calculating the transient solutions of an open WSR [START_REF] Cuoci | [END_REF] . Geometrical characteristics were the same as those of the experimental reactors, i.e. a volume of 113 cm 3 and an internal surface of 113 cm 2 for the Napoli reactor and a volume of 85 cm 3 and an internal surface of 93 cm 2 for the Nancy one.

Experimental mixture inlet conditions reported in Tab. 1 are declared in the input files. Since the two experimental systems are non-adiabatic, Newtonian-type heat losses were enabled.

The reactor overall heat transfer coefficient (U) in the JSFRs was calculated by means of an empirical correlation available in the literature 34 and fixed at 1.3 x 10 -3 cal/cm 2 s K. Reactor walls were assumed to be infinitely thin with negligible temperature gradients. The wall temperature was set equal to the inlet one.

A recently updated detailed gas-phase chemical kinetic mechanism "POLIMI" (C1C3LT_1412) was used. It was developed by Ranzi et al. 35 for the pyrolysis and combustion of a large variety of fuels. The mechanism consists of 2642 elementary reactions with 107 species. Thermodynamic properties are included in the database. A detailed analysis and validation of the reliability of this scheme is discussed elsewhere 35 . This kinetic mechanism was chosen as reference scheme because it was one of the most reliable in predicting oscillatory behaviors under MILD operating conditions, for both methane 2,24 and propane 36,37 mixtures.

Results

Experiments were carried out from fuel lean to fuel rich mixtures for a wide range of inlet temperatures, at atmospheric pressure. The bath gas was changed along with the fuel mole fraction and the overall molar fraction of diluent species (hereafter indicated as dilution "d").

The experimental conditions are summarized in Tab. 1. Data for cases S1, S2 and S3 are from Bagheri et al. 28 . Furthermore, the data collected for each dataset are provided in the Supporting Information. In order to classify the combustion regimes as a function of Tin, as exemplifying case, the reactor temperature (T), and the CH4 and O2 concentration, for a stoichiometric CH4/O2 mixture diluted in N2 at 90% (Napoli tests case S1) are reported in Fig. 1.

Case Φ τ (s) X CH4 X O2 X He X N2 X CO2 X
In this figure, vertical lines identify the transition conditions from one state to another in response to a stepwise change of Tin. At low temperatures (T < 1000 K), non-reactive conditions (I), characterized by zero temperature increase and zero fuel conversion, are identified. For increasing Tin, a secondary stationary state (II) with very small extents of oxidation is recognized for 1000 K < Tin < 1100 K with a maximum temperature increase of For Tin > 1170 K, the mixtures ignite and the system reaches a stationary state (IV) with temperature increase between 70 and 80 K and almost complete conversion. This steady state regime can be attained either monotonically or via damped oscillations.

The experimental condition reported in Fig. 1 is used as a reference case because it includes all the observed critical phenomena. However, the stationary state IV has been observed only for relatively short residence times (namely in the Napoli reactor). This can be explained considering the relative importance of sensible heat removal by system losses at the vessel walls and the heat convective transport associated with the outflow by calculating the ratio between two characteristic times, the mixture mean residence time τ and the Newtonian cooling time (tc). The latter depends on the vessel characteristic dimensions V/S (V and S are respectively the vessel volume and surface), and on the ratio ρcp/U (ρ and cp represent the mixture density and specific heat capacity), which does not vary significantly with the temperature. When an exothermic reaction is in progress, if τ > tc, heat losses to the walls predominate over heat transport in the outflow. Therefore, the longer the residence times, the lower the system temperature. For both reactors, tc can be assumed to be equal to 0.07 ± 0.01 s. Considering the residence times used in the two JSFRs, τ/tc is about 7 for the Napoli reactor and 30 for the Nancy one.

Thus, both systems are far from the adiabatic condition. In particular, in the Nancy configuration, temperature increases due to oxidation reactions are not measurable at stationary states because the system approaches the isothermal condition. However, interestingly, temperature oscillations were still observed.

The analysis of the complex waveform periodic events could elucidate the role of parameters in triggering oscillations, and hence it could provide some useful hints for the study of oscillatory phenomena in practical systems.

Based on the above classifications, Tin -Φ maps were constructed (Figs. 2 and3) considering all experimental conditions investigated in both reactors (see Tab. Fig. 2a summarizes the Napoli experimental tests carried out for mixtures of CH4/O2 diluted in N2 at 90% for several equivalence ratios (Φ = 0.5, 1, 1.5). For Φ = 0.5 (L1), temperature oscillations occur for Tin between 1055 K and 1100 K. For Φ = 1 (S1), the dynamic behavior starts for Tin > 1100 K and ends at Tin = 1180 K. Finally, for the fuel rich condition (R7), the dynamic behavior was not observed in the whole investigated Tin range.

An analogous identification of regimes can also be made when H2O or CO2 are used as diluent species.

Fig. 2b shows experimental results for a mixture of CH4/O2 diluted in N2 and H2O with a relative concentration of 55% and 45% and a global mixture dilution level of 90% (L3, S3

and R9). It can be seen that steam leads to a shift of the oscillations regions towards higher Tin with respect to the N2-diluted system. This shift is even more evident when mixtures are diluted in CO2 (Fig. 2c, L2, S2 and R8). In this case, temperature oscillations also appeared for the fuel rich condition. involve a wider range of Tin (Tin > 1100 K). For the fuel rich conditions (R1), temperature oscillatory behaviors were not detected in the whole temperature range considered.

As the methane mole fraction increases (Figs. 3c,3e), the oscillatory region on the map widens under stoichiometric conditions (S6 and S8), and shrinks for lean ones. In particular, at Φ = 0.75, this region covers a temperature range from 1100 K to 1125 K when XCH4 = 0.05 (L8) and vanishes when XCH4 = 0.025 (L6). Oscillations were never observed under rich conditions, for mixtures diluted in He. A similar identification of regimes can also be made when CO2 is used as bath gas (Figs. 3b, 3d, 3f). Oscillatory regions are clearly wider than in the He diluted cases. More specifically, in agreement with the data obtained at a lower residence time in the Napoli reactor, the oscillatory region extends up to fuel rich conditions, indicating that the nature of temperature oscillations with CO2 bath gas under fuel rich conditions is more complex, in agreement with the observations of Sabia et al. 26 The performance of the chosen detailed kinetic scheme to predict experimental data was tested by means of systematic runs. As it can be seen from the location of the grey rectangles on the maps (Figs. 2 and3), the adopted kinetic model is able to reproduce the experimental Tin -Φ with relatively minor differences on the location of boundaries between regions of different behavior. Accuracy checks of the numerical procedures ensured that the predicted oscillations are not due to numerical integration problems. For validation, the simulations were also performed using ChemKin PRO [START_REF]CHEMKIN-PRO 15131[END_REF] , yielding the same results.

In particular, the detailed kinetic scheme is able to reproduce quite accurately the influence of the equivalence ratio and of the bath gas on the oscillation region boundaries. Moreover, the kinetic model correctly predicts the absence of oscillations for fuel rich conditions when using N2, He, or H2O as bath gases. It has to be pointed out that in the case of Napoli experiments, the kinetic model predicts the occurrence of oscillations under fuel rich conditions when N2 and N2-H2O are used as bath gases, but the Tin range is very narrow (less than 5 K). However, disagreements exist in model predictions for several conditions. In particular, in Figs. 2c and 3f it is possible to observe that the POLIMI model predicts a discontinuity in the oscillation Tin range under fuel lean conditions. Such a behavior was not experimentally observed. Furthermore, in Figs. 3d and3e, for fuel lean conditions, the numerically predicted Tin range for temperature oscillations is wider compared to the experimental one, and, in Fig. 3c, the model predicts oscillations under fuel lean conditions that were not experimentally detected. In addition, in Fig. 3f oscillations are not predicted under rich conditions whereas they were detected experimentally.

In order to show some characteristic features of temperature oscillations as a function of Tin, four typical temporal temperature profiles are displayed in Figs. 4 and5.

Fig. 4 refers to a stoichiometric CH4/O2 diluted in N2 at 90% (S1) preheated to Tin = 1100 K (Fig. 4a) and Tin = 1165 K (Fig. 4b). The periodic profiles are reported after ignition, thus the initial value is arbitrary. In Fig. 4a, temperature increases slowly at first then abruptly up to a maximum value, followed by an instantaneous temperature decrease. The oscillation period is about 6.7 s with an amplitude of about 60 K. The minimum temperature is about 25 K higher than Tin. As Tin is increased (Fig. 4b), the temporal temperature profile has a different shape with a lower period (less than 1 s). These two profiles delineate a characteristic feature of the oscillatory behaviors: for increasing Tin, temporal profiles become more symmetric, while the oscillation amplitude decreases and the frequency increases. Furthermore, the difference between the minimum temperature and the inlet one also increases. These observations rule out the possibility of classifying the identified periodic phenomenon as multiple ignitions. The periodicity should be a feature of the kinetic mechanism itself and is not an artifact of flow operation. In Fig. 5a, the temporal profile is symmetric and the oscillation period is about 6 s with an amplitude of about 25 K. The minimum temperature is equal to the inlet one. As the equivalence ratio is increased (Fig. 5b), the temporal temperature profile has the same period, but a different shape. In particular, two maxima can be identified within each oscillation cycle with higher amplitude for the first peak. This confirms that the nature of temperature oscillations with CO2 under fuel rich conditions is more complex, in agreement with Sabia et al. 26 . Other examples of temporal temperature profiles are provided as Supporting Information.

The computed temporal temperature profiles are not reported in Figs. 4 and5 because, although the overall features of the computed temperature history resemble the experimental measurements very closely, the predicted temperature increase is substantially larger than the experimental (usually within 200 K). This occurs because the predicted temperature oscillations are characterized by spike variations in less than 10 ms and, even though the response time of the thermocouples is very short (≅ 30 ms), the temperature increase cannot be properly resolved in time. A detailed analysis of numerical temperature profiles can be found in Bagheri et al. 28 . For all bath gases, frequencies are relatively low at first. As Tin increases, oscillations frequencies increase. In both figures (Figs. 6a and6b), CO2 leads to a lower frequency with respect to other bath gases. The data related to N2-H2O dilution (Fig. 6a) lie between those obtained using N2 bath gas and those related to CO2 bath gas for Tin < 1140 K. For Tin > 1140 other two mixtures. Overall Fig. 6 suggests that the oscillation frequencies are significantly modified by the presence of different bath gases and they are strongly dependent on Tin. The POLIMI kinetic model satisfactorily predicts the experimental trends. However, depending on the parameters, a significantly longer period is predicted by the kinetic model at the highest Tin compared to experimental results, with a maximum deviation of 1.33 s.

Discussion

The factors that govern the stability of stationary states and the establishment of temperature oscillations need to be understood using numerical analysis. The reference model gives a satisfactory, qualitative and quantitative account for the Φ -Tin regimes along with oscillatory frequencies, as Tin and mixture composition are changed. Therefore, numerical analyses of complex dynamic behavior can be used to understand the role of elementary reactions in determining the observed behaviors.

It is well known that a central feature of oscillatory regimes is a negative feedback mechanism (i.e. the reaction products influence the rates of earlier reaction steps in the mechanism) that restricts the runaway acceleration of an intermediate concentration and differentiates the oscillatory reaction from an explosion 4 , where the feedback is always positive. Negative feedback can be caused by a number of different phenomena, some purely physical, others of chemical nature 4 . As reported in the introduction, oscillatory regimes may be significantly influenced by heat losses. Therefore, as a first step it is appropriate to analyze the influence of the heat loss on the observed phenomenon (Fig. 7).

Figure 7a shows the effect of modifications to the heat transfer coefficient (U) on the steadystate reaction temperatures as a function of Tin for a stoichiometric CH4/O2/N2 mixture. The different conditions analyzed allow to sensibly vary the system working temperature, thus representing an indirect way to parametrically evaluate the operating conditions that boost the onset of oscillatory behavior. In addition, the system was modeled under the asymptotic assumptions of adiabatic (Fig. 7b) and isothermal (Fig. 7c) conditions in order to highlight crucial aspects of the problem. Figure 7a shows the reactor temperature for a stoichiometric mixture diluted in N2 at 90% by varying the heat transfer coefficient (U) from 0.0013 to 0.13 cal/cm 2 sec K in steps of one order of magnitude. As expected, the temperature values are inversely proportional to the heat loss coefficients. For U = 0.0013 cal/cm 2 sec K, the system temperature remains equal to Tin up to Tin = 1080 K, then for 1080 < Tin < 1090 K, the temperature increases very slightly, and finally, at about Tin = 1090 K, the mixture ignites and temperature oscillations are predicted. The ignition condition is identified on the figure by an open circle. In Fig. 7a, the shadowed area under the Tmax curves identifies the maximum and minimum values of temperature oscillations. For these inlet conditions, the minimum reactor temperature during oscillations is equal to the inlet one. Oscillatory behavior is identifiable up to Tin = 1180 K (square symbol). For Tin > 1180 K, the system shows a steady stationary state behavior with a temperature increase of about 80 K, in agreement with the experimental results reported in Fig. 1.

For U = 0.013 cal/cm 2 sec K the behavior is very similar but the mixture ignition occurs at higher Tin (about 1120 K), then temperature oscillations establish. They have lower amplitudes with respect to the previous case. Temperature oscillations die out at Tin = 1250 K, when a steady state condition is reached. Similar observations apply for the case U = 0.13 cal/cm 2 sec K. Ignition occurs for Tin = 1130 K with the onset of temperature oscillations that cease at Tin = 1260 K. Interestingly, these oscillations are observed when the temperature decreases below a fixed value, i.e. T = 1260 K, reported in Fig. 7 with a horizontal dotdashed line, thus suggesting that kinetics is controlling the suppression of oscillations.

Fig. 7b shows the behavior of the ideal PSR under adiabatic conditions for a stoichiometric CH4/O2 mixture diluted in N2 at different levels (from 90% to 99.5%).

For N2 = 90%, the ignition occurs at about Tin = 1060 K, then the system reaches a steady stationary state. No oscillatory behavior is predicted. As the dilution level of the mixture increases (dN2 = 94% and 98%) the auto-ignition event occurs at higher Tin with lower temperature increase. None of the system exhibits a periodic dynamic behavior because of the larger temperature increase. In fact, when ignition occurs, the working temperature T is above 1260 K, thus no oscillations are observed, in agreement with Figure 7a. In the case of dN2 = 99.5 %, the ignition occurs at 1260 K with a temperature increase of 35 K. In this case the ignition occurs when Tin exceeds the identified threshold temperature, thus no oscillations are observed. However, for the latter system a further simulation has been considered. In particular, the auto-ignition event was forced by doping the inlet mixture with a low concentration of OH (200 ppm). For this case, the auto-ignition occurs for Tin = 1200 K and oscillations are observed. The periodic dynamic behavior can no longer be sustained when the system temperature exceeds 1260 K (horizontal dot-dashed line). It is worth noting that the addition of OH does not significantly change (1 K) the system temperature in the high temperature stationary solution.

On the other hand, transition from low-temperature steady states to periodic oscillatory states is attainable when the mixture residence time is higher than the mixture ignition delay time.

Since an increase of the heat transfer coefficient or of XN2 leads to an increase of the ignition delay time, a shift of the appearance of oscillations with Tin is observed in Figs. 7a and7b.

In addition, it was noted that an oscillatory dynamic regime is also observed under isothermal conditions, as shown in Fig. 7c. Although the temperature is not allowed to vary in this latter case (the shadow area is replaced by a dashed line), oscillations in methane mole fraction over time are highlighted in the insert of Fig. 7c, referring to inlet temperature conditions of 1200 K. Interestingly, methane concentration oscillations occur when Tin is in the temperature range 1120 -1260 K. These simulations showed some peculiarity of the oxidation process at intermediate temperature for stoichiometric CH4/O2 mixtures:

 the system is prone to exhibit periodic dynamic behavior when the working temperature is approximately in the range 1120 -1260 K, even under adiabatic or isothermal conditions.

 the high temperature steady state always occurs for T ≥ 1260 K.

Further simulations were performed for lean mixtures (Φ = 0.5) diluted in N2 at 90%. They confirmed the presence of an invariant temperature that marks the shift from oscillatory behaviors to high temperature steady states. This temperature differs from the one identified for a stoichiometric mixture, thus further highlighting the chemical nature of the oscillatory phenomenon.

As discussed in the introduction, oscillatory cool flames in LTC conditions arise from the interaction between heat release, complex chemical kinetics, and heat losses, therefore they could not be identified in adiabatic or isothermal conditions. Moreover, cool flame regions for higher hydrocarbons occur at lower temperatures, where peroxide species play the major As mentioned in the introduction to this study the oscillations observed in the oxidation of hydrogen, or of syngas mixtures, originate from a kinetic interaction between free radical chain branching and termination reactions (i.e. H + O2 = OH + O vs. H + O2 + M = HO2 + M). The last reaction is promoted by the high collisional efficiency of water formed during mixture oxidation. In particular, the POLIMI mechanism proposes a collisional efficiency for H2O 14 times higher than that of nitrogen. To decouple possible relevant effects of reactions belonging to the H2/O2 sub-system, further simulations were performed setting H2O and CO2 third body efficiencies equal to the nitrogen one. Interestingly, temperature oscillations were still predicted, thus suggesting that oscillations observed for methane mixtures arise from a different kinetic subset.

Rate of Reaction analyses (RR), flux diagrams, and first-order temperature sensitivity analyses were used to determine the rate limiting reactions for different Tin, considering a stoichiometric mixture diluted at 90% in N2. The numerical analysis has been performed for simplicity under isothermal conditions.

The most significant and representative results (RR analyses) are reported in Fig. 8, showing the rate of the controlling reactions as a function of the inlet temperature. The aim of Fig. 8 is not to show the whole kinetic routes that lead to methane conversion as a function of Tin, but to show how kinetics change from slow combustion to high temperature conditions. In previous papers 2,26 , the authors highlighted the strong competition between the methyl oxidation channels and methyl recombination to ethane, occurring at intermediate temperatures. Given this information, only the reaction rates involving methyl radicals are reported in Figure 8a.

The highlighted grey area represents the Tin range in which oscillations are predicted. The main reaction paths change substantially moving from the region of slow reaction to the oscillatory region and from the oscillatory region to the high temperature region.

At low temperatures (1000 K < Tin < 1100 K), methane is converted to methyl radicals by means of H-abstraction reactions (mostly by OH and HO2 not reported in Figure 8) which is subsequently oxidized to CH3O, formaldehyde, formyl radical and CO. Alternatively, two methyl radicals can either react to form ethane (CH3 + CH3 + M = C2H6 + M) terminating the radical chain or enhance reactivity through propagation reactions such as CH3 + HO2 = CH3O

+ OH. This latter reaction in particular, despite being formally a propagation step, converts two relatively unreactive radicals (CH3 and HO2) to very reactive radicals (CH3O and OH). For completeness, the same reaction rates are reported in Fig. 9 for one oscillation cycle at Tin = 1200 K, for the same reference mixture reported in Fig. 7c. Such operating conditions lead to the establishment of species concentration oscillations in time. Note that methane profile is also reported on the secondary axis.

It is possible to note that during an oscillation period, when methane concentration is relatively high, the methyl recombination to ethane and the oxidation to CH3O by HO2 radicals are the most important reactions. The reactions between CH3 and OH radicals (CH3 + OH + M = CH3OH + M, CH3 + OH + M = CH2(S) + H2O) become important only when methane is fully converted, due to higher availability of radicals. After methane reaches its minimum value, species accumulated within the reactor mix with the incoming unburned mixture. However, the methyl recombination reaction is still predominant, and subtracts CH3 radicals to the oxidation pathways thus suppressing the conversion of methane. It accumulates in the reactor restarting the new cycle. Summarizing, the RR analyses suggest that the oscillations occur when the high temperature branching pathway overcomes the hydrogen peroxide one, thus producing a large amount of radical species that can promote methane oxidation reactions. On the other hand, the recombination channel is very active and subtracts CH3 radicals to the oxidative channel, thus decreasing the system reactivity.

The dynamic behavior dies out when the reactions between methyl radicals and OH (CH3 + OH + M = CH3OH + M, CH3 + OH + M = CH2(S) + H2O) are fast enough to relieve the system by the inhibiting effect of the recombination reaction.

Therefore, it is clear that diluted combustion conditions force the reactive system to work in a very narrow temperature window where the competition between oxidation and recombination is stressed. This competition results in the peculiar behavior discussed above which is strongly dependent on temperature, equivalence ratio and bath gas composition.

As a matter of fact, ad hoc experimental measurements aiming at testing the effect of equivalence ratio, inlet temperature, and bath gas provide direct evidence for the kinetic foundations discussed above.

In particular, the displacement of the oscillations to a higher inlet temperature at a constant equivalence ratio, when the diluent is changed from to He or N2 to H2O or CO2, is consistent with an increase of the overall third-body efficiencies. In fact, the role of termolecular reactions is reinforced by the high collisional efficiencies of such species postponing the shift between the H2 sub-system branching reactions (H2O2 + M = OH + OH + M vs. H + O2 = OH + O), thus the onset of oscillations, to higher inlet temperatures, and reinforcing the role of methyl recombination to ethane, thus forcing the end of the oscillatory behavior at higher inlet temperatures. In the case of non-isothermal system, it has also to be considered that such species have higher heat capacities with respect to N2 or He, thus in the same operating conditions, also the reactor working temperature is lower. This can affect the relative weight of competitive reactions.

Despite the discussion here presented mostly refers to chemical kinetics aspects, none of the oscillatory phenomena can be solely attributed to isothermal kinetics. Heat release and dissipation are integral to the occurrence of such events, and are of particular relevance in real devices [START_REF] Afarin | A. 7 th Mediterranean Combustion Symposium[END_REF] . Moreover, the competition between exothermic oxidation and endothermic pyrolysis must be also taken into account.

The kinetic model here used to interpret the experimental measurements shows good agreement with boundaries on the Φ -Tin maps and for the oscillatory frequencies. However, some disagreements were observed. It has to be pointed out that the evaluation of the heat transfer coefficient used in the adopted model is affected by uncertainties. However, it was observed that such an uncertainty does not significantly change the predictions in the explored experimental conditions. Indeed, a similar satisfactory agreement is observed for isothermal conditions. However, under isothermal conditions the predicted boundaries diverge towards slightly higher inlet temperatures compared to the measurements. This further confirms the predominant nature of chemical effects in the occurrence of such phenomena.

An indirect proof of this statement is given by the fact that simulations using different kinetic mechanisms available in the literature (see Supporting Information) lead to noteworthy disagreement among model predictions themselves. This suggests that some additional effort should be dedicated to the development of detailed kinetic schemes to correctly describe fuels oxidation processes under diluted combustion conditions, at intermediate temperature.

The general importance of recombination is particularly significant and represents a major limitation in comparing experimental and modeling results since the rate constants for such reactions can vary with the type of bath gas used. Therefore, in order to match modeling to experiment with any accuracy it seems that the third body collisional efficiencies have to be carefully evaluated. Such a consideration is often either very difficult or is neglected.

For all the above reasons, there is still much scope in developing and refining models for the simulation of the large variety of observable dynamic behaviors of methane combustion.

In this regard, the present work provides a safe guidance to reliabilities of modeling, opportunities for further research, and rational employment of efforts.

Conclusion

Two Jet Stirred Flow Reactors were used to characterize the oscillatory behavior of CH4/O2 mixtures as a function of system external parameters, such as inlet temperature, equivalence ratio, residence time, methane mole fraction, and the chemical structure of the bath gas (He, N2, CO2 or H2O), at atmospheric pressure. The agreement between results obtained in the two employed experimental reactors is certainly highly valuable and inspires considerable confidence of the reliability of results.

The experimental observations revealed several novel features of methane oxidation at intermediate temperatures in presence of high concentration of different bath gases:

(i) Oscillatory regimes were established for almost all the CH4/O2 considered mixtures in both the reactors, when inlet temperatures are approximately between 1050 and 1250 K.

(ii)

The critical inlet temperature for the onset of temperature oscillations diminishes as the mixture equivalence ratio is gradually reduced. As the methane mole fraction is increased, the oscillation region on Φ -Tin maps widens under stoichiometric conditions, but shrinks under fuel lean conditions.

(iii) When the equivalence ratio is increased above the stoichiometric value, temperature oscillations are no longer observed for He, N2 and N2-H2O as bath gases, whereas for CO2 diluted mixtures, they also appear for fuel rich conditions.

(iv) H2O and CO2 lead to a general shift of the oscillation regions toward higher inlet temperatures with respect to N2-and He-diluted mixtures.

(v) The waveforms, frequency and amplitude of oscillations are strongly dependent on inlet experimental conditions.

Summarizing, temperature oscillations occur when the reactor working temperature accesses a determined range. They strongly depend on system intensive and extensive inlet parameters as well as reactor adiabatic features. The nature of the dynamic behavior is essentially related to the chemistry promoted under intermediate working temperatures.

Numerical analyses identified such kinetic features. They originate from the competition between methyl oxidation and recombination kinetic pathways and between low and high temperature branching reactions of the H2/O2 sub-system. The end of oscillations occurs at high temperatures when the inhibiting role of methyl recombination pathways is weakened by the promotion of different methane oxidation routes. It was found that bath gases can play a crucial role to the establishment of oscillations because of different heat capacities, thus different reactor working temperatures, and because of high collisional efficiencies in termolecular reactions, thus altering the relative importance of the identified competitive pathways responsible for temperature oscillations. A more fundamental definition of collisional efficiencies and their dependence on temperature and pressure would also be beneficial to improvements of existing kinetic models.

Tab. 1

 1 Experimental conditions studied in the JSFRs for CH4 mixtures. 790 K < Tin < 1225 K, P = 1.1 atm when τ = 0.5 s (Napoli) and P = 1.05 atm when τ = 2 s (Nancy).

about 5 K

 5 . An oscillatory region (III) establishes in the Tin range 1100 -1170 K with periodic temperature variations. For this oxidation regime, two T values are reported for the same Tin in Fig. 1 corresponding to the maximum and minimum values of temperature T detected within the oscillation cycle. In this case, the measured concentration values are averaged over time because of the time elapsed for the gas sampling before the chemical analyses, thus only one concentration value is reported in Fig. 1. To better identify these conditions, the closed symbols are substituted with open ones.

Fig. 1

 1 Fig. 1 Experimental reactor temperature (T) and reactants concentrations versus Tin. Φ = 1, τ = 0.5 s, P = 1.1 atm, N2 bath gas, XCH4 = 0.033. Closed symbols represent steady states, open symbols oscillatory regimes. See text for description of I, II, III, IV zones.

  1). In particular, each point on a map corresponds to an experimental condition. Given that the present work is mainly concerned with oscillatory states, the demarcation of the different reaction modes is only provided on the maps for oscillatory and non-oscillatory regimes. Therefore, regimes (I), (II), and (IV) are all identified as "no oscillation" regimes with full circles, whereas dynamic conditions are represented by open circles. Maps are also provided with numerical results, where the temperature ranges of oscillating conditions are indicated by grey rectangles.

Fig. 2

 2 Fig. 2 Experimental (symbols) and simulated (rectangles) Tin -Φ maps. τ = 0.5 s, P = 1.1 atm, Xbath gas = 0.9. (a) N2 bath gas, (b) N2-H2O bath gas, (c) CO2 bath gas. Closed circles: "no oscillation"; open circles: "oscillation"; grey rectangles: numerical oscillatory conditions simulated with the "POLIMI" model.

Fig. 3

 3 Fig. 3 shows the Tin -Φ maps for the systems diluted in He (a, c, e) and CO2 (b, d, f). Tests were repeated for three specific methane mole fractions, i.e. XCH4 = 0.01 (a, b), 0.025 (c, d), 0.05 (e, f), at τ = 2 s and for three Φ values (0.75, 1, 1.25) with Tin stepwise of 25 K (Nancy tests). Fig. 3a identifies oscillatory regimes in the Tin -Φ plane for XCH4 = 0.01 and He bath gas. The "oscillation" region covers a quite narrow temperature range of 50 K (1075 K < Tin < 1125 K) for the fuel lean condition (L4). Oscillations for stoichiometric mixtures (S4)

Fig. 3

 3 Fig. 3 Experimental (symbols) and simulated (rectangles) Tin -Φ maps. τ = 2 s, P = 1.05 atm. (a) He bath gas and XCH4 = 0.01, (b) CO2 bath gas and XCH4 = 0.01, (c) He bath gas and XCH4 = 0.025, (d) CO2 bath gas and XCH4 = 0.025, (e) He bath gas and XCH4 = 0.05, (f) CO2 bath gas and XCH4 = 0.05. Closed circles: "no oscillation"; open circles: "oscillation"; grey rectangles: numerical oscillatory conditions simulated with the "POLIMI" model.

Fig. 4

 4 Fig. 4 Representative experimental temporal temperature profiles obtained at Φ = 1 in N2 bath gas with XCH4 = 0.033 and τ = 0.5 s at (a) Tin = 1100 K and (b) Tin = 1165 K. The initial time is arbitrary.

Fig. 5

 5 Fig. 5 refers to CH4/O2 mixtures diluted in CO2, obtained at τ = 2 s. In particular, Fig. 5a is for a stoichiometric mixture preheated to Tin = 1100 K (S7), while Fig. 5b refers to a fuel rich mixture (Φ = 1.25) preheated to Tin = 1125 K (R4).

Fig. 5

 5 Fig. 5 Representative experimental temporal temperature profiles obtained in CO2 bath gas with XCH4 = 0.025 and τ = 2 s at (a) Φ = 1 and Tin = 1100 K and (b) Φ = 1.25 and Tin = 1125 K. The initial time is arbitrary.

Fig. 6

 6 Fig. 6 shows the experimental (symbols) and numerical (lines) oscillations frequencies as a function of Tin, in two representative sets of experiments. In particular, Fig. 6a refers to stoichiometric mixtures diluted in N2, N2-H2O and CO2, with XCH4 = 0.033 and τ = 0.5 s (S1, S2 and S3), whereas Fig. 6b is related to stoichiometric mixtures diluted in He and CO2, with XCH4 = 0.05 and τ = 2 s (S8 and S9).

Fig. 6

 6 Fig. 6 Experimental (symbols) and numerical (lines) oscillation frequency as a function of Tin at Φ = 1 in different bath gases with (a) XCH4 = 0.033 and τ = 0.5 s and (b) with XCH4 = 0.05 and τ = 2 s.

Fig. 7

 7 Fig. 7 Computed reactor temperature (T) versus inlet temperature (Tin) under a) nonadiabatic conditions for different values of the heat transfer coefficient; b) adiabatic conditions for different dilution; c) isothermal conditions. Φ = 1, τ = 0.5 s, P = 1.1 atm, N2 bath gas. Open circles: ignition condition. Open squares: end of oscillations. Shadowed area: maximum and minimum values of temperature oscillations.

  role in favoring branching or termination. The results discussed in the present study indicate a rather different type of phenomenon occurring in diluted mixtures, at intermediate inlet temperatures (1050 -1200 K).

Fig. 8

 8 Fig. 8 Rate of reaction analysis for key reactions involving (a) methyl radical and (b) H2/O2 subset as a function of Tin under isothermal conditions. Φ = 1, τ = 0.5 s, P = 1.1 atm, N2 bath gas. The grey area represents the Tin range in which oscillations are predicted.

Fig. 9

 9 Fig. 9 Rate of reaction analysis for key reactions involving methyl radical (left axis) and methane mole fraction (right axis) as a function of time for one oscillation cycle underisothermal conditions. Φ = 1, τ = 0.5 s, P = 1.1 atm, N2 bath gas, Tin = 1200 K.

Acknowledgments: This work has received funding from the European Union H2020 (H2020-SPIRE-04-2016) under grant agreement n°723706 and from the COST Action CM1404 "Chemistry of smart energy carriers and technologies".

Supporting Information.

 Photos of the JSFRs used and of the thermocouple located inside the reactor.

 Experimental data.

 Comparisons of model predictions with experimental measurements using different kinetic schemes.


Roy Pea

Michael Eisenbergt

Franklyn Turbak

A Computational System to Enhance and Illuminate the Development of Scientific Thinking

Keywords: microworlds, dynamical systems, scientific reasoning

Creatures of Habit is a computer-based microworld designed to engage middlttehigh school students in the process of scientific inquiry. The system depicts a universe of interacting programmable "creaturesn whose individual behavior is guided by simple rules that may model naive psychology, physical laws, chemical affinities, and other domains. Students can create or revise creature rules and explore the resulting (and often surprising) emergent behaviors within 'artificial ecosystemsn; or they may employ predesigned ecosyetems in undertaking more structured problem-solving activities. Our objective ia for students to use these ecosystem simulations as an enjoyable introduction to a variety of scientific domains, especially the area of dpomicul rystems, a field of science where experiments with such simulations often leads theory. The system encourages a wide range of reasoning and learning central to scientific methodologypattern obervation, hypothesie formation, experimentation, data collection and analysis, and deduction. We describe the rationale behind the system, discus8 wme sample activities, and outline the system's potential both aa a learning environment and as a research laboratory for empiricial studies of scientific thought. Finally, we briefly describe the present state of our prototype Creatures of Habit system.

Introduction

uDoing science" involves learning to use complex techniques and skillsmaking observations, noticing interesting p a t t e r s , forming hypothesee and the-*Thii mearch and development project was supported by a v t frwn the Spencer Foundation awarded to Roy Pea, under a pmject entitled 'Intelligent TOOL for Education."

'Bell Labs Ph.D. Scholar ories, making conjectures, and designing and running experiments. If introducing students to scientific content were sufficient, education would be hard enough; studies show children bring ill-formed models of scientific domains such as mechanics, electricity, chemistry, and biology to school (Driver, 1985; Osborne, 1985).

Students need to learn not only a body of structured beliefs, but how t o participate in the processes of science, which are intimately related to content. Students must understand that scientific ideas are motivated and supported by theory, experiment, and argumentationnot authority. Beyond the content and proceol, of Mience, students should also enjoy doing science. They often view science as a mysterious, unapproachable culture in which they observe rather than participate. A truly effective science education should dispel this image by giving students opportunities for designing and refining inquiries. We describe a computational system under develop ment -"Creatures of Habitn (henceforth, Creatures) designed to address these science education issues by providing rich, exploratory, and enjoyable scientific experiences for the middle-to-high school years.

Creatures ie a microworld of programmable interacting "creaturesn whose behaviors are based on rules. In exploring this microworld, students can be introduced to important content in various Mientific domains, notably in the area of dynomicol rystems; moreover, the program allows for a wide range of activities central to doing eciencesuch M conjecture, theory formation, experimentation, deduction, and communication of results. Perhaps most important: Creatures is intended to lead students toward original and creative workt o have them participate in science as fledgling researchers driven by curiosity, rather than onlookers motivated by wignment.

Besides its utility for science education, we see roles for Creatures in empirical studies of scientific thought. Because it f i o r d s activities engaging reasoning skills P E A , ELSENBERG, & TURBAK such u analysis, qnth&s, and evaluation (see below), it may be used to examine how these different rkills develop. The system can thus provide a uniform environment in which to dedgn specific instructional and experimental t u b .

In the remainder of thin paper, we fimt motivate the design of Creatures by discussing the subject of dynunicd systems; we then describe the system1, outline mome leuning activities that might be undertaken with it, sketch how it reflects key leuning go& for science education, and propae a design for an experiment in which Creatures may be d to illuminate the development of scientific thought.

D ynamical Systems:

Analytic Science, Synthetic Science, and "New Wave" Science

Since Creatures is designed to provide introduction t o concepts and methods involved in the study of complex dynamical systems (Thompson k [START_REF] Thornpeon | Nonlineu Dynamics and Chaos[END_REF], we briefly describe the growing importance of such u ience. Certainly there has been an exploeion of interest in this subject within the scientific community [START_REF] Gleick | C h m[END_REF]; and this phenomenon has ramifications both for sdentific methodology and for science education.

Historically, scientific thought h a been c h u u t e rized u either 'analyticn or 'synthetic" [START_REF] Oldroyd | The Arcb of Knowledge[END_REF]. In this classic formulation, the analytic method b observing phenomena and reeking laws to ucount for them; the synthetic method involves confirming the validity of laws via prediction and experimentation. The advent of scientific computing has added a new texture to thin division. In studying a rystem, a u ientist may ucount for ita behavior by conrtructing an abstract model that can be realized M a progrun; that model may then be simulated by running the progrun as a test of its applicability. Thus, there b an anolytic ride to computer modeling (observing red-world phenomena and designing model s y s t e m that might fluminate them); and a rynthetic side (changing parameters in the model system to match oboervations or we new phenomena). Hut and S w m m (1987) d t &be this approach aa 'analysis by synth&"; while Farmer and Packard (1986) call it "new wave science" a methodology "characterized by attempts at synthesis rather than reduction, cutting u r o m conventional disciplinuy boundruies .... New 'effecb' u e d i s covered through a combination of insight and serendip ity, m d more often than not experiment leads theory .... Simulations u e frequently Iwd to develop qualitative insight, often by rtudying highly simplified models which u e nonethela complicated enough to poasew universal properties found in more complicated systems."

Developing computer simulations ae experimental systems h u led to a blossoming scientific literature 6eyond cognitive science. ~imuiations are employed to discover robust properties of intergalactic collisions, 'Crutura exuu in prototype fonn on a Hewlett PIcksrd Seria 300 Model 320.

heat flow in rolids, kinetics of chemical reaction mechd a m s , evolutionary adaptation, and many other dynunical systems. These simulations enable the study of incre&ngly complex phenomena; but they .Lo t& skill in developing formal model, and p r e n t studentr with new phenomena to underatand and methods to use. Techniques of developing, rtudying, and documenting computer models have become an integral element of rdentific method. In thin context we see Creatures as having p u t i c u l u potential d u e , since it provid-an engaging ''introduction to complexity" and an environment in which to learn-bydoing these skills of "andysb by rynthesis." As we nill show below, studenta can construct systems of creatures that exemplify concepts such rs stability, oscillations, and bifurcations. Because these systems may be of their own devising, the concepts become personalized in a way that canned demonstrations could not.

Creatures of Habit: The Basic Elements

Creatures b an environment in which students can explore "artificial ecwystemsw composed of interacting programmed entities c d e d 'creatures." Inspired by [START_REF] Braitenberg | Vehides. Cunbridge: MIT Press & S a r a , A[END_REF] Vehicles, these creatures b e have and interact on the computer screen according to ub of simple rules. The rules might be chogcn to reflect physical laws, naive pychology, chemical &nities, and K) on. A simple example of a group of creat u r a and how they might interact illustrates the system'r componenta.

Consider the scenario illustrated by the sequence of " m a p h o b " in Figures 1 to 5. ' Here we see three different creature types (distinguished by geometric rhape) with the following behavioral rules:

1. Squanr u e attracted to oquues, but indifferent to triangles and circles;

2. Circlu u e repelled by triangles and circles, but attracted to squues;

3. Trianglu u e attracted to circles, but indifferent to squues and triangles.

Figure 1 showa an initial "ecosystem" with five creatures, and Figurea 2 -5 demonstrate how thin ecosystem evolves over time according to these rules (note that creatures leave a visible path). Even in this simple c w we m many types of behavior: mutual attrution (between squares), mutual repulsion (between circles), and predator-prey relationships (the triangle 'chws" the circle that is 'running awayn from it). Note how the interutions can be pused into 'episodes" in a naive model of animal behavior:

The oqoues like each other A d begin to move together. The a r d e s , interested in the squares but oblivious to each other (and the triangle), tag along. The triangle, spying a tasty circle, s t a m sneaking up on it; 2Th= figura were generated by our prototype Crcat u r t r ryrtem.

P E A , EISENBERG, k TmLBAK

Afkr notiang one another, the frightened citda dec in oppodte directions; The triangle chuer the topmoot d t d e while the other d r d e , regaining ita composure, he& back towud the r q u u a . The maturer r e d a find 'equilibrium" configuration when the triangle 'catchan the topmoat a r d e and the other citde joins the two q a u e s .

Thia dmple v c n u i o indicates how even a small number of creatara governed by rtraightforrud rules can l e d to interesting, vuied behavior. Under slightly different situations, it t poodble to obsene more complex behavior and 'emergent phenomena"3. Even this scenuio ruggestr intelbting questionr to explore, e.g., how do initid pbsitionr of the creatures influence the epiaodic nature of the interution or the creature "fates"? How would adding one creature to the start state d e c t the interaction? Given these species, is it pousible to design a configuration that exhibits rtable d a t i o n r ?

We now preacnt a fuller exposition of the basic c l e menta of the Creatures of Habit ryrtem:

Creature Morphology

Creatures u e rmdl mobile k n t ~y creatures that live" on the m e n .

Creatures come in many 'species,' identified by a set of discrete chuacteristia. Above, creatures were dirtingubhed only by rhape, but more mimal-like m a t u r e r are conceivable (a species might be identified u purple, with bobbling eyes, pointy earn, rquiggly t d ) . Species may d m be endowed with l e u visible propertia such u m u or birth rate.

Specisr Rules

Creatures interact with d other m d their environment by obeying a r m d set of rpeaesrpedfic rules I linking perception to action. Typically, species rules indicate which properties u e "attractive" or "repellent" to a epecies. Once a rpedea ruleset t defined, every member of that rpecia u r a those rules to govern ita behavior.

Ecoeyetem Rules

The complete dynamics of a population of creatures b determined by a higher-level u t of ecoryrtem rules which specify how the creatures employ the notions of 'attraction" and 'repulsionn when moving. Ecorye tem rules u e p e r h a p b a t illustrated by exlmple; the p u t i c u l u exampla that ue provide rhould be r e d merely u indicating a rprce of p k b i l i t i e s . Vuiations on these r u l e u t a u e of course pordble, and other r u l e u t a u e computationally viable. specifying that collirionr between creatures u e elmtic and momentum-connerving4. Note that this rule-net implicitly r u u m e r that m u s will be included among the rpecia properties. Creatures here u e regarded This mleset embodies a simple 'naive psychologicalw model for creature interactions; it w u the ecooystem rule-set implicitly used in the scenario depicted in the figures. To determine how much creatnre CJ would "likelnot-likew creature C2 if it were one unit away, look a t Cl's species rules and subtract the number of repellent properties exhibited by C2 from the number of attracting properties. To determine how much C1 actually likes C2, divide the previous result by the di* tance; this inverse distance dependence models a dtuation in which creatnres pay more attention to d o ~r objects. Thus, if creature C l finds three properties of creature CZ attractive, and one repellent, then the unit-distance strength of Cl's attraction for C2 is 2; if C1 is 10 units of distance from C2 at that moment, then Cl's current attraction toward C2 hrs a strength of 0.2. To determine how creature CJ rhould move, the rystem finds the creature about which C1 f& most strongly (like/notlike); creature C1 rhould then move at a constant (species-specific) speed towards that creature (if liked) or away (if not). Enwunter-Dependent Birth and Death R a t u . In thin r u l t s e t b u e d on population biology, we usume that creatures have speciebspecific properties of 'birth and death rates." Creature movement might be governed by the ruleact above, and in addition new rpecies members may be born or die. The birth rules might depend upon encounter rates between plemberr of the u m e species. The birth rule for rpeaes A might be the following: on any time step, an A-creature hrs a 2% chance of giving birth to another A-creature, except when the previous time step resulted in contact with m A-creature, in which c u e there is a 50% chance of giving birth. The death rate of a species might be defined amilarly (e.g., on any p u t i d u time r k p , an

I d d Gor

A n e a t u r e h u a 3% chance of dying, except when the previous step m a l t e d in cont.ct with a 'predatorw rpecier B, with a 60% chance of dying).

Interface to the System; Additional 'hob

We have dacribed the brric elementa of Creatures but thur far ignored the interface. In this s t i o n , we describe several upectr of the interface design for Creatures. Many h u a murt be resolved to make Creatures both a c c d b l e to fimt-time werr and extensible for w by expert ocientista. We rim to fimt provide -me beginning ecoryrtem de-sets for e u y ucess to Creatures; given a choice of ecwystems, the user may edit the properties and rules for individual species with iconic menus. The basic interface should alm provide a fair amount of control over simulations: initial configurations may be specified by directly moving creatures to their desired starting paitions, and individual runs may be paused, single-stepped, replayed, and stored.

The description above does not address the needs of more experienced users who wbh to create or edit the ecooystem mle-acts. The range of ecosystem rule-sets may be constrained no that they might be customized via a menu interface; but more likely, a more elaborate interface to ecooystem rulcaets (verging on a specialpurpose programming language) will be required.

We u e working on other rpecial-purpoee tools for the Creatures system for meů modifying, and experimenting with configurations and rules. For ex-A p l e , we aim to provid;multiple screen windows ao the rtudent can make two runs at once, comparing results from two systems by varying only one parameter say the initid poeition of a creature. Experimentation will r b o be supported by 'unotation windowsw where students enter information obtained from the results of previow rum; e.g. to examine a sequence of runs and make deductions about the species-rules of individnal creatnres, structured tables help a student codify the information obtained and highlight the information needed for a mlution. The system will indude graphical memuring devices that enable the rtudent to take precise me~urementr of distance, angle, and time; in our experience with the Creatures p r o t e type, these variables have proven crucial in describing important features of simulations.

Sample A c t i v i t i e s Exploratory Activities

The moot fundamentd activity we imagine is e x p b ration. Here, students u e provided with seta of creatures in an ecorystem, and o b u n e how the creatures behave and interact. Since both the species rulc-nets and the ecwystem mle-acts u e uceaaible, rtudenb can investigate the relationship between the behavior they observe and the rulesets governing that behavior. And because the system facilitates change a t a wide variety of levels (from the number u d positioning of the creatnres t o the species and ecosystem mle-bets), students can readily explore the effects of 'perturbing" a given s c e n d o in diverse ways. Even inexperienced students could examine how the qualitative behavior of a configuration depended on the initial position of one of the creatures; more advanced rtudents could explore the effecta of "trrr.kingn rpeciu rule-seta; yet more advanced dndentr could work at the ecorystem rule-act level.

Even within thia purely exploratory frame, a number of importmt acientific idem u e introduced. To highlight, it quickly becoma apparent in working with Creatures that complcz rptemr moy orue even from vcty rimple ruler. Often these r y s t e m exhibit emergent propertiesrtability, d a t i o n r , inevenibiditythat r a b t explanation in terms of low-level rules. Moreover, by playing uound with these ryat e m , students can develop important rkills that u e naually alighted by classroom science: making conjectures, warching for patkrxu in data, and generating qualitative or statistical descriptions. We emphasize that although exploration in a low-threshold activity, it ia not merely m introductory one; rather, it b the essential point of the design of Creatures. Indeed, although Creatures is targeted for students in middle-t+ high school, issues r a i d by 'mere exploration" could evolve into highly challenging problems at or beyond the college level.

Problem-Solving Activities

Although Creatura lends itself t o exploration, it supports more structured Uproblem-mlving" activities m well. We have used the prototype to play a "My* tery Creatures" game where usen are provided with a set of creatnres whose rpeaes rules are unknown and must be determined from their behaviors. The goal of the game is to design experimental runs m d thereby deduce the rultsets of the creatures. Sometimes, one well-chosen experiment can illuminate the behavior of a number of mystery creatures all at once; in other situations, the user must design a sequence of experiments. A related game, Invisible Creatures, involves several visible creatures (whoae mle-sets u e ucessible) that interact with UI "invisible" creature. in this caw, the goal is to determine the identity and location of the invisible creature b a d on the observed actions of the visible ones. These kinds of activities stress skills of experimentation rnd logical deduction.

D e s i g n Activities

Many u t i v i t i a supported by Creatures have an important design focus. Here the emphasis is not on why a given rystem evolves m it does; rather, the goal k to wnrtruct a ryrtem exemplifying some d t aired behavior. For example, given a aet of creatures, one might mk: ' L there an initial configuration for t h e creatures much that all the creatures will collide a t once?" Design can take place at many levels, from constructing configurations of creatures to ucreating new worlds" with new ecosystem rules.

A Hypothetical S c e n a r i o

The following scenario shows how r rtudent might use Creatures. We use the same rpeaeil and ecosystem rules M in the ~e n a r i o shown in the figures earlier: b The student begins with a screen on which a dozen creaturalome mqnuer, triangles, and circlermove in nrionr directionr, creating intricate patterna. The rtndent tries similar runs uring the r u n e nnmbem of creatnrea, but with different initial configurations, just to get a feeling for the kin& of phenomena that might occur. b The rtudent decider to figure out what the rpecies rules for e u h of the three creature types might be. She rum a Kqnence of experiments in which two a e l tuxes of the u m e type u e placed a t a alight distance from each other in the m e e n center. She oberves that when triangla u e rued, they rtay where they began; q u u u move towud each other; and circles move away from e u h other. She concludes that triangle creatures u e indifferent to other triangles, while quares u e attracted to other squares and circles are mutndiy repelled. b When she places a square and two circles in a particular initial confignration, the two circles move toward the square until, a t a later time, they drift apart; then, shortly after, they move towud the square .gain. The student hypothesizes that arcles find the q u a r e attractive, but when they are too dose to each other, their mutual repukon becomes stronger than their attraction for the square. She replays the earlier run, thia time occuionally pausing the dmulation and mearuring dirtances and headings of the creatures to test her conjecture. b Replaying the previous run, the student notices the square remained stationary throughout. She concludes that squues u e indifferent to circles. In fact, she r e c d e that in her eulier imulations, whenever there was only one square on the acreen, that square did not move a t dl regudleaa of how many circles and trirngles were about; so she concludes that equares must be indifierent to triangles, too. b During her experimenb with circles, the student remembem noticing that the two creatures ended up standing rtill in o p p i t e halves of the screen.' She decide8 to investigate: Does m y initial confignration of m y number of circles always end up in a "stationary" state; and if so, do the find states reveal some pattern?

Clearly, the scenario above could be extended in different ways. For example, the species rules c o d d involve features besides shape; or the creatures c o d d move at vuying speeds; or the rtndent might choose to exunine whether any ustationary configurationw r e mains stationary if one other creature in added to the acreen; or the rtndent could design a new creatnre that moves away from squares, and see what happerm when it b placed in a crowd of (mutually attracting) squues.

From this very rimple beginning, many projectanome touching on very sopbticated questionsmay arbs.

Creatures as Part of A S c i e n c e Curriculum

A perennial tension in science education elcisb between teuhing content knowledge m d scientific It ahould perhapa be mentioned that, for this particulsr scenario, we -ume that the .cmn on whir31 the c n a t u r u move 'naps around" in both the x-and y-directions.

PEA, EEENBERG, & TURBAK method. The flaw in oreremphaaizing content t that it k n d r to l e d to m k leuning, with little senre of how uientistr work to develop theories. T h u approach t y p i c d y friL to recognize rtudentr' preexisting conceptionr of uientific domainn (diSusa, 1987); rtudents u e t u i t l y encouraged to memorize laws m d rerults without ever engaging their own reuoning powerr to m how they were urived at, or could be tested. In contrast u 'dkovery learning", which e m p h u i z a experimentation m d the scientific process. The difficulty with thL approach u twofold. F h t , the notion of "di.coveryn is often realized in prutice by cmned experiments; the 'diecovery," whome resulb are known beforehand, is ruperfiaal. In this kind of setting, rtud e n b spend much of their time trying to find the 'right mrwcr," rather than what the experiment means. Sec- ond, when discovery learning is realized mwt sincerely, by rimply letting the student explore on his or her own, the teacher is often left with question8 about what the student's explorations were, m d what was learned (Hawkins, 1987).

We feel that Creatures can serve as a medium for guided, exploratory discovery-type learning (d White m d Horwitz, 1987). The worlds students can explore with Creatures are brand-new; there are few "right answera," and real discoveries can be m d e . But while making these discoveries, rtudentr c m obtain a rtructared introduction to the role of conjecture, measure ment, and experimentation in the Mientific process. Activities such as the 'Mystery Creatures" game place m emphasis on logical reasoning; these puzzlolike activities can provide the benefits of occasional settings in which a right answer doen in fact exiet.

An Experimental Design for Studying the Development of Scientific Thought

We have stressed the role of Creatures in sdence education. Here we illustrate how it may .Iso be used for rtudying the development of raentific thought, propwing an experimental design we aim to undertake using the next iteration of the system. Consider the following three t a s l , corresponding to three common modes of scientific thought, that a rtudent might be asked to perform using a particular ecosystem and starting creature configuration:

1. The rtudent watches the system run (i.e. creatures move on the screen), m d ir anked to infer the r p e d a r u l a of each creatnre. This onolytic trek involva working backward from observed behavior to underlying rules.

2. The rtudent ir rhown the species rules for the aeaturer, then asked to predict the result of running the ryrtem uring the rtartup configuration. Thir ia a rynthetic task of predicting b e havior from a known set of deterministic rules.

3. The rtudent ir provided the set of species rules for the creatures, and shown the results of running the ryrtem using the startup configuration; his task would dmply be t o explain events ob-~r v e d in the ecenario in terms of the rpecics rules. Thir is m eualuotiu, or p e r h a p esplonotory, t u k . Now the very u m e ecosystem m d rtuting configuration could be uned for e u h of t h e three krk, u given to three different groupr of rtudentr; but the kindr of dercriptioru generated for each of the three t u k rm d indeed, the featurw of a p u t i c u l u uc nario attended tomight r u y widely. By obtaining think-aloud protocoh from that different subject group, n wold begm to develop a coherent picture of how the several m o d a of rdentific thinkinn differ: rpedcally, we could begin to distinguish t h e k mod& according to the sort8 of phenomena to which they typically apply. h a p o ~i b i i t y , we might u k whether the phenomena uoed to deduce rules in the 'analytic" tasks are the Mme phenomena for which explanations u e offered in the 'evaluative" task. In other words, are the phenomena that people find rdient for the deduction process the u m e u those that people find representative in the explanatory process? For example, in looking for underlying ruler, rubjecta might tend to focus on events in which creatures reverse direction (e.g., a situation in which creatnre A, in moving toward an "attractive" creatnre B, fin& itself too n e u a Yrepellentn creature C and turns uound). In contrast, rubjecta given the task of explaining a given ecenario might indrt on a chronologically faithful nurative; i.e., they might pay as much attention to explaining creature A'r initial direction rs they do to explaining i b change of direction.

Or we might try to churcterize the kinds of 'epirodic groupingr" that different task groups a d g n to the given scenuiw: e.g., u e there certain t a ~h in which ecenariw u e typically viewed u "punctuated" by changes in direction? Are there certain tub in which the find rtates (e.g., the achievement of equilibrium) u e attributed more importance? Under what arcumrtances do people reuon a t the "ryskrn level," talking about behaviors of luger g r o u p of creatures, rather t h m a t the 'atomic level" of individud creatures' hirtories?

Yet another h u e involves the possibility of m order effect between different activity modes. For example, in a diAerent rtudy we could rek rubjecb who have just predicted a scenario in the "synthetic" t u k t o watch the u t u d scenario and perform the "explmatory" t u k ; their protocob could be compued with thwe generated by rubjecb given the explanatory task alone. We could then begin to characterize how a r u b ject'r previous predictions affect the kindr of explanatiom generated for a p u t i c u l u ecenario.

Finally, it is worth mentioning the issue of noticing "creature indifferencen in the context of thin experimental design. In our attempts a t watching rcenuion and deducing the rules underlying them, we found that a great deal of information ir conveyed by noting which creatures u e indifferent to which others. Indifference between creatures seems to be a rubtler notion that attraction or repulsion; when watching a scenario, there is a tendency to attribute a creature's movements solely t o attraction or repulsion without

PEA, EISENBERG, & TURBAK

recogniring that the obrena of motion towud or away from wme other crutnre im rLo informative. Thir kndency is reminiwcent of the disinclination to use negative information in dentific reuoning (Mynatt ct. d., 1977;Wuon, 1977). Thlu another quation dd r a u b l e within thu experimental daign im whether any of the three t u b im more likely to l e d mbjecta to confront the h u e of creature indifference.

Conclusion:

Prerent Status, Future Directionr We now have two prototype Creatures of Habit ryetern: one written in Scheme (a Lkp dialect), one in C. Both run on a Hewlett-PLclrud Series 300 Model 320 computer. The Scheme version containn facilities for developing and editing uwystem and species rule ~b , and induden featurea for maintaining histories of individual rune. But it t slow and a t the moment unsuitable for running ryrtems with more than 4 5 creatures. The C program b less elaborate or interactive, but it. running speed b two ordern of magnitude b t e r than the Scheme program; using the C prototype, simd a t i o ~ employing 16 creatures have been run at a ntkfactory rpeed. We continue work with both rye k m s , using Scheme .r, a medium for developing new tools and trying out ecwystems in the small, and the C progrun to run rimulationr quickly .nd experiment with properties of larger creaturoworldn.

Much of our current daign effort t developing a suitable interface to the system, and delineating a ''.tarter act" of sample ecosystems appropriate for middlctc+high school students. A eecond theme noted above is using the system M a laboratory tool for conducting experimenb in the development of Kientific thought. It rhould also be porrible to follow up with work in student modeling and intelligent scorching" efforta to promote the varioun inquiry skills in ocience we have outlined. Ultimately, we hope that Cree tor-of Habit will provide an environment in which students c m acquire the concepb, the methodological techniques, andtoo frequently neglectedthe intellectual curiosity of the working Kientlt.

Figure 1

 1 Figure 1: Initial Configuration

  Rule-Set. k thia simple model, all creatures u e regarded u indifferent to one mother (i.e. p r o p erties that u e attrutive/repubim to a given ere* ture'r rpedes do not affect its motion). The ecwyrtem ruleset doer indude, however, a default collision rule =Far btrpcr, ur entire durta of cruturu might move u one unit due to the internal pattenu of attractim m d rcpukion between the duster's individual.

  u i n h a t e b i u d -b d -l i k e objects, each with a ( r p e c i ~p e c i f i c ) m u . Inuerre-Square Force Rule-Set. This mle-net d a b rates the last one. Creatures u e partickr upon which inveraesquued forces act. To determine the individual force of creature C2 on creature C1, we rum the attractions (+) and npulsionr (-) that CJ h u for C2's properties. The resulting number divided by the square of the distance between CJ and C2 is the magnitude of this individnal force; the direction between CJ and C2 determines the direction of the individnal force. The total force on creatnre C1 is the vector sum of the individual forces from all other creatnres. Dividing this total force by creatnre C 2 s maas determines the acceleration, or velocity change, of Cl. (Again we ursume creatnres have maaa M a species-specific p r o p ertyand we might add mme other properties, ruch .B M well.) Dirtance-Dependent "Animd Attmction" Rule-Set.

'It ia rl.o possible to collide with the boundaria of the ureen; these colliaiona &odd be w e d u elmtic collieions with a rbtionary wall.