Development and Application of a ReaxFF Reactive Force Field for

ARTICLE pubs.acs.org/JPCA

Development and Application of a ReaxFF Reactive Force Field for Hydrogen Combustion Satyam Agrawalla and Adri C. T. van Duin* Department of Mechanical and Nuclear Engineering, The Pennsylvania State University, University Park, Pennsylvania 16802, United States

bS Supporting Information ABSTRACT: To investigate the reaction kinetics of hydrogen combustion at high-pressure and high-temperature conditions, we constructed a ReaxFF training set to include reaction energies and transition states relevant to hydrogen combustion and optimized the ReaxFF force field parameters against training data obtained from quantum mechanical calculations and experimental values. The optimized ReaxFF potential functions were used to run NVT MD (i.e., molecular dynamics simulation with fixed number of atoms, volume, and temperature) simulations for various H2/O2 mixtures. We observed that the hydroperoxyl (HO2) radical plays a key role in the reaction kinetics at our input conditions (T g 3000 K, P > 400 atm). The reaction mechanism observed is in good agreement with predictions of existing continuum-scale kinetic models for hydrogen combustion, and a transition of reaction mechanism is observed as we move from high pressure, low temperature to low pressure, high temperature. Since ReaxFF derives its parameters from quantum mechanical data and can simulate reaction pathways without any preconditioning, we believe that atomistic simulations through ReaxFF could be a useful tool in enhancing existing continuum-scale kinetic models for prediction of hydrogen combustion kinetics at high-pressure and high-temperature conditions, which otherwise is difficult to attain through experiments.

1. INTRODUCTION Hydrogen is currently considered as the cleanest renewable energy source. Unlike fossil fuels, there are plentiful resources of hydrogen everywhere on this planet. Its major advantages are its high heat of combustion and zero CO2 emission, which is important from a global warming perspective.1,2 However, there are several issues, like storage3,4 and safety,5 which need to be resolved before hydrogen can be used as a fuel for practical purposes. The prospect of a hydrogen-based economy has increased interest in the use of hydrogen as a fuel. In order to extensively use hydrogen as fuel, a detailed understanding of the reaction kinetics of the hydrogen combustion process is necessary at various input conditions. The reaction kinetics of H2/O2 has been studied extensively in the past.6-11 The reaction kinetics is well-understood for regions below the classical extended second explosion limit,8 which comprises low pressures and high temperatures and is dominated by chain-branching reactions shown in reactions R1 and R2. H þ O2 f OH þ O

ðR1Þ

O þ H2 f H þ OH

ðR2Þ

r 2011 American Chemical Society

However, the kinetics above the extended second explosion limit, which comprises high-pressure and low-temperature conditions, is less well-understood. Mueller et al.8 have suggested that the kinetics in this zone involves the formation and consumption reactions of HO2 and H2O2, which are not extensively studied as of now. Among the current kinetic models for H2/O2 combustion, Li et al.9 have validated their model with experimental data for a wide range of conditions (298-3000 K, 0.3-87 atm, φ = 0.25-5.0) and found excellent agreement. Strohle and Myhrvold10 have evaluated several existing kinetic models of H2/O2 combustion under gas turbine conditions (pressure up to 33 atm) and concluded that the O’Conaire mechanism11 and Li mechanism9 are most appropriate for predicting reaction kinetics under these conditions. However, at very high pressures (P > 400 atm) it is difficult to obtain experimental data and none of the existing continuum-scale kinetic models have been validated for such elevated pressures. Received: September 1, 2010 Revised: December 17, 2010 Published: January 25, 2011 960

dx.doi.org/10.1021/jp108325e | J. Phys. Chem. A 2011, 115, 960–972

The Journal of Physical Chemistry A

ARTICLE

High-pressure combustion has applications in rocket engines and shock-induced reactions. Therefore, an accurate, computationally feasible, atomistic method would be useful to describe the complete reaction kinetics and enhance the understanding of the reaction kinetics under such elevated pressures. Quantum mechanical (QM) methods are the most accurate atomistic methods known for providing barriers and reaction energies for individual reactions. However, QM-based methods are computationally expensive; although they can be applied in dynamical simulations, these simulations are limited to small systems and short time scales,12 making it difficult to use them for getting a detailed, statistically relevant description of the H2/O2 reaction chemistry. Empirical methods, including tightbinding13,14 and force field15-23 (FF)-based approaches, provide a computationally inexpensive alternative to QM methods. Although FF methods are mainly used for nonreactive systems,15-19 currently a number of reactive FF methods are available20-23 that enable nanosecond-scale simulations on systems .1000 atoms. Here, we present the development of an application of a ReaxFF reactive force field24 for hydrogen combustion reactions. ReaxFF has earlier been applied to a wide range of materials, including organic reactions, 24 energetic materials under extreme conduction,25,26 decomposition of explosives,27 thermal decomposition of polymers,28 heterogeneous catalysts,29 fuel cells,30 crack propagation in silicon crystals,31 dissociation of H2 on Pt surfaces,32 storage of H2 in Mg nanoclusters,33 catalytic formation of carbon nanotubes,34 tribology of metal-metal oxide interfaces,35 and hydrocarbon oxidation.36 ReaxFF is a quantumderived force field which retains nearly the accuracy of results from QM for both reaction energies and transition states. It was developed in order to simulate the formation and dissociation of bonds correctly and does not require predefined reactive sites or reaction pathways to simulate a chemical reaction. Since the FF parameters are derived solely from QM, ReaxFF can be directly applied to systems that cannot be studied experimentally. In this paper we will describe the development of a ReaxFF force field to carry out molecular dynamics (MD) simulations for investigating the chemical processes involved in hydrogen combustion. We present here the optimization of ReaxFF force field parameters by training against QM-derived transition states and chemical reactions that are relevant to hydrogen oxidation (section 3.1). We will then present results from NVT MD (i.e., molecular dynamics simulation with fixed number of atoms, volume, and temperature) and NVE MD (i.e., molecular dynamics simulation with fixed number of atoms, volume, and total energy) simulations in ReaxFF of gas-phase reactions of hydrogen/oxygen mixtures at high temperature (T g 3000 K), high pressure (P > 400 atm), and equivalence ratio (φ) ranging from 0.25 to 4.5 (section 3.2). A detailed analysis of the MD simulations is performed to determine initiation reactions, important intermediate reactions in the simulation, and overall energy barrier of the hydrogen combustion reaction. These simulations demonstrate the ability of ReaxFF as a computational method to describe the chemical events involved in the hydrogen combustion process and should provide an accurate tool for analyzing the reaction mechanism under a wide range of input conditions.

of atoms by solving Newton’s equations of motion; this methodology is popularly referred as molecular dynamics. ReaxFF provides accurate description of bond formation and dissociation between atoms during a chemical reaction. Unlike nonreactive force fields, ReaxFF determines connectivity between a pair of atoms through their bond order, which it calculates at every MD step based on interatomic distance. ReaxFF accounts for nonbonded interactions like van der Waals and Coloumbic interactions for a reactive system by calculating interactions between all pair of atoms, irrespective of connectivity. Excessively close range nonbonded interactions are avoided by using a shielding term. ReaxFF also accounts for polarization effects through a geometry-dependent charge calculation scheme.37,38 A list of the potential functions used by ReaxFF can be accessed through the supporting material of Chenoweth et al.36 The potential functions used in this application of ReaxFF are the same as reported previously36 with one exception: we modified the lone pair potential function to enable ReaxFF to distinguish the heat of reaction for reactions R3 and R4. ðR3Þ H2 O f H þ OH OH f O þ H

ðR4Þ

Both the reactions R3 and R4 involve the breaking of an O-H bond. However, the O-H bond in R4 is easier to break than that in R3. This is due to lone pair stabilization on the O atom in reaction R4. Earlier ReaxFF was not able to distinguish this, because it was only allowing up to two lone pairs on an oxygen atom, which is not connected to any other atom. In order to account for this effect, we added a new parameter vpar2 in the lone pair potential function of ReaxFF, which ensures that the optimal number of lone pairs on an oxygen atom without connectivity is three instead of two. The modified lone pair potential function is shown through eqs 1-5. Ulp ¼

plp2 Δlp i lp

1 þ expð - 75Δi Þ

ð1Þ

where lp

Δi ¼ nlp, opt - nlp, i

ð2Þ

nlp, opt ¼ 0:5ðValei - nval, i Þ þ vpar2

ð3Þ

e e 2 ! Δi Δ e nlp, i ¼ int þ exp - plp1 2 þ Δi - 2 int i 2 2 ð4Þ Δei ¼ - Valei þ

neighborsðiÞ X

BOij

ð5Þ

j¼1

Equation 1 represents the lone pair potential function for an atom i. plp1, plp2, and vpar2 are parameters of the ReaxFF force field. Valei is the total number of free electrons in the outer shell of atom i and nval,i is the valency of atom i. BOij is the bond order between atom i and its neighboring atoms j, calculated based on interatomic distance. A penalty energy given by eq 1 is applied to an atom i, if the number of available lone pairs on the atom given by nlp,i is different from the optimal number of lone pairs nlp,opt.

2. METHOD 2.1. ReaxFF Reactive Force Field Method. ReaxFF is a reactive force field which can simulate chemical reactions among a set of atoms. It can be applied to simulate dynamics of a system 961

dx.doi.org/10.1021/jp108325e |J. Phys. Chem. A 2011, 115, 960–972


ARTICLE

Table 1. Temperature, Density, and Equivalence Ratio of the Systems Studied ensemble

temp (K)

density (kg/dm3)

equivalence ratio (φ)

NVT

3000-4000

0.03-0.25

0.25-4.50

NVE

2000

0.03-0.25

1

The ReaxFF force field parameters have been optimized against an extensive training data set obtained from our QM results and QM data from literature, consisting of transition states and reaction energies relevant to hydrogen oxidation reactions. The parameters of the ReaxFF potential functions were optimized through a single-parameter search optimization39 performed for parameters relevant to hydrogen oxidation, to minimize the following sum of squares: " #2 n X ðxi, QM - xi, ReaxFF Þ ð6Þ error ¼ σ i¼1

Figure 1. Cubic periodic box of side length 25 Å containing 67 H2 and 33 O2 molecules (left). The red balls portray O atoms, and white ones represent H atoms. NVT MD simulations are run at 3500 K for 0.5 ns after which the content of the periodic box changes to 57 H2O, 9 H2, 4 O2, 1 OH, and 1 H (right).

the parameters. We have added reactivity and transition states of all the important reactions in H2/O2 chemistry in the training set. Our training set consists of standard heat of formation values of key reaction components of hydrogen combustion (i.e., H, O, OH, HO2, H2O2, and H2O). It also consists of the heat of reactions and energy barriers of the key reactions of hydrogen combustion. In addition, the training set contains data for bond stretching energy and valence angle distortion energy between all combinations of H and O atoms in the training set. The parameters of the ReaxFF force field were obtained by training the force field against these data, using the method described in section 2.1. The following section describes the performance of the optimized ReaxFF force field against the training data. 3.1.1. Bond Stretching Curves. Bond stretching energies for different combinations of bonds containing H and O atoms were calculated using DFT/B3LYP/6-311G**. The ground-state energy of each bond was calculated through full geometry optimization. Bond stretching curves were calculated by constraining the bond length of interest and optimizing the rest of the geometry. Figure 2 shows the bond stretching curves obtained for the H-H bond in H2, O-H bond in H2O, OdO double bond in O2, O-O single bond in H2O2, and O-O single bond in HO2, respectively, for ReaxFF and QM. As seen in Figure 2, the points at zero energy represent the ground-state energy of the corresponding bond. The portion where the curve starts to flatten out represents the dissociation limit for the bond. For the stretching curve of the O-O bond in HO2 (in Figure 2) we have calculated multiple spin states in QM. This is because the quartet state calculations yield lower energy values compared to doublet state when the O-O bond is near the dissociation limit. Since ReaxFF does not include the concept of multiple spin states and is parametrized to reproduce the energy corresponding to the lowest energy state, we have added quartet-state energy values for the large separation distance of the O-O bond and doublet energy values for the near-equilibrium bond distance. Overall, we can say that ReaxFF captures the trend of bond stretching quite well, especially near the bond equilibrium distance. At extreme bond compression ReaxFF starts to deviate from the QM energies; the Morse potential employed in the ReaxFF formulation underestimates the compression energy. These highly compressed bonds are highly uncommon, even at the elevated temperatures encountered in combustion environments. As such, we believe that this discrepancy should not

where xi,QM represents the QM-calculated value and xi,ReaxFF represents the ReaxFF-calculated value of the transition states and reaction energies fed in the training set. σ represents the accuracy specified in the training set. 2.2. QM Calculations. For obtaining training data, QM calculations were performed using a B3LYP hybrid DFT (density functional theory) functional40,41 and the Pople 6-311G** basis set42 as implemented in Jaguar (version 7.5).43 Heat of formation data was obtained from G3 energy values derived by Pople et al.44-47 2.3. MD Simulations. MD simulations were performed using ReaxFF potential functions. In our simulations we have studied H2/O2 mixtures for a wide range of input conditions. The simulations were performed with fixed number of atoms (N), fixed volume (V), and a fixed temperature (T). We designate these conditions as NVT. The temperature is controlled using a Berendsen thermostat48 (damping constant equals 0.5 ps). We have also studied mixtures of H2/O2 seeded with one OH radical at low initial temperature. Simulations of such mixtures were performed with fixed number of atoms (N), fixed volume (V), and fixed total energy (E) (referred as NVE). A summary of all the input conditions can be seen in Table 1. The simulation procedure begins with creating input structures in a box with periodic boundaries. This is followed by an energy minimization of the structures (using the Shanno conjugate gradient method49). Using a time step of 0.1 fs, MD simulations were started by assigning random velocities to each atom in the input structure. ReaxFF then uses the velocity-Verlet scheme50,51 to calculate the trajectory of atoms. Simulations are performed for a period of 0.5 ns. Figure 1 shows the input structure and output structure of an NVT MD simulation run at 3500 K with 67 H2, 33 O2 (N = 200) molecules in a cubic periodic box of side length 25 Å. The analysis of the simulations was performed by using a bond-order cutoff of 0.3 Å. The bondorder cutoff does not affect the simulations but only the interpretation in terms of chemical components.

3. RESULTS AND DISCUSSION 3.1. Optimization of the ReaxFF Force Field. In order to obtain a ReaxFF force field which predicts accurately all the aspects of H2/O2 chemistry, we required adequate data to train 962



ARTICLE

Figure 3. QM(B3LYP/6-311G**) and ReaxFF energies for distortion of the H-O-H angle (in H2O) and H-O-OH angle (in H2O2) calculated by constraining the corresponding bond angles and optimizing the rest of the geometry. Figure 2. QM(B3LYP/6-311G**) and ReaxFF bond stretching energies for the H-H bond (in H2), HO-H bond (in H2O), OdO bond (in O2), HO-OH bond (in H2O2), and HO-O bond (in HO2) calculated by constraining the corresponding bond lengths and optimizing the rest of the geometry.

Table 2. Comparison of Standard Heats of Formation at 0 K of Key Species of Hydrogen Combustion in QM and ReaxFF ΔH0f (0 K) (kcal)

greatly affect the ReaxFF accuracy. For applications requiring highly compressed bonds (e.g., shock compression) we recommend adding a more repulsive van der Waals term to the ReaxFF potential. 3.1.2. Angle Distortion Curve. Similar to bond stretching curves, angle distortion curves were calculated in QM for the H-O-H angle in H2O and H-O-O angle in H2O2. The ground-state energy of each angle was calculated through full geometry optimization. Other points on the angle distortion curve were obtained by restraining the internal coordinates of interest and optimizing the rest of the geometry. Figure 3 shows the angle distortion energy curves for the H-O-H angle in water and H-O-O angle in hydrogen peroxide. ReaxFF captures the trend of the H-O-O angle distortion as predicted by QM; however, the ground-state angle in ReaxFF shifts from 100° to 106°. We could have improved the ReaxFF prediction by increasing the required accuracy for the H-O-O angle in the force field training as mentioned in eq 6. However, making it more accurate for this angle strain curve resulted in loss of accuracy of the transition state energy of the H atom of the HO2 radical switching from one O atom to another. Given our application target we deemed that this transition state energy was more important than the H-O-O angle equilibrium value; as such, we accepted this deviation between ReaxFF and QM.

ReaxFF

QM

H

51.54

51.89

O

59.19

58.43

OH

7.35

8.4

H2O

-51.59

-56.8

H2O2 HO2

-35.28 1.67

-29.9 3.13

Overall we can say that ReaxFF captures the trend of angle distortion well, especially near the equilibrium valence angle. 3.1.3. Heats of Formation. The standard heat of formation (at 0 K) values of key reaction components of hydrogen combustion were calculated from G3 energy values.44-47 These values were added in the training set for ReaxFF. Table 2 shows the comparison between QM and ReaxFF. The heats of formation of ReaxFF and QM match well with a maximum deviation of ∼6 kcal for hydrogen peroxide. 3.1.4. Heats of Reaction and Energy Barriers. We employed a list of key reactions associated with hydrogen combustion from Li et al.9 in the ReaxFF training set. This same list of reactions has been quoted as the important reactions in hydrogen combustion by Mueller et al.8 and earlier by Yetter et al.52 Table 3 shows the comparison of heats of reaction and energy barriers of these reactions between ReaxFF and as obtained from Li et al.9 and 963



ARTICLE

Table 3. Comparison of Standard Heats of Reaction and Reaction Energy Barriers (Ea) of Key Reactions of Hydrogen Combustion in Literature and ReaxFFa reactions

ΔHr°0K (ReaxFF)

ΔHr°298K (Baulch et al.)b

Ea (ReaxFF)

Ea (Li et al.)c

16.6

Chain Reactions 1

H þ O2 f O þ OH

2 3

O þ H2 f H þ OH H2 þ OH f H2O þ H

4

H2O þ O f OH þ OH

15.00

16.31

19.76

-0.30 -7.40

1.41 -14.56

6.72 0

7.11

15.95

11.81

13.4

104.38

6.29 3.43

H2/O2 Dissociation/Recombination Reactions 5

H2 þ M f H þ H þ M

103.08

104.11

105.28

6

O þ O þ M f O2 þ M

-118.37

-119.1d

7.91

0e

7

O þ H þ M f OH þ M

-103.38

-102.3

8

H þ OH þ M f H2O þ M

-110.48

-118.68

9 10

H þ O2 þ M f HO2 þ M HO2 þ H f H2 þ O2

-49.87 -53.21

-48.342 -55.562

0 5.36

11

HO2 þ H f OH þ OH

-38.51

-37.832

0

0.3

12

HO2 þ O f OH þ O2

-53.51

-53.362

0

0

13

HO2 þ OH f H2O þ O2

-60.62

-70.362

0

14

HO2 þ HO2 f H2O2 þ O2

-38.62

-39.052

15

OH þ OH þ M f H2O2 þ M

-49.98

-50.282

1.23

16

H2O2 þ H f H2O þ OH

-60.51

-68.412

0

3.97

17 18

H2O2 þ H f H2 þ HO2 H2O2 þ O f OH þ HO2

-14.59 -14.89

-16.032 -14.312

3.99 3.73

7.95 3.97

19

H2O2 þ OH f H2O þ HO2

-21.99

-30.332

0

d

0

0

0

0

0.82

Formation/Consumption HO2

Formation/Consumption H2O2

a

0

All energies are in kcal. b Ref 53. c Ref 9. d Refers to values taken from Mueller et al. (ref 8). e M 6¼ Ar, He.

Baulch et al.53 We observe an overall good accuracy between ReaxFF and literature values. The maximum deviation between ReaxFF and literature values is within 10 kcal/mol for heats of reaction and 8 kcal/mol for energy barriers. Hydrogen exchange reactions R5 and R6 were also added in the training set. ðR5Þ HA - HB þ HC f HA þ HB - HC H A - HB þ HC - HD f HA - H C þ HB - HD

with QM predictions, especially for the lowest energy path with a maximum deviation within 20 kcal/mol. A similar study was done for reaction R8, i.e. H2 þ O2 f HO2 þ H

ðR8Þ

A hydrogen molecule was made to approach an oxygen molecule from several angles (i.e., O-O-HH angle varying from 80° to 180°), and the reaction energy curves were calculated in QM for each angle by constraining the internal coordinates to keep the angle fixed. The QM results were added in the ReaxFF training set for optimization. Figure 5 shows the comparison between QM calculations and ReaxFF predictions of reaction energy barriers. ReaxFF predictions are matching well with QM predictions, especially for the lowest energy path. Maximum deviation was observed to be approximately 22 kcal/mol. A more comprehensive study was done for reaction R2. Reaction R2 involves the formation of an O-H bond and dissociation of an H-H bond. For this reaction we considered a 2D grid of bond lengths. Considering the H-H bond distance and O-H bond distance as the coordinates of the grid (Figure 6), we created geometries for each grid point and energy-minimized them in QM (by constraining the O-H and H-H bond lengths and optimizing the rest of the geometry). The left top grid point in Figure 6 represents the starting point of the reaction (H2 and O). This point is taken as the zero-energy point. The right bottom point of the grid is considered the end point of reaction (H and OH). Figure 6 shows the reaction energy comparison in QM and ReaxFF. We can see that ReaxFF

ðR6Þ

The energy barriers for reactions R5 and R6 are in good agreement with QM calculations. These results are further discussed in the Supporting Information of this manuscript. Besides predicting the heat of formation, heat of reaction, and energy barriers with good accuracy we also need to predict the whole potential energy surface for a reaction with good accuracy. This is necessary because we are simulating the dynamics of a system, where molecules in a reaction do not always approach along the lowest energy path. For reaction R7, i.e. ðR7Þ H2 O þ O f 2OH we calculated the reaction energy barrier as a function of the angle of approach between O and H2O. We made the O atom approach the water molecule from several angles (i.e., O-H-OH angle varying from 80 to 180°) and calculated the reaction energy curves for each angle by constraining the internal coordinates to keep the angle fixed. The QM predictions were added to the training set for ReaxFF. Figure 4 shows that the ReaxFF predictions of the reaction energy barrier match well 964



ARTICLE

the QM and ReaxFF plots show the corresponding lowest energy paths observed for the reaction. QM predicts a barrier of 9.5 kcal/ mol along the lowest energy path, whereas ReaxFF predicts a barrier of 6.7 kcal/mol. 3.1.5. Force Field Validation from MD Results. We conducted an NVT MD simulation (N = 200 (φ = 1), V = 25 Å 25 Å 25 Å, T = 3500 K) using the optimized ReaxFF force field. From this MD simulation trajectory, we randomly picked molecule geometries of key reaction components of hydrogen combustion (i.e., H2, O2, H2O, OH, and HO2) from various instants of the simulation. Using these geometries we calculated their standard heats of formation (at 0 K) in ReaxFF and QM based on their single-point energies. Figure 7 shows a comparison between QM and ReaxFF results. We observed that the heat of formation of ground-state molecules like H2, O2, and H2O are predicted with good accuracy, compared to QM, whereas there is more deviation in case of reaction intermediates like OH and HO2 radicals. This is likely because our training data chiefly consists of bond stretching and angle distortion data related to ground-state molecules, whereas the reaction intermediates are only described in the transition state data. This analysis enables us to identify areas where the force field is less accurate, which enables us to further improve the quality of the force field in the future by adding corresponding structures in the training set with QM-calculated values. The final ReaxFF force field parameters can be found in the Supporting Information of this manuscript. A detailed description of the parameters and the corresponding potential functions can be found in the supporting material of Chenoweth et al.36 3.2. MD Simulations of Hydrogen Combustion 3.2.1. Water Formation Rate as a Function of Temperature, Volume, and Mixture Composition. We analyzed the reactivity of a H2/O2 mixture through NVT MD simulations by observing the water formation as a function of temperature, volume, and mixture composition of the simulation. Figure 8 shows the number of water molecules formed against the simulation time for a single NVT MD simulation (N = 200 (67 H2, 33 O2), V = 25 Å 25 Å 25 Å, T = 3500 K). We can observe from the figure that the water formation occurs in steps. This is expected because we are simulating a very small system and reaction occurs through rare events. In order to obtain continuum type behavior we need to take an average of multiple simulations. We ran NVT MD simulations at (1) N = 200 (67 H2, 33 O2), V = 25 Å 25 Å 25 Å, T = 3000-4000 K, (2) N = 200 (67 H2, 33 O2), V = 20 Å 20 Å 20 Å to 40 Å 40 Å 40 Å, T = 3500 K, and (3) N = 200 (φ = 0.25-4.5), V = 25 Å 25 Å 25 Å, T = 3500 K to study the effect of temperature, volume, and mixture composition, respectively. All the simulation were run for 0.5 ns. Figure 9 shows the water formation profile obtained as a function of temperature, volume, and mixture composition. Each plot in this figure has been obtained from an average of five independent simulations. We observe that, in general, the higher the temperature is, the earlier is the initiation and higher reaction rate is observed. Similarly, the smaller the volume (i.e., higher density), the earlier is the reaction initiation and higher is the reaction rate. The stoichiometric mixture (φ = 1) has the fastest reaction kinetics, and reaction is faster for fuel-lean mixtures as opposed to fuelrich mixtures, until the fuel-lean mixtures exhaust the entire fuel (for example, observe φ = 0.5 and φ = 2.0 in Figure 9). 3.2.2. Initiation Reaction. For all our input conditions (Table 1) we observed that the initiation reaction is predominantly

Figure 4. QM(B3LYP/6-311G**) (top) and ReaxFF (bottom) energies for distortion of the O-H-O angle when an O radical approaches the H atom in a H2O molecule.

Figure 5. QM(B3LYP/6-311G**) (top) and ReaxFF (bottom) energy for distortion of the H-O-O angle when a H2 molecule approaches an O2 molecule.

predicts the whole potential energy surface for reaction R2 with reasonably good accuracy as compared to QM. The black lines in

H2 þ O2 þ M f H þ HO2 þ M 965

ðR9Þ



ARTICLE

Figure 6. ReaxFF and QM(B3LYP,6-311G**) energies (in kcal) for the reaction path R2, i.e., H2 þ O f OH þ H. The figure represents a 2D grid of H-H and O-H bond lengths in reaction R2 with energy calculated at each grid point by constraining the corresponding bond lengths and optimizing the rest of the geometry. The top left point is the start point of the reaction, and bottom right is the end point. The black lines on the contours represent the lowest energy paths.

Figure 7. QM(B3LYP,6-311G**) and ReaxFF calculated standard heats of formation at 0 K for molecules obtained randomly from an NVT MD simulation run for 0.5 ns at 3500 K in a cubic periodic box of side length 25 Å containing 67 H2 and 33 O2 molecules (F = 0.13 kg/dm3).

However, at very high temperature (4000 K) and low density conditions (0.03 kg/dm3) we sometimes observe the initiation reaction as reaction R10, whereas at high density (0.25 kg/dm3) and low temperature (3000 K) we sometimes observe reaction R11. ðR10Þ H2 þ M f H þ H þ M H2 þ O2 þ O2 f 2HO2

before it finally reacts with one O2 molecule. The process is shown illustratively for one of the NVT MD simulations in Figure 10. It must be noted that not all collisions necessarily lead to a reaction. Reaction R9, without third-body dependence, has been suggested as the initiation reaction by earlier studies.54 3.2.3. Reaction Mechanism. The reaction pathways observed at our simulation conditions (Table 1) have predominantly consisted of reactions involving HO2 and H2O2. However, the percentage of chain-branching reactions R1 and R2 was found to increase when we increase the temperature or increase the volume of the simulation box. The reaction intermediates observed at our input conditions are primarily H, OH, HO2, O, and H2O2 molecules before leading to water as the final product.

ðR11Þ

In order to evaluate the role of the third body M in reaction R9, we investigated the trajectory of the H2 molecule involved in the reaction. We observed that the H2 molecule experiences several collisions with other molecules in the mixture (other H2 and O2 molecules in our case), resulting in vigorous vibrations and rotations 966



ARTICLE

Figure 8. Number of water molecules formed as a function of time for an NVT MD simulation run for 0.5 ns.

Figure 9. Water formation vs time for varying temperatures (volume and mixture composition fixed), varying volumes (temperature and mixture composition fixed), and varying mixture compositions (temperature and volume fixed). The results were obtained from an average of five independent simulations in each case.

The reaction pathway is shown illustratively for one of the NVT MD simulations (N = 200 (φ = 1), V = 25 Å 25 Å 25 Å, T = 3500 K) in Figure 11. After the initiation reaction, the reaction pathway primarily consists of consumption of H radicals by O2 molecules in the mixture to form HO2 radicals. The HO2 radicals then react to get converted to OH radicals, either directly or via H2O2 formation. And finally the OH radicals consume H2 molecules to form water and H radicals which continue the chain of reactions. Mueller et al.8 predicted a similar reaction pathway for input conditions lying above the extended second explosion limit. 3.2.4. Extended Second Explosion Limit at Our Input Conditions. In order to gain perspective of our input conditions w.r.t. the classical explosion limit diagram for mixtures of hydrogen/oxygen,

we analytically derived the extended second explosion limit using the steady-state approximation for radicals as shown in eq 7. 2k1 P ð7Þ ¼ ½M ¼ Ru T k12 k1 and k12 are rate constants of reactions R1 and R12. i.e. H þ O2 f OH þ O

ðR1Þ

H þ O2 þ M f HO2 þ M

ðR12Þ

and have been obtained from Li et al.9 [M] is the concentration of the third body, which in our case is the whole system and is expressed in terms of pressure P and temperature T using the ideal 967



ARTICLE

Figure 10. H-H bond distance for the H2 molecule in reaction with O2 at T = 3000 K, just before the reaction occurs. We observe collision of the H2 molecule with a third body M (= O2 in this case) resulting in an increase in amplitude of the H-H bond vibration and subsequent reaction.

Figure 12. Extended second explosion limits calculated using eq 7, assuming third-body (M) collisions only from H2, O2, and H2O, respectively. The collision efficiency values were taken from Li et al. (ref 9).

gas equation. The extended second limit is derived for three cases, i.e., (1) hydrogen is taken as third body M (fuel-rich system), (2) oxygen is taken as third body M (fuel-lean system), (3) water is taken as third body M (system at high conversion). The extended second limit derived is shown in Figure 12 and shows the input conditions of our study on the figure. For comparison sake we have also included input conditions of an earlier study done using experiments and continuum-scale computational methods. We observe from the figure that our input conditions lie close to the analytically derived extended second explosion limit, and we could very well observe a transition of reaction mechanism from the HO2/H2O2 reaction pathway to chain-branching reactions involving OH/O radicals in reactions R1 and R2.

Figure 11. O2 reaction pathway, for the input conditions N = 200 (φ = 1), V = 25 Å 25 Å 25 Å, T = 3500 K, over a simulation time of 0.5 ns. The thickness of the arrows is proportional to the number of molecules consumed of the species at the rear end of the arrow to form the front end species (normalized by the OH f H2O case). The table below lists the reactions involved in the paths shown along with their percent contribution to the destruction of a given species. 968



ARTICLE

Figure 13. H radical formation vs time for varying volumes of a periodic box. Temperature (3500 K) and mixture composition (N = 200, φ = 1) are kept fixed. The results were obtained from an average of five independent simulations at each volume.

3.2.5. Transition of Reaction Kinetics. In order to better understand the reaction mechanism, we analyzed the formation of intermediate radicals as a function of varying volume, similar to that done for water formation in section 3.2.1. Figure 13 shows the number of H radicals formed against simulation time as a function of density. We can observe from the figure that, at higher densities, the concentration of H radicals gradually increases and then gets consumed, whereas at lower densities we observe sharp peaks in H radical concentration. This indicates an increase in branching factor as we move toward lower densities, leading to formation of a large number of H radicals at a time. This behavior implies a transition in reaction mechanism from straight-chain to chain-branching reactions as we move from higher densities toward lower densities. In order to explore this further we compared the reaction mechanism of two extreme points from our input conditions, i.e., (1) temperature = 3000 K, density = 0.25 kg/dm3 (high pressure, low temperature) and (2) temperature = 4000 K, density = 0.03 kg/dm3 (low pressure, high temperature). We observed that, at high pressure and low temperature, the reaction pathway is predominantly via formation and consumption of HO2 and H2O2. At low pressure and high temperature, we observe another reaction pathway involving chain-branching reactions R1 and R2. However, the HO2/H2O2 reaction pathway is predominant in both cases. This is shown illustratively in Figure 14. This suggests that, even at our input conditions, the extended second explosion limit remains an important boundary separating regions of different reaction mechanisms. 3.2.6. Overall Reaction Energy Barrier of the Hydrogen Oxidation Reaction. The overall water formation reaction is simply ðR13Þ H2 þ 1=2O2 f H2 O

Figure 14. Transition of reaction kinetics from T = 3000 K, density = 0.25 kg/dm3 (top) to T = 4000 K, density = 0.03 kg/dm3 (bottom).

where [H2] and [O2] are the average concentrations of H2 and O2, respectively, obtained from five different simulations. d[H2]/dt is calculated by fitting a curve to the [H2] concentration profile and evaluating its slope. The rate constant is evaluated for each time point in the simulation, and an average is taken over the whole simulation time. The reaction energy barrier is then calculated using the Arrhenius rate equation (eq 9), i.e.

NVT MD simulations of H2/O2 mixtures were run at two different densities (N = 200 (φ = 1), V = 25 Å 25 Å 25 Å and V = 32 Å 32 Å 32 Å) for temperatures ranging from 3000 to 4000 K. For each temperature, five simulations with different initial configurations were run for 0.5 ns. The rate constant of the overall reaction R13 is calculated using eq 8, i.e. d½H2 1 pffiffiffiffiffiffiffiffi ð8Þ kov ¼ dt ½H2 ½O2

kov 969

Ea ¼ A exp RuT

ð9Þ



ARTICLE

Figure 16. Water formation vs time for NVE MD simulation of a hydrogen/oxygen mixture (66 H2 and 33 O2 molecules) seeded with one OH radical with initial temperature (T = 2000 K) and density (F = 0.25 kg/dm3). The figure also shows the consumption rate of H2 and O2 and temperature profile.

4. CONCLUSION We have described an application of ReaxFF for getting a more detailed understanding of the reaction kinetics of hydrogen combustion under high-pressure conditions. The ReaxFF force field parameters were optimized against a training set containing reaction energy and transition state data of important reactions of hydrogen combustion. Quantum mechanics derived bond stretching profiles, angle distortion profiles, potential energy surface, and heats of formation were also added to the training set. The optimized force field was used to perform a range of NVT MD and NVE MD simulations to elucidate the reaction kinetics of hydrogen combustion at our input conditions. Under all our input conditions (3000 K e T e4000 K, 0.25 e φ e 4.5, 0.03 kg/dm3 e density e 0.25 kg/dm3), the initiation reaction was found to be predominantly reaction R9, i.e.

Figure 15. Arrhenius plot for the overall rate constant derived from NVT MD simulations N = 200 (φ = 1), V = 32 Å 32 Å 32 Å (density = 0.06 kg/dm3), T = 3000-4000 K and N = 200 (φ = 1), V = 25 Å 25 Å 25 Å (density = 0.13 kg/dm3), T = 3000-4000 K. The black line in each case is a least-squares fit of the form k = A exp(-Ea/RT). From the plot we get a reaction energy barrier Ea = 55.73 ( 4.97 kcal/mol for density = 0.06 kg/dm3 and Ea = 34.26 ( 3.28 kcal/mol for density = 0.13 kg/dm3.

We obtained energy barrier values of 34.26 ( 3.28 kcal/mol for V = 25 Å 25 Å 25 Å and 55.73 ( 4.97 kcal/mol for V = 32 Å 32 Å 32 Å as shown in Figure 15. Mueller et al.8 had conducted a similar study for fuel-lean H2/O2/N2 mixtures (P = 6.5 atm, T = 884-934 K, φ = 0.3) and found the overall reaction energy barrier to be 61 ( 10 kcal/mol. We had no existing experimental/ continuum model data at our input conditions to compare our results with. However, the values obtained by us are higher than 16 kcal/mol, which is the overall barrier value for input conditions in the zone below the extended second explosion limit. This suggests that our input conditions do not lie completely in the zone below the extended second explosion limit. 3.2.7. Mixture Seeded with One OH Radical. We performed an NVE MD simulation for a H2/O2 mixture seeded with one OH radical (N = 200 (66 H2, 33 O2, 1 OH), V = 20 Å 20 Å 20 Å). Starting with a temperature of 2000 K, we performed the simulation for a period of 0.5 ns. Figure 16 shows the rate of water formation against simulation time for this case. We observed that water formation occurs very early in the simulation due to the presence of OH radical in the mixture, by reaction R14, i.e. H2 þ OH f H2 O þ H

H2 þ O2 þ MðH2 =O2 Þ f H þ HO2 þ M

ðR9Þ

As opposed to earlier literature,8,9,54 we have observed a thirdbody effect in reaction R9. The rate of water formation for a fixed density and fixed equivalence ratio was found to increase with temperature. Similarly for fixed temperature and equivalence ratio, rate of water formation increased with increase in mixture density. The rate of water formation was found to be maximum for a stoichiometric mixture of hydrogen and oxygen, at a fixed temperature and fixed density. It decreased for both fuel-lean and fuel-rich mixtures, with fuelrich mixtures having a lower rate of formation. Reaction initiation was found to be faster for higher temperature and higher density. The intermediate reactions were found to be dominantly consisting of formation and consumption reactions of HO2 and H2O2. Most of the oxidizer O2 gets converted into HO2 radical instead of choosing the chain-branching reaction pathways R1 and R2. However, a transition in reaction mechanism from the HO2/H2O2 pathway to chain-branching reactions R1 and R2 was observed with increasing temperature and decreasing density. Earlier experimental and continuum-scale kinetic studies as summarized by Mueller et al.8 have discovered these same intermediate reactions involving HO2/H2O2 for low-temperature and high-pressure input conditions and reactions involving

ðR14Þ

making this simulation equivalent to simulating combustion in a H2/O2 mixture after initiation. We observe that initially the reaction proceeds slowly, until around 200-250 ps there is a sudden rise in the water molecule concentration as well as the temperature and pressure of the mixture. This suggests the occurrence of explosion in the mixture. This type of analysis will be useful to study shock behavior in hydrogen combustion and is part of our future work. 970



ARTICLE

reactions R1 and R2 for high-temperature and low-pressure conditions. This suggests that the extended second explosion limit remains an important boundary even at our range of input conditions and there exist two distinct reaction mechanisms of water formation depending on the input pressure and temperature. The overall reaction energy barrier of hydrogen combustion found was much higher than 16 kcal/mol, which suggests that our input conditions do not lie in the chain-explosive region (i.e., the zone below the extended second explosion limit). From these results, we can say that ReaxFF can be a useful tool to elucidate complex reaction mechanisms and obtain reaction kinetic data. ReaxFF can be used to simulate reactions at a variety of input conditions, which are otherwise difficult to perform experimentally, and its results can be used to improve existing continuum-scale kinetic models. With our newly added reaction mechanism analysis tool, it is possible to study in greater detail the initiation and intermediate reactions.

(13) Slater, J. C.; Koster, G. F. Phys. Rev. 1954, 94, 1498. (14) Andersen, O. K.; Jepsen, O. Phys. Rev. Lett. 1984, 53, 2571. (15) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179–5197. (16) Cygan, R. T.; Liang, J. J.; Kalinichev, A. G. J. Phys. Chem. B 2004, 108, 1255–1266. (17) Rappe, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A.; Skiff, W. M. J. Am. Chem. Soc. 1992, 114, 10024–10035. (18) Mayo, S. L.; Olafson, B. D.; Goddard, W. A. J. Phys. Chem. 1990, 94, 8897–8909. (19) Allinger, N. L.; Yuh, Y. H.; Lii, J. H. J. Am. Chem. Soc. 1989, 111, 8551–8566. (20) Stuart, S. J.; Tutein, A. B.; Harrison, J. A. J. Chem. Phys. 2000, 112, 6472–6486. (21) Brenner, D. W. Phys. Rev. B 1990, 42, 9458–9471. (22) Polzer, T.; Kiefer, W. J. Organomet. Chem. 1996, 508, 153–157. (23) Tersoff, J. Phys. Rev. B 1988, 37, 6991–7000. (24) van Duin, A. C. T.; Dasgupta, S.; Lorant, F.; Goddard, W. A. J. Phys. Chem. A 2001, 105, 9396–9409. (25) Strachan, A.; van Duin, A. C. T.; Chakraborty, D.; Dasgupta, S.; Goddard, W. A. Phys. Rev. Lett. 2003, 91, 098301. (26) Strachan, A.; Kober, E. M.; van Duin, A. C. T.; Oxgaard, J.; Goddard, W. A. J. Chem. Phys. 2005, 122, 054502. (27) van Duin, A. C. T.; Zeiri, Y.; Dubnikova, F.; Kosloff, R.; Goddard, W. A. J. Am. Chem. Soc. 2005, 127, 11053–11062. (28) Chenoweth, K.; Cheung, S.; van Duin, A. C. T.; Goddard, W. A.; Kober, E. M. J. Am. Chem. Soc. 2005, 127, 7192–7202. (29) Goddard, W. A.; van Duin, A.; Chenoweth, K.; Cheng, M. J.; Pudar, S.; Oxgaard, J.; Merinov, B.; Jang, Y. H.; Persson, P. Top. Catal. 2006, 38, 93–103. (30) Goddard, W. A.; Merinov, B.; Van Duin, A.; Jacob, T.; Blanco, M.; Molinero, V.; Jang, S. S.; Jang, Y. H. Mol. Simul. 2006, 32, 251–268. (31) Buehler, M. J.; van Duin, A. C. T.; Goddard, W. A. Phys. Rev. Lett. 2006, 96, 095505. (32) Ludwig, J.; Vlachos, D. G.; van Duin, A. C. T.; Goddard, W. A. J. Phys. Chem. B 2006, 110, 4274–4282. (33) Cheung, S.; Deng, W.-Q.; van Duin, A. C. T.; Goddard, W. A. J. Phys. Chem. A 2005, 109, 851–859. (34) Nielson, K. D.; van Duin, A. C. T.; Oxgaard, J.; Deng, W. Q.; Goddard, W. A. J. Phys. Chem. A 2005, 109, 493–499. (35) Zhang, Q.; Cagin, T.; van Duin, A.; Goddard, W. A.; Qi, Y.; Hector, L. G. Phys. Rev. B 2004, 69, 045423. (36) Chenoweth, K.; van Duin, A. C. T.; Goddard, W. A. J. Phys. Chem. A 2008, 112, 1040–1053. (37) Mortier, W. J.; Ghosh, S. K.; Shankar, S. J. Am. Chem. Soc. 1986, 108, 4315–4320. (38) Rappe, A. K.; Goddard, W. A. J. Phys. Chem. 1991, 95, 3358– 3363. (39) van Duin, A. C. T.; Baas, J. M. A.; van De Graaf, B. J. Chem. Soc., Faraday Trans. 1994, 90, 2881–2895. (40) Becke, A. D. J. Chem. Phys. 1993, 98, 5648–5652. (41) Lee, C. T.; Yang, W. T.; Parr, R. G. Phys. Rev. B 1988, 37, 785– 789. (42) Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A. J. Chem. Phys. 1980, 72, 650–654. (43) Jaguar, 7.5 ed.; Schrodinger,LLC: New York, 2008. (44) Curtiss, L. A.; Raghavachari, K.; Redfern, P. C.; Rassolov, V.; Pople, J. A. J. Chem. Phys. 1998, 109, 7764–7776. (45) Curtiss, L. A.; Raghavachari, K.; Redfern, P. C.; Pople, J. A. J. Chem. Phys. 2000, 112, 7374–7383. (46) Curtiss, L. A. G3 Heats of Formation. http://www.cse.anl.gov/ OldCHMwebsiteContent/compmat/g3energies/g3neut.htm (accessed Feb 1, 2009). (47) Curtiss, L. A. G3 Energies of Auxiliary Neutral Species. http:// www.cse.anl.gov/OldCHMwebsiteContent/compmat/g3energies/ G3aux.htm (accessed Feb 1, 2009).

’ ASSOCIATED CONTENT

bS Supporting Information. Optimized force field, which can be run on the standalone ReaxFF program (available from the authors on request) and requires modification in the lone pair potential function to be run with LAMMPS, as indicated in this manuscript, and figures showing the energy barrier prediction of ReaxFF and QM for hydrogen exchange reactions R5 and R6. This material is available free of charge via the Internet at http:// pubs.acs.org. ’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

’ ACKNOWLEDGMENT This research was supported by KISK Startup Grant No. C000032472. We also acknowledge valuable discussions with Professor Richard Yetter of Mechanical and Nuclear Engineering at Penn State University. ’ REFERENCES (1) Midilli, A.; Dincer, I. Int. J. Hydrogen Energy 2008, 33, 4209– 4222. (2) Marban, G.; Vales-Solis, T. Int. J. Hydrogen Energy 2007, 32, 1625–1637. (3) Schlapbach, L.; Zuttel, A. Nature 2001, 414, 353–358. (4) Palumbo, O.; Paolone, A.; Rispoli, P.; Cantelli, R. Mater. Sci. Eng., A 2009, 521-522, 134–138. (5) Petukhov, V. A.; Naboko, I. M.; Fortov, V. E. Int. J. Hydrogen Energy 2009, 34, 5924–5931. (6) Lewis, B.; von Elbe, G. Combustion,Flames and Explosions of Gases, 3rd ed.; Academic Press: Orlando, FL, 1987. (7) Dixon-Lewis, G.; Williams, D. J. Compr. Chem. Kinet. 1977, 17, 1-248. (8) Mueller, M. A.; Kim, T. J.; Yetter, R. A.; Dryer, F. L. Int. J. Chem. Kinet. 1999, 31, 113–125. (9) Li, J.; Zhao, Z.; Kazakov, A.; Dryer, F. L. Int. J. Chem. Kinet. 2004, 36, 566–575. (10) Strohle, J.; Myhrvold, T. Int. J. Hydrogen Energy 2007, 32, 125– 135. (11) O’Conaire, M.; Curran, H. J.; Simmie, J. M.; Pitz, W. J.; Westbrook, C. K. Int. J. Chem. Kinet. 2004, 36, 603–622. (12) Raty, J. Y.; Gygi, F.; Galli, G. Phys. Rev. Lett. 2005, 95, 096103. 971



ARTICLE

(48) Berendsen, H. J. C.; Postma, J. P. M.; Vangunsteren, W. F.; Dinola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684–3690. (49) Shanno, D. F. Math. Oper. Res. 1978, 3, 244–256. (50) Verlet, L. Phys. Rev. 1967, 159, 98 ff. (51) Swope, W. C.; Andersen, H. C.; Berens, P. H.; Wilson, K. R. J. Chem. Phys. 1982, 76, 637–649. (52) Yetter, R. A.; Rabitz, H.; Hedges, R. M. Int. J. Chem. Kinet. 1991, 23, 251–278. (53) Baulch, D. L.; Bowman, C. T.; Cobos, C. J.; Cox, R. A.; Just, T.; Kerr, J. A.; Pilling, M. J.; Stocker, D.; Troe, J.; Tsang, W.; Walker, R. W.; Warnatz, J. J. Phys. Chem. Ref. Data 2005, 34, 757–1397. (54) Westbrook, C. K. Combust. Sci. Technol. 1982, 29, 67–81.

972


Development and Application of a ReaxFF Reactive Force Field for

Recommend Documents