- Email: [email protected]
Exothermic and endothermic chemical reactions involving very many particles modeled with molecular dynamics
Physica D 146 (2000) 261–274
Exothermic and endothermic chemical reactions involving very many particles modeled with molecular dynamics Witold Alda a , David A. Yuen b,∗ , Hans-Peter Lüthi c , James R. Rustad d a
Institute of Computer Science, University of Mining and Metallurgy AGH, 30-059 Krakow, Poland b Department of Geology and Geophysics and Minnesota Supercomputing Institute, University of Minnesota, 1200 Washington Ave. South, Minneapolis, MN 55415-1227, USA c Swiss Center for Scientific Computing, ETH, Zurich CH-8092, Switzerland d The Pacific Northwest National Lab, W.R. Wiley EMSL, Richland, WA 99352, USA Received 5 August 1999; accepted 3 July 2000 Communicated by R.P. Behringer
Abstract The traditional continuum approach of modeling chemical reactions with specified kinetic rates suffers from numerical difficulties in reactive flows and other highly non-equilibrium situations due to the stiffness of the differential equations in both space and time. These drawbacks can be eliminated within the framework of the discrete-particle approach in which the chemical reactions are modeled by means of two-body interparticle potentials with classical molecular dynamics. Here we present a simple model for prescribing the binary endothermic and exothermic reactions of the type A + B → C, in the presence of many (greater than 104 ) reacting particles. Two-body Lennard–Jones based potentials, Vij , have been utilized, in which the reaction takes place via an attractive potential VAB with a deep enough potential well. The other potentials VAC and VBC have deeper wells than VAB for creating a suitable situation for inducing exothermic reactions, while a high-energy plateau for VAC and VBC is used for simulating endothermic reactions. We have deployed 20 000 to 100 000 particles in a computational domain with an area of around 1 million Å2 . These ensembles have all been integrated out to around the order of nanoseconds. We have examined the bimolecular reactions between two layers initially consisting of A particles lying atop B particles, equally divided in the total population. Reactions take place first along this boundary. There is a nonlinear threshold phenomenon in the exothermic reactions, which are greatly accelerated by the heat liberated from the neighboring reactions. In contrast, the endothermic reactions are more subdued, as there is no positive feedback. Exothermic reactions occur in a spatially heterogeneous manner, especially at the lower temperatures, in which no single reaction rate can be assigned to the entire reaction front. More coherent spatial structures are found with decreasing temperature. These spatial structures can influence the overall flow field associated with the motions of the reacting particles. Timescales of the reactions are greatly lengthened by a reduction in the temperature in a threshold manner. © 2000 Published by Elsevier Science B.V. Keywords: Exothermic and endothermic reactions; Molecular dynamics; Lennard–Jones potential
1. Introduction
∗ Corresponding author. Tel.: +1-612-625-1818; fax: +1-612-624-8861 E-mail address: [email protected] (D.A. Yuen).
Today chemical reacting flows form indeed important and fundamental foci for many scientific and engineering endeavors. Numerous examples can be drawn from fields as diverse as metamorphic reacting fronts
0167-2789/00/$ – see front matter © 2000 Published by Elsevier Science B.V. PII: S 0 1 6 7 - 2 7 8 9 ( 0 0 ) 0 0 1 5 0 - 0
262
W. Alda et al. / Physica D 146 (2000) 261–274
imprinted in Alpine rocks and the ionization processes accompanying the reentry of hypersonic vehicles descending back into the earth’s atmosphere. Most of the emphasis on the study of chemically reacting flows has been directed toward the continuum approach based on solving nonlinear partial differential equations with the chemical reactions added in as a set of nonlinear ordinary differential equations in time. These chemical conservation relations are generally prescribed by given kinetic-rate laws, which are thermally activated, and hence strongly temperature-dependent in the Arrhenius sense. However, solving these coupled species-energy-momentum equations is beset by the multiple-scale nature of the coupled equations, which baffle even the best of numerical methods, especially at the chemically reacting fronts with sharp gradients. The stiffness of the partial differential equations arising from the use of both the spatial and temporal operators in the continuum approach can be ameliorated by using the discrete methods, such as cellular automata, lattice-gas, lattice-Boltzmann and molecular dynamics (MD), as they operate over much shorter timescales, many orders in magnitude smaller than the time steps taken in the continuum approach. The above methods are based on very simplified models and often they do not give accurate or even realistic results, but according to Binder et al. [1], “. . . the goal of such simulations is not to predict numbers but to provide better understanding”. In our case, this means a better understanding of the dynamics of chemical reaction processes with many particles in highly non-equilibrium situations. Cellular-automaton methods have been used to study diffusion and reaction processes by Mai and von Niessen [2], Chopard and Drozd [3] as well as Weimar and Boon [4]. There have been some previous modeling of chemical reactions at mineral surfaces by lattice gas automata (LGA) model by Boon et al. [5], Wells et al. [6] as well as Rothman and Zaléski [7]. They found that LGA calculations appeared to be well suited for problems in the mesoscopic scale as in the range of pore-fluid domain, of the order of tens of microns. Classical MD along with simple two-body potentials has been used to study phase transitions in
two-dimensional Lennard–Jones (LJ) binary systems by Li et al. [8]. The efficacy of this simplified approach for long-time integration, involving 1000 to 10 000 particles for over 107 time steps, was emphasized because of the need to monitor the onset of melting, a relatively long timescale phenomenon. Chemical reaction simulations on the atomic level have been performed by Wilson and coworkers using classical MD (e.g. [9]) or transition state theory [10]. Stillinger and Weber [11] have used two-body and three-body potentials in MD simulations to observe spontaneous forming of diatomic molecules. Ortoleva and Yip [12] used a hard disks system for binary reactions, where particles react by collisions (forward reactions) and by spontaneous transitions (backward reactions), giving different steady states depending on the reaction rates. In contrast to our approach, all these simulations were performed in equilibrium systems with exchange reactions only and constant number of molecules. On the other hand, direct integration of the coupled quantum-mechanical self-consistent molecular dynamics (SCMD) problem [13–15] for several hundred atoms and for timescales, up to the nanosecond range, remains somewhat out of reach with the present generation of computers. Thus, the time is now for some middle-ground solution for intermediate timescales by means of MD, which can address this problem of bimolecular reactions with a large swarm of heterogeneous atoms which can display their collective fluid-like nature [16]. Some macroscopic fluid-dynamical phenomenon can be manifested, as in the MD study of the fluid-dynamical instabilities by Mo´sci´nski et al. [17], Alda et al. [18] and Dzwinel et al. [19], in which as many as 107 particles were integrated on massively parallel computers for around 106 time steps, in order to simulate mesoscale fluid-dynamical instabilities, such as the Rayleigh–Taylor instability. In this same spirit, we present in this paper a simpler conceptual and computational model for investigating binary chemical reactions, based on classical MD and two-body interparticle potentials. With the MD simulations, we can systematically investigate the essential characteristics of chemically reacting hydrodynamics with far more participating atoms, exceeding by
W. Alda et al. / Physica D 146 (2000) 261–274
several orders in magnitude what would be feasible by the SCMD technique, and for much longer timescales. The physical scale in which we will be operating still falls within the microscale range of under 1 m but the spatial range can be extended to the mesoscale domain of Wells et al. [6] by combining the dissipative particle dynamics (DPD) method [20,21], where particles represent clusters of atoms, with classical atomic MD particles. In the above approach, one can think of a system built up of large DPD particles representing non-reacting solvent with dissolved reacting MD atoms. The method described here will potentially constitute a powerful new tool for the computational treatment of chemical reacting hydrodynamics in the microscale and also eventually be connected to the mesoscale of the order of microns by DPD methods [22]. Our approach stands in sharp contrast to MD studies in which the chemical kinetics are assumed a priori for the reacting particles [23].
2. MD model Our simulations have been carried out in two dimensions. We have employed a two-body LJ-based interparticle potential, for both exothermic and endothermic reactions, modified for certain interactions: " # σαβ 6 σαβ 12 LJ − , (1) Vαβ (rij ) = 4αβ rij rij where α and β denote the particle species, αβ the well depth of the potential, σαβ a parameter determining the nearest neighbor distance in the interaction, and rij the distance between the particles i and j . All of the potentials are truncated at certain rcut and shifted up, so as to reach zero potential at cut-off and achieve potential continuity at the cut-off radius rcut . In this model, three different types of particles, A, B and C, which may differ in mass and interaction potential parameters, are employed. All interactions are assumed to have identical σαβ = σ parameters; thus all particles have the same size. The total number of particles in the system changes, as the reactions develop. The VAB potential for both types of reactions has been modified in this study by including a hump
263
(Figs. 1a and 2a) in the overall potential, which is given by LJ VAB (rij ) = VAB + Vhump .
(2)
The hump is modeled by piecing together two cubic LJ (r ), where r = ( 26 )1/6 σ , splines, which start at VAB sh sh 7 go through the top of the hump at rh and end at LJ (r ), where r = 2r − r . V VAB AB is truncated at eh eh h sh rcut = 2.5σ ; thus the hump lies within the interaction range. The idea of using this artificial feature is to enable the reactions to take place only for particles which have a sufficiently high kinetic energy. The VAA , VBB and VCC potentials, which are also identical for both the exothermic and endothermic reactions (Figs. 1c and 2c), are truncated at rcut = LJ , and shifted up by , 1.12σαβ , a minimum of Vαβ αβ thus making the potential purely repulsive in this study. The VBC and VAC potentials differ for both types of reactions. For exothermic reactions, the potential well depth, , in VBC and VAC is deeper than in VAB (Fig. 1b). This results in smaller average potential energy of A–C and B–C interactions than A–B that occurred before the reaction. Assuming that the total energy of the system is kept constant, the average kinetic energy increases. Changing the shape of the VAC and VBC functions by moving them up, increases the average A–C and B–C potential energy, eventually allowing the endothermic reactions to occur. For this purpose, we have chosen a very simple criterion, by introducing the repulsive potential with a flat plateau (Fig. 2b), where particles in A–C and B–C interactions can lose their kinetic energy. In this way it is possible to find the A–C and B–C potential for which there will be zero heat generated by the reaction. With A–B and A–A, B–B, C–C potentials used, this can be achieved with a very shallow potential well depth, = 20 J/K, in the A–C and B–C interactions. In addition to the exothermic and endothermic runs, we have carried out a series of zero heat runs at different temperatures in order to demonstrate the temperature dependence of the reaction speed. These results are shown at the end of this paper. The parameters of the exothermic and endothermic potentials are listed in Table 1.
264
W. Alda et al. / Physica D 146 (2000) 261–274
Fig. 1. The form of the three different kinds of energy potentials used in producing the binary exothermic reaction: (a) VAB represents the potential which produces the reaction A + B → C. Reaction takes place when both A and B remain together within one time step inside the well. The hump in front of the well controls the rate of the reaction; (b) VCB and VCA have deeper potential than VAB to enable exothermic nature of the reaction; (c) VAA , VBB , VCC are repulsive in nature. Thus, there are no clusters involving the formed A, B and C particles.
Fig. 2. The form of the three different kinds of potentials used in creating the binary endothermic reaction: (a) VAB represents the potential which produces the reaction A + B → C, identical to that in Fig. 1; (b) VCB and VCA have the plateau on which the endothermic reaction takes place. (c) VAA , VBB , VCC are repulsive in nature, identical to those in Fig. 1.
W. Alda et al. / Physica D 146 (2000) 261–274
265
Table 1 Potential parameters used in the three-component system with exothermic and endothermic reactions Interaction
αβ (J/kB )
σαβ (Å)
rcut (σαβ /Å)
Plateau cut-off (σαβ /Å)
Potentials equal for exothermic and endothermic reactions A–B 119.8 3.4 2.5 – A–A 119.8 3.4 1.12 – B–B 119.8 3.4 1.12 – C–C 119.8 3.4 1.12 – Exothermic only C–B 239.6 3.4 2.5 – C–A 239.6 3.4 2.5 – Endothermic only C–B 119.6 3.4 1.12 0.9 C–A 119.6 3.4 1.12 0.9
The mass of A and B particles has been set to 40 amu, while the mass of C particle, a product of A + B reaction, to 80 amu. The particles are confined within a rectangular box, and initially are placed on a grid forming a mesh built up from 2D regular squares. With the value of σ assumed, the size of such a cell is about 3.53 Å × 3.53 Å, which yields a total area of the computational box greater than 106 Å2 (Fig. 3) for a simulation with 105 particles. The initial configuration is not close packed, so the particles can diffuse more easily. This involves the initial stage of the
Hump height (J/kB )
Hump distance (σαβ /Å)
20.0 – – –
1.6 – – –
– –
– –
– –
– –
simulation, as the particles can start reacting earlier. In the course of simulation the structure is not kept and the reacting particles can form a collective fluid-like particle flow. The size of the box is kept constant. Periodic boundary conditions have been applied along the vertical edges, while for the horizontal boundary we have introduced reflecting walls. The time step, 1t, has been set to 10−14 s, which puts an upper-bound on the relative velocity of the interacting particles. With this model, we can achieve a speed of around two time steps per second for 105
Fig. 3. Schematic diagram of the computational domain, where the exothermic reactions take place (see colliding particles depicted). The initial grid layout with the regular square configuration of particles A and B is displayed.
266
W. Alda et al. / Physica D 146 (2000) 261–274
particles on a 600 MHz workstation. We can then complete 106 time step run (or 10 nanoseconds with the given 1t) within 1 day.
3. Algorithms for inducing binary endothermic and exothermic reactions For simulating chemical reactions with both endothermic and exothermic nature in MD, we have implemented simple algorithms to a model of A + B → C chemical reactions to a classical MD scheme based on LJ particles [24]. By a reaction, we mean here an “instantaneous” or an extremely short time-instant involving the transformation of particles A and B into a single particle C, when certain conditions concerning the duration of time spent inside the well of VAB and velocity magnitudes of particles are fulfilled. The occurrence and reaction speed depend on the four criteria: 1. the potential depths and widths applied, 2. the average energies or speed of the ensemble of particles, 3. the global thermodynamic conditions, such as the density and the temperature, and 4. the initial boundary conditions. In order for a reaction to occur between particles A and B, the VAB potential must have a potential well deep enough. The “hump” added in VAB can be especially important in low-density systems, where mean rAB goes beyond the hump. In such a case, the hump acts as a threshold. The non-reacting A–A, B–B and C–C pairs, interact by repulsive potential only (Fig. 1b). Repulsion prevents the forming of clusters of same-kind particles and enables better mixing. Particles A–C and B–C attract each other; however, the potential well depth is deeper than in the A–B interaction. Moreover, we do not allow C to react with A and B in this study, i.e. no back reaction. In this model, the principal criterion for the occurrence of a chemical reaction, is the interparticle distance. We assume that a reaction takes place for the condition rAB < rcrit , where rcrit is a certain critical distance. We also assume that the rAB < rcrit state should be kept for a certain number of time steps, con-
tinuously. The reaction rate can be slowed down then, by decreasing rcrit or increasing the time interval. In the simulations performed, we have assumed that the radius for reactive occurrence, rAB < rcrit , in a single time step, O(10−14 ) s, is sufficient for the reaction. In general, rcrit is chosen to lie within 0.8 . . . 0.9 of σ in A–B interaction. To avoid very long runs with very few reactions, we have chosen a relatively large value of rcrit = 0.87 to speed up the reactions. The model reactions do fulfill the conservation laws. Conservation of mass and momentum is straightforward, as mC = mA + mB and pEC = pEA + pEB . The center of mass of C is equal to the center of mass of A + B. The total energy, however, is not automatically conserved because each time the reaction occurs, two particles A and B are replaced with one particle C, which interacts with its nearest neighbors by a different potential. The particle C has no internal dynamical structure and cannot keep any hidden potential energy. To keep the total energy constant, we can only change the kinetic energy of some particles whenever a reaction takes place. We accomplish this task at each time step by adjusting the kinetic energy of all particles within 2.5σαβ radius around C, so as to conserve the total energy of the neighboring area at each reaction. 4. Results We have used a constant energy scaling in conducting the simulations of a system involving 20 000 particles for both exothermic and endothermic reactions for a range of temperatures in order to see some noticeable effects. The changes of potential, kinetic and total energies with time are displayed in Figs. 4 and 5. Exothermic runs have been performed in three initial dimensionless temperatures, T = 0.1, 0.3 and 0.7 in reduced LJ units. As shown in Fig. 4, for T = 0.1 and 0.7 runs, the total energy is conserved, while the potential and kinetic energies vary with time. The kinetic energy increases because of the exothermic nature of the reaction. Due to the increase in the number of particles C and mixing at the interface area, the potential energy now has a tendency to decrease, while the temperature of the entire system increases. A similar, but reversed situation is shown in Fig. 5,
W. Alda et al. / Physica D 146 (2000) 261–274
Fig. 4. Conservation of the total energy with time for exothermic reactions. Dimensionless initial temperature is 0.1–0.5 by the end of the run. 20 000 particles consisting of A and B initially were employed. The total number of particles changes in time due to the reactions.
in which with endothermic reactions the potential energy increases, while the temperature decreases monotonically. Figs. 6–8 show the snapshots of the particle locations for the exothermic reactions in three different initial temperatures, T = 0.7, 0.3 and 0.1, and in three different instants of time. Gray areas represent atoms A and B in the lower and upper parts of the box, while
Fig. 5. Conservation of the total energy with time for endothermic reactions. Dimensionless initial temperature is 0.2 decreasing to 0.08 by the end of the run. 20 000 particles consisting of A and B initially were employed. The total number of particles varies with time due to the chemical reactions.
267
the atoms C, the product of the reaction, are shown as darker dots. The latter are displayed in larger format for better visualization. Initially the system consists only of atoms A and B occupying the upper and lower halves of the box, separated by a straight line of contact. It can be observed from the snapshots, that in higher temperature the reaction starts much faster and quickly spreads over the entire interface. At the highest temperature (T = 0.7) the particle motion is more chaotic and the reactions take place along the entire front. On the other hand, reactions at T = 0.1 are sporadic and induce a heterogeneous character to the interface. This is due to the exothermic character of the reactions and the fact that, once the reaction occurs, the next reaction becomes more probable to occur in the neighborhood. Timescales are faster at the higher temperatures. One can also see that the reaction front forms a curve line, which is rugged but its shape remains periodic with a dominant wavelength. The curvature of the reaction interface can be explained by the exothermic character of the reactions and the tendency of the heat liberated to stretch the interface boundary, which is energetically adjusted. On the other hand, the reactions rate is slowed down with time due to the growing thickness of C particles layer which separates the reacting A and B species. In these three runs, starting from T = 0.7, 0.3 and 0.1, the temperatures rise asymptotically to T = 0.98, 0.60 and 0.42, respectively. Similar snapshots showing the endothermic reactions with the initial temperatures, T = 0.6, 0.4 and 0.2, are shown in Figs. 9–11. We observe that the image of the initial phase of the reactions in Figs. 9 and 10 is similar to that in the exothermic model. Furthermore, the interface line remains rather straight, with only short-wavelength fluctuations. Particles C are more scattered (especially in Fig. 10) due to the smaller probability of the next reaction in the vicinity of an already reacted site. A larger domain has been established by elongating the computational domain by fivefold. Two simulations for both exothermic and endothermic reactions have been performed for the larger box in the initial temperature T = 0.3. We observe in Fig. 12, that the spatial heterogeneity of the reactions is much
268
W. Alda et al. / Physica D 146 (2000) 261–274
Fig. 6. Global views of the particles for exothermic reactions. Darker dots in the central part denote C, while A and B are placed in the upper and lower halves of the box, respectively. Dimensionless initial temperature is set to 0.7.
Fig. 7. Global snapshots of the exothermically reacting particles for dimensionless initial temperature of T = 0.3.
Fig. 8. Global snapshots of the exothermically reacting particles for dimensionless initial temperature of T = 0.1.
W. Alda et al. / Physica D 146 (2000) 261–274
269
Fig. 9. Global views of the particle distribution for endothermic reactions. Darker dots in the central part denote C, while A and B are placed in the upper and lower halves of the box, respectively. The dimensionless initial temperature is set to 0.6.
Fig. 10. Global snapshots of the endothermically reacting particles for dimensionless initial temperature of T = 0.4.
Fig. 11. Global snapshots of the endothermically reacting particles for a dimensionless initial temperature of T = 0.2.
270
W. Alda et al. / Physica D 146 (2000) 261–274
Fig. 12. Global views of run with 100 000 initial particles for exothermic reactions. Dimensionless temperature is set to 0.3. The configuration has an aspect ratio of 5 and an area of 500 Å × 2500 Å.
stronger than in the smaller boxes. This tendency is probably due to greater length of the reaction front and smaller influence of periodic boundary conditions. The occurrence of the reaction generates heat, which is then distributed over the neighboring particles. The energy distribution along the interface is slower, when not constrained by the end conditions. This also shows that the speed and the range of the kinetic energy being distributed about, can greatly influence the overall dynamics of the reaction, thus indicating some positive feedback tendency. The shape of the growing C particle cluster is shown in Fig. 13. In contrast, for
the endothermic reaction shown in Fig. 14, in spite of the much lower reaction rate, the C particles remain in separated spots, uniformly distributed along the interfacial line. Fig. 15 shows the increasing number of C atoms with the development of the exothermic reactions. The results come from three 20 000-atom runs in temperatures T = 0.1, 0.3 and 0.7. As expected, the reaction proceeds much faster at the higher temperature. In T = 0.7 and 0.3 the speed of reaction monotonically decreases with time, as the contact line becomes filled with C atoms which separate A and B atoms. On the
W. Alda et al. / Physica D 146 (2000) 261–274
271
Fig. 13. Zoomed-in shot for the aspect ratio 5 run shown in Fig. 12. A and B particles are placed in the upper and lower halves of the box, respectively, while C particles are shown as darker dots.
Fig. 14. Global views of run with 100 000 initial particles for endothermic reactions. Dimensionless temperature is set to 0.3. The configuration has an aspect ratio of 5 and an area of 500 Å × 2500 Å.
272
W. Alda et al. / Physica D 146 (2000) 261–274
Fig. 15. Number of exothermic reactions taken globally as a function of time for various temperatures. Three temperatures have been considered, T = 0.1, 0.3 and 0.7.
other hand, at the lowest temperature, T = 0.1, the reaction speed is slow initially, then increases, and after roughly 0.4 ns decreases. This phenomenon is caused by the fact, that at low temperature only a few particles can react, only after a relatively long incubation time. Once a reaction takes place, a local increase of the temperature enhances the rate of the reactions and acts as a source of positive feedback. We do not observe this effect with endothermic reactions. The number of reactions with time (Fig. 16) shows that the speed of reaction decreases monotonically for all temperatures. Fig. 17 displays the temperature dependence of
Fig. 17. Number of reactions as a function of temperature for approximately a zero heat of reaction. Curves indicate subsequent instants.
number of reactions obtained in 10 runs with zero heat of the reaction. Runs have been performed with the same constant energy, as in the previous cases. Due to the potential chosen, the temperature remained roughly constant during the entire run, fluctuating only, with no drift. Two regions consisting of faster and slower growth of number of reactions varying with the temperature can be observed. There exists a threshold temperature dividing the two regions. It is roughly equal to 0.25 for short simulation time and moves down to T = 0.2 for t = 1.0 ns. This phenomenon can be explained by the scarcity of reactions at lower temperature due to the absence of any heat liberated by the reactions. Thus, the interfacial growth involving the C particles does not decelerate the global reaction rates.
5. Concluding remarks
Fig. 16. Number of endothermic reactions taken globally as a function of time for several temperatures. Five temperatures have been considered, T = 0.2, 0.4, 0.6, 0.8 and 1.0.
In this work, we have introduced for the first time a simple model for simulating chemical reactions, based on the formalism of classical MD. This technique has great potential of being utilized for investigating the collective nonlinear interaction involving millions to billions of reacting species being coupled to a larger-scale flow, thus allowing the possibilities for upscaling of chemical kinetics. Some of our prelim-
W. Alda et al. / Physica D 146 (2000) 261–274
inary results already offer some interesting insights, especially those concerning the enhancing effects of increasing the temperature, the marked heterogeneities produced at the reacting interface, the threshold nature of increasing the temperature and the fundamental issue regarding the concepts of local versus global rates of reaction. Our approach differs fundamentally from previous studies using discrete particles, by means of cellular automata or lattice-gas techniques, because we do not employ rules for updating the lattice dynamics but only two-body, LJ-like potentials with some chemical–physical meaning. This method is also different from other MD approaches, employing already assumed kinetic-rate equations [23]. Previous works in this area of MD with reactions have been concerned with equilibrium systems of rare gas solutions consisting of hard-spheres [25] and with a single LJ potential but relatively few particles, O(100), [26]. We are aware of the simplifications and drawbacks in this model, but our primary interest lies in building a computationally efficient scheme to study the collective interaction of millions to billions of reacting particles in the hydrodynamic mode and to understand some basic chemical–physical principles in many-body reactions. First, we do not account for quantum-mechanical effects, which would slow down considerably the computations and limit the number of participating atoms to at most several hundred. Second, this model lacks a three-body interaction term [11] for ensuring the formation of chemical bonds between the correct number of atoms. Third, we have not allowed for the possibilities of back reactions, involving C + A and C + B. Back reactions would certainly have a substantial nonlinear feedback effects and studies are currently under way to address this important issue. Lastly, there is the issue of longer timescales, necessary for achieving both hydrodynamic effects and more realistic reaction time rates. Extending these calculations to the temporal domain of microseconds would definitely require the use of massively parallel processors with large shared memory architecture of 10 GB or more. By elongating the simulation time, it is possible to observe the reactions along with their interactions with mixing
273
processes induced by larger-scale external flow fields, pressure or velocities applied to particles, or by the reaction itself as well. In this case, we can expect collective hydrodynamic-type motion of particles with a wide-range of length-scales. In the very near future, we will complete simulations of up to 107 reacting particles. The initial conditions we have employed is only an end-member. Other initial conditions involving, e.g., a prepared state of well-mixed reactants would also be of considerable interest in studies of chemical combustion. This method may be extended to DPD for studying at the mesoscale the interaction between mixing and chemical reactions in the novel cross-scale sense [24,27]. Many of the analytical tools developed for mixing in the macroscale (e.g. [28]), can be brought to bear on these newer cross-scale issues in the micro- to mesoscale. By sacrificing some of the details of fine-scale structures found in quantum-chemical reactions, we may now be able to connect, with this approach, chemical reacting flows across considerable spatial–temporal scales, which will be of direct relevance in the environmental and engineering disciplines [29,30], as well as in the chaos control of chemical reaction systems [31].
Acknowledgements We thank Witold Dzwinel, Rashid Zia and Julio Ottino for stimulating discussions. This research has been supported by the Geosciences Program of the Department of Energy (DOE). References [1] K. Binder, M. Müller, F. Schmid, Comput. Sci. Eng. 3 (1999) 10. [2] J. Mai, W. von Niessen, J. Chem. Phys. 98 (1993) 2032. [3] B. Chopard, M. Drozd, Europhys. Lett. (Switzerland) 15 (1991) 459. [4] J.R. Weimar, J.P. Boon, Phys. Rev. E 49 (1994) 1749. [5] J.P. Boon, D. Dab, R. Kapral, A. Lawniczak, Phys. Rep. 273 (1996) 556. [6] J.T. Wells, D.R. Janecky, B.J. Travis, Physica D 47 (1991) 115. [7] D. Rothman, S. Zaléski, Rev. Mod. Phys. 66 (1994) 1417. [8] M. Li, W.L. Johnson, W.A. Goddard III, Phys. Rev. B 54 (1996) 12067.
274
W. Alda et al. / Physica D 146 (2000) 261–274
[9] L.L. Lee, Y.S. Li, K.R. Wilson, J. Chem. Phys. 95 (1991) 2458. [10] J.P. Bergsma, J.T. Hynes, J.R. Reimers, K.R. Wilson, J. Chem. Phys. 85 (1986) 5625. [11] F.H. Stillinger, Th.A. Weber, J. Chem. Phys. 88 (1988) 5123. [12] P. Ortoleva, S. Yip, J. Chem. Phys. 65 (1976) 2045. [13] I. Stich, R. Car, M. Parinello, Phys. Rev. B 44 (1991) 4262. [14] J.R. Chelikowsky, N. Troullier, Y. Saad, Phys. Rev. Lett. 72 (1994) 1240. [15] K.C. Hass, W.F. Schneider, A. Curioni, W. Andreoni, Science 282 (1998) 265. [16] J.P. Boon, S. Yip, Molecular Hydrodynamics, McGraw-Hill, New York, 1980. [17] J. Mo´sci´nski, W. Alda, M. Bubak, W. Dzwinel, J. Kitowski, M. Pogoda, D.A. Yuen, in: D. Stauffer (Ed.), Annual Review of Computational Physics, Vol. V, World Scientific, Singapore, 1997, p. 97. [18] W. Alda, W. Dzwinel, J. Kitowski, J. Mo´sci´nski, M. Pogoda, D.A. Yuen, Comput. Phys. 12 (1998) 595.
[19] W. Dzwinel, W. Alda, M. Pogoda, D.A. Yuen, Physica D 137 (2000) 157. [20] P.J. Hoogerbrugge, J.M.V.A. Koelman, Europhys. Lett. 19 (1992) 155. [21] P. Español, Phys. Rev. E 57 (1998) 2930. [22] W. Dzwinel, D.A. Yuen, Mol. Sim. 22 (1999) 369. [23] K. Geisshirt, E. Praestgaard, S. Toxvaerd, J. Chem. Phys. 107 (1997) 9406. [24] W. Dzwinel, W. Alda, D.A. Yuen, Mol. Sim. 22 (1999) 397. [25] J. Gorecki, J. Gryko, J. Stat. Phys. 48 (1987) 329. [26] J. Gorecki, J. Gryko, Comput. Phys. Commun. 54 (1989) 245. [27] W. Dzwinel, D.A. Yuen, J. Colloid Interf. Sci. 225 (2000) 179. [28] J.M. Ottino, Scient. Am. 260 (1989) 40. [29] A.A. Ten, Y.Y. Podladchikov, D.A. Yuen, T.B. Larsen, A.V. Malevsky, Geophys. Res. Lett. 25 (1998) 3205. [30] J.M. Ottino, A. Souvaliotis, G. Metcalfe, Chaos, Solitons and Fractals 6 (1995) 425. [31] A.A. Hoff, H.H. Diebner, G. Baier, Z. Naturforsch 50a (1995) 1141.
Exothermic and endothermic chemical reactions involving very many particles modeled with molecular dynamics Witold Alda a , David A. Yuen b,∗ , Hans-Peter Lüthi c , James R. Rustad d a
Institute of Computer Science, University of Mining and Metallurgy AGH, 30-059 Krakow, Poland b Department of Geology and Geophysics and Minnesota Supercomputing Institute, University of Minnesota, 1200 Washington Ave. South, Minneapolis, MN 55415-1227, USA c Swiss Center for Scientific Computing, ETH, Zurich CH-8092, Switzerland d The Pacific Northwest National Lab, W.R. Wiley EMSL, Richland, WA 99352, USA Received 5 August 1999; accepted 3 July 2000 Communicated by R.P. Behringer
Abstract The traditional continuum approach of modeling chemical reactions with specified kinetic rates suffers from numerical difficulties in reactive flows and other highly non-equilibrium situations due to the stiffness of the differential equations in both space and time. These drawbacks can be eliminated within the framework of the discrete-particle approach in which the chemical reactions are modeled by means of two-body interparticle potentials with classical molecular dynamics. Here we present a simple model for prescribing the binary endothermic and exothermic reactions of the type A + B → C, in the presence of many (greater than 104 ) reacting particles. Two-body Lennard–Jones based potentials, Vij , have been utilized, in which the reaction takes place via an attractive potential VAB with a deep enough potential well. The other potentials VAC and VBC have deeper wells than VAB for creating a suitable situation for inducing exothermic reactions, while a high-energy plateau for VAC and VBC is used for simulating endothermic reactions. We have deployed 20 000 to 100 000 particles in a computational domain with an area of around 1 million Å2 . These ensembles have all been integrated out to around the order of nanoseconds. We have examined the bimolecular reactions between two layers initially consisting of A particles lying atop B particles, equally divided in the total population. Reactions take place first along this boundary. There is a nonlinear threshold phenomenon in the exothermic reactions, which are greatly accelerated by the heat liberated from the neighboring reactions. In contrast, the endothermic reactions are more subdued, as there is no positive feedback. Exothermic reactions occur in a spatially heterogeneous manner, especially at the lower temperatures, in which no single reaction rate can be assigned to the entire reaction front. More coherent spatial structures are found with decreasing temperature. These spatial structures can influence the overall flow field associated with the motions of the reacting particles. Timescales of the reactions are greatly lengthened by a reduction in the temperature in a threshold manner. © 2000 Published by Elsevier Science B.V. Keywords: Exothermic and endothermic reactions; Molecular dynamics; Lennard–Jones potential
1. Introduction
∗ Corresponding author. Tel.: +1-612-625-1818; fax: +1-612-624-8861 E-mail address: [email protected] (D.A. Yuen).
Today chemical reacting flows form indeed important and fundamental foci for many scientific and engineering endeavors. Numerous examples can be drawn from fields as diverse as metamorphic reacting fronts
0167-2789/00/$ – see front matter © 2000 Published by Elsevier Science B.V. PII: S 0 1 6 7 - 2 7 8 9 ( 0 0 ) 0 0 1 5 0 - 0
262
W. Alda et al. / Physica D 146 (2000) 261–274
imprinted in Alpine rocks and the ionization processes accompanying the reentry of hypersonic vehicles descending back into the earth’s atmosphere. Most of the emphasis on the study of chemically reacting flows has been directed toward the continuum approach based on solving nonlinear partial differential equations with the chemical reactions added in as a set of nonlinear ordinary differential equations in time. These chemical conservation relations are generally prescribed by given kinetic-rate laws, which are thermally activated, and hence strongly temperature-dependent in the Arrhenius sense. However, solving these coupled species-energy-momentum equations is beset by the multiple-scale nature of the coupled equations, which baffle even the best of numerical methods, especially at the chemically reacting fronts with sharp gradients. The stiffness of the partial differential equations arising from the use of both the spatial and temporal operators in the continuum approach can be ameliorated by using the discrete methods, such as cellular automata, lattice-gas, lattice-Boltzmann and molecular dynamics (MD), as they operate over much shorter timescales, many orders in magnitude smaller than the time steps taken in the continuum approach. The above methods are based on very simplified models and often they do not give accurate or even realistic results, but according to Binder et al. [1], “. . . the goal of such simulations is not to predict numbers but to provide better understanding”. In our case, this means a better understanding of the dynamics of chemical reaction processes with many particles in highly non-equilibrium situations. Cellular-automaton methods have been used to study diffusion and reaction processes by Mai and von Niessen [2], Chopard and Drozd [3] as well as Weimar and Boon [4]. There have been some previous modeling of chemical reactions at mineral surfaces by lattice gas automata (LGA) model by Boon et al. [5], Wells et al. [6] as well as Rothman and Zaléski [7]. They found that LGA calculations appeared to be well suited for problems in the mesoscopic scale as in the range of pore-fluid domain, of the order of tens of microns. Classical MD along with simple two-body potentials has been used to study phase transitions in
two-dimensional Lennard–Jones (LJ) binary systems by Li et al. [8]. The efficacy of this simplified approach for long-time integration, involving 1000 to 10 000 particles for over 107 time steps, was emphasized because of the need to monitor the onset of melting, a relatively long timescale phenomenon. Chemical reaction simulations on the atomic level have been performed by Wilson and coworkers using classical MD (e.g. [9]) or transition state theory [10]. Stillinger and Weber [11] have used two-body and three-body potentials in MD simulations to observe spontaneous forming of diatomic molecules. Ortoleva and Yip [12] used a hard disks system for binary reactions, where particles react by collisions (forward reactions) and by spontaneous transitions (backward reactions), giving different steady states depending on the reaction rates. In contrast to our approach, all these simulations were performed in equilibrium systems with exchange reactions only and constant number of molecules. On the other hand, direct integration of the coupled quantum-mechanical self-consistent molecular dynamics (SCMD) problem [13–15] for several hundred atoms and for timescales, up to the nanosecond range, remains somewhat out of reach with the present generation of computers. Thus, the time is now for some middle-ground solution for intermediate timescales by means of MD, which can address this problem of bimolecular reactions with a large swarm of heterogeneous atoms which can display their collective fluid-like nature [16]. Some macroscopic fluid-dynamical phenomenon can be manifested, as in the MD study of the fluid-dynamical instabilities by Mo´sci´nski et al. [17], Alda et al. [18] and Dzwinel et al. [19], in which as many as 107 particles were integrated on massively parallel computers for around 106 time steps, in order to simulate mesoscale fluid-dynamical instabilities, such as the Rayleigh–Taylor instability. In this same spirit, we present in this paper a simpler conceptual and computational model for investigating binary chemical reactions, based on classical MD and two-body interparticle potentials. With the MD simulations, we can systematically investigate the essential characteristics of chemically reacting hydrodynamics with far more participating atoms, exceeding by
W. Alda et al. / Physica D 146 (2000) 261–274
several orders in magnitude what would be feasible by the SCMD technique, and for much longer timescales. The physical scale in which we will be operating still falls within the microscale range of under 1 m but the spatial range can be extended to the mesoscale domain of Wells et al. [6] by combining the dissipative particle dynamics (DPD) method [20,21], where particles represent clusters of atoms, with classical atomic MD particles. In the above approach, one can think of a system built up of large DPD particles representing non-reacting solvent with dissolved reacting MD atoms. The method described here will potentially constitute a powerful new tool for the computational treatment of chemical reacting hydrodynamics in the microscale and also eventually be connected to the mesoscale of the order of microns by DPD methods [22]. Our approach stands in sharp contrast to MD studies in which the chemical kinetics are assumed a priori for the reacting particles [23].
2. MD model Our simulations have been carried out in two dimensions. We have employed a two-body LJ-based interparticle potential, for both exothermic and endothermic reactions, modified for certain interactions: " # σαβ 6 σαβ 12 LJ − , (1) Vαβ (rij ) = 4αβ rij rij where α and β denote the particle species, αβ the well depth of the potential, σαβ a parameter determining the nearest neighbor distance in the interaction, and rij the distance between the particles i and j . All of the potentials are truncated at certain rcut and shifted up, so as to reach zero potential at cut-off and achieve potential continuity at the cut-off radius rcut . In this model, three different types of particles, A, B and C, which may differ in mass and interaction potential parameters, are employed. All interactions are assumed to have identical σαβ = σ parameters; thus all particles have the same size. The total number of particles in the system changes, as the reactions develop. The VAB potential for both types of reactions has been modified in this study by including a hump
263
(Figs. 1a and 2a) in the overall potential, which is given by LJ VAB (rij ) = VAB + Vhump .
(2)
The hump is modeled by piecing together two cubic LJ (r ), where r = ( 26 )1/6 σ , splines, which start at VAB sh sh 7 go through the top of the hump at rh and end at LJ (r ), where r = 2r − r . V VAB AB is truncated at eh eh h sh rcut = 2.5σ ; thus the hump lies within the interaction range. The idea of using this artificial feature is to enable the reactions to take place only for particles which have a sufficiently high kinetic energy. The VAA , VBB and VCC potentials, which are also identical for both the exothermic and endothermic reactions (Figs. 1c and 2c), are truncated at rcut = LJ , and shifted up by , 1.12σαβ , a minimum of Vαβ αβ thus making the potential purely repulsive in this study. The VBC and VAC potentials differ for both types of reactions. For exothermic reactions, the potential well depth, , in VBC and VAC is deeper than in VAB (Fig. 1b). This results in smaller average potential energy of A–C and B–C interactions than A–B that occurred before the reaction. Assuming that the total energy of the system is kept constant, the average kinetic energy increases. Changing the shape of the VAC and VBC functions by moving them up, increases the average A–C and B–C potential energy, eventually allowing the endothermic reactions to occur. For this purpose, we have chosen a very simple criterion, by introducing the repulsive potential with a flat plateau (Fig. 2b), where particles in A–C and B–C interactions can lose their kinetic energy. In this way it is possible to find the A–C and B–C potential for which there will be zero heat generated by the reaction. With A–B and A–A, B–B, C–C potentials used, this can be achieved with a very shallow potential well depth, = 20 J/K, in the A–C and B–C interactions. In addition to the exothermic and endothermic runs, we have carried out a series of zero heat runs at different temperatures in order to demonstrate the temperature dependence of the reaction speed. These results are shown at the end of this paper. The parameters of the exothermic and endothermic potentials are listed in Table 1.
264
W. Alda et al. / Physica D 146 (2000) 261–274
Fig. 1. The form of the three different kinds of energy potentials used in producing the binary exothermic reaction: (a) VAB represents the potential which produces the reaction A + B → C. Reaction takes place when both A and B remain together within one time step inside the well. The hump in front of the well controls the rate of the reaction; (b) VCB and VCA have deeper potential than VAB to enable exothermic nature of the reaction; (c) VAA , VBB , VCC are repulsive in nature. Thus, there are no clusters involving the formed A, B and C particles.
Fig. 2. The form of the three different kinds of potentials used in creating the binary endothermic reaction: (a) VAB represents the potential which produces the reaction A + B → C, identical to that in Fig. 1; (b) VCB and VCA have the plateau on which the endothermic reaction takes place. (c) VAA , VBB , VCC are repulsive in nature, identical to those in Fig. 1.
W. Alda et al. / Physica D 146 (2000) 261–274
265
Table 1 Potential parameters used in the three-component system with exothermic and endothermic reactions Interaction
αβ (J/kB )
σαβ (Å)
rcut (σαβ /Å)
Plateau cut-off (σαβ /Å)
Potentials equal for exothermic and endothermic reactions A–B 119.8 3.4 2.5 – A–A 119.8 3.4 1.12 – B–B 119.8 3.4 1.12 – C–C 119.8 3.4 1.12 – Exothermic only C–B 239.6 3.4 2.5 – C–A 239.6 3.4 2.5 – Endothermic only C–B 119.6 3.4 1.12 0.9 C–A 119.6 3.4 1.12 0.9
The mass of A and B particles has been set to 40 amu, while the mass of C particle, a product of A + B reaction, to 80 amu. The particles are confined within a rectangular box, and initially are placed on a grid forming a mesh built up from 2D regular squares. With the value of σ assumed, the size of such a cell is about 3.53 Å × 3.53 Å, which yields a total area of the computational box greater than 106 Å2 (Fig. 3) for a simulation with 105 particles. The initial configuration is not close packed, so the particles can diffuse more easily. This involves the initial stage of the
Hump height (J/kB )
Hump distance (σαβ /Å)
20.0 – – –
1.6 – – –
– –
– –
– –
– –
simulation, as the particles can start reacting earlier. In the course of simulation the structure is not kept and the reacting particles can form a collective fluid-like particle flow. The size of the box is kept constant. Periodic boundary conditions have been applied along the vertical edges, while for the horizontal boundary we have introduced reflecting walls. The time step, 1t, has been set to 10−14 s, which puts an upper-bound on the relative velocity of the interacting particles. With this model, we can achieve a speed of around two time steps per second for 105
Fig. 3. Schematic diagram of the computational domain, where the exothermic reactions take place (see colliding particles depicted). The initial grid layout with the regular square configuration of particles A and B is displayed.
266
W. Alda et al. / Physica D 146 (2000) 261–274
particles on a 600 MHz workstation. We can then complete 106 time step run (or 10 nanoseconds with the given 1t) within 1 day.
3. Algorithms for inducing binary endothermic and exothermic reactions For simulating chemical reactions with both endothermic and exothermic nature in MD, we have implemented simple algorithms to a model of A + B → C chemical reactions to a classical MD scheme based on LJ particles [24]. By a reaction, we mean here an “instantaneous” or an extremely short time-instant involving the transformation of particles A and B into a single particle C, when certain conditions concerning the duration of time spent inside the well of VAB and velocity magnitudes of particles are fulfilled. The occurrence and reaction speed depend on the four criteria: 1. the potential depths and widths applied, 2. the average energies or speed of the ensemble of particles, 3. the global thermodynamic conditions, such as the density and the temperature, and 4. the initial boundary conditions. In order for a reaction to occur between particles A and B, the VAB potential must have a potential well deep enough. The “hump” added in VAB can be especially important in low-density systems, where mean rAB goes beyond the hump. In such a case, the hump acts as a threshold. The non-reacting A–A, B–B and C–C pairs, interact by repulsive potential only (Fig. 1b). Repulsion prevents the forming of clusters of same-kind particles and enables better mixing. Particles A–C and B–C attract each other; however, the potential well depth is deeper than in the A–B interaction. Moreover, we do not allow C to react with A and B in this study, i.e. no back reaction. In this model, the principal criterion for the occurrence of a chemical reaction, is the interparticle distance. We assume that a reaction takes place for the condition rAB < rcrit , where rcrit is a certain critical distance. We also assume that the rAB < rcrit state should be kept for a certain number of time steps, con-
tinuously. The reaction rate can be slowed down then, by decreasing rcrit or increasing the time interval. In the simulations performed, we have assumed that the radius for reactive occurrence, rAB < rcrit , in a single time step, O(10−14 ) s, is sufficient for the reaction. In general, rcrit is chosen to lie within 0.8 . . . 0.9 of σ in A–B interaction. To avoid very long runs with very few reactions, we have chosen a relatively large value of rcrit = 0.87 to speed up the reactions. The model reactions do fulfill the conservation laws. Conservation of mass and momentum is straightforward, as mC = mA + mB and pEC = pEA + pEB . The center of mass of C is equal to the center of mass of A + B. The total energy, however, is not automatically conserved because each time the reaction occurs, two particles A and B are replaced with one particle C, which interacts with its nearest neighbors by a different potential. The particle C has no internal dynamical structure and cannot keep any hidden potential energy. To keep the total energy constant, we can only change the kinetic energy of some particles whenever a reaction takes place. We accomplish this task at each time step by adjusting the kinetic energy of all particles within 2.5σαβ radius around C, so as to conserve the total energy of the neighboring area at each reaction. 4. Results We have used a constant energy scaling in conducting the simulations of a system involving 20 000 particles for both exothermic and endothermic reactions for a range of temperatures in order to see some noticeable effects. The changes of potential, kinetic and total energies with time are displayed in Figs. 4 and 5. Exothermic runs have been performed in three initial dimensionless temperatures, T = 0.1, 0.3 and 0.7 in reduced LJ units. As shown in Fig. 4, for T = 0.1 and 0.7 runs, the total energy is conserved, while the potential and kinetic energies vary with time. The kinetic energy increases because of the exothermic nature of the reaction. Due to the increase in the number of particles C and mixing at the interface area, the potential energy now has a tendency to decrease, while the temperature of the entire system increases. A similar, but reversed situation is shown in Fig. 5,
W. Alda et al. / Physica D 146 (2000) 261–274
Fig. 4. Conservation of the total energy with time for exothermic reactions. Dimensionless initial temperature is 0.1–0.5 by the end of the run. 20 000 particles consisting of A and B initially were employed. The total number of particles changes in time due to the reactions.
in which with endothermic reactions the potential energy increases, while the temperature decreases monotonically. Figs. 6–8 show the snapshots of the particle locations for the exothermic reactions in three different initial temperatures, T = 0.7, 0.3 and 0.1, and in three different instants of time. Gray areas represent atoms A and B in the lower and upper parts of the box, while
Fig. 5. Conservation of the total energy with time for endothermic reactions. Dimensionless initial temperature is 0.2 decreasing to 0.08 by the end of the run. 20 000 particles consisting of A and B initially were employed. The total number of particles varies with time due to the chemical reactions.
267
the atoms C, the product of the reaction, are shown as darker dots. The latter are displayed in larger format for better visualization. Initially the system consists only of atoms A and B occupying the upper and lower halves of the box, separated by a straight line of contact. It can be observed from the snapshots, that in higher temperature the reaction starts much faster and quickly spreads over the entire interface. At the highest temperature (T = 0.7) the particle motion is more chaotic and the reactions take place along the entire front. On the other hand, reactions at T = 0.1 are sporadic and induce a heterogeneous character to the interface. This is due to the exothermic character of the reactions and the fact that, once the reaction occurs, the next reaction becomes more probable to occur in the neighborhood. Timescales are faster at the higher temperatures. One can also see that the reaction front forms a curve line, which is rugged but its shape remains periodic with a dominant wavelength. The curvature of the reaction interface can be explained by the exothermic character of the reactions and the tendency of the heat liberated to stretch the interface boundary, which is energetically adjusted. On the other hand, the reactions rate is slowed down with time due to the growing thickness of C particles layer which separates the reacting A and B species. In these three runs, starting from T = 0.7, 0.3 and 0.1, the temperatures rise asymptotically to T = 0.98, 0.60 and 0.42, respectively. Similar snapshots showing the endothermic reactions with the initial temperatures, T = 0.6, 0.4 and 0.2, are shown in Figs. 9–11. We observe that the image of the initial phase of the reactions in Figs. 9 and 10 is similar to that in the exothermic model. Furthermore, the interface line remains rather straight, with only short-wavelength fluctuations. Particles C are more scattered (especially in Fig. 10) due to the smaller probability of the next reaction in the vicinity of an already reacted site. A larger domain has been established by elongating the computational domain by fivefold. Two simulations for both exothermic and endothermic reactions have been performed for the larger box in the initial temperature T = 0.3. We observe in Fig. 12, that the spatial heterogeneity of the reactions is much
268
W. Alda et al. / Physica D 146 (2000) 261–274
Fig. 6. Global views of the particles for exothermic reactions. Darker dots in the central part denote C, while A and B are placed in the upper and lower halves of the box, respectively. Dimensionless initial temperature is set to 0.7.
Fig. 7. Global snapshots of the exothermically reacting particles for dimensionless initial temperature of T = 0.3.
Fig. 8. Global snapshots of the exothermically reacting particles for dimensionless initial temperature of T = 0.1.
W. Alda et al. / Physica D 146 (2000) 261–274
269
Fig. 9. Global views of the particle distribution for endothermic reactions. Darker dots in the central part denote C, while A and B are placed in the upper and lower halves of the box, respectively. The dimensionless initial temperature is set to 0.6.
Fig. 10. Global snapshots of the endothermically reacting particles for dimensionless initial temperature of T = 0.4.
Fig. 11. Global snapshots of the endothermically reacting particles for a dimensionless initial temperature of T = 0.2.
270
W. Alda et al. / Physica D 146 (2000) 261–274
Fig. 12. Global views of run with 100 000 initial particles for exothermic reactions. Dimensionless temperature is set to 0.3. The configuration has an aspect ratio of 5 and an area of 500 Å × 2500 Å.
stronger than in the smaller boxes. This tendency is probably due to greater length of the reaction front and smaller influence of periodic boundary conditions. The occurrence of the reaction generates heat, which is then distributed over the neighboring particles. The energy distribution along the interface is slower, when not constrained by the end conditions. This also shows that the speed and the range of the kinetic energy being distributed about, can greatly influence the overall dynamics of the reaction, thus indicating some positive feedback tendency. The shape of the growing C particle cluster is shown in Fig. 13. In contrast, for
the endothermic reaction shown in Fig. 14, in spite of the much lower reaction rate, the C particles remain in separated spots, uniformly distributed along the interfacial line. Fig. 15 shows the increasing number of C atoms with the development of the exothermic reactions. The results come from three 20 000-atom runs in temperatures T = 0.1, 0.3 and 0.7. As expected, the reaction proceeds much faster at the higher temperature. In T = 0.7 and 0.3 the speed of reaction monotonically decreases with time, as the contact line becomes filled with C atoms which separate A and B atoms. On the
W. Alda et al. / Physica D 146 (2000) 261–274
271
Fig. 13. Zoomed-in shot for the aspect ratio 5 run shown in Fig. 12. A and B particles are placed in the upper and lower halves of the box, respectively, while C particles are shown as darker dots.
Fig. 14. Global views of run with 100 000 initial particles for endothermic reactions. Dimensionless temperature is set to 0.3. The configuration has an aspect ratio of 5 and an area of 500 Å × 2500 Å.
272
W. Alda et al. / Physica D 146 (2000) 261–274
Fig. 15. Number of exothermic reactions taken globally as a function of time for various temperatures. Three temperatures have been considered, T = 0.1, 0.3 and 0.7.
other hand, at the lowest temperature, T = 0.1, the reaction speed is slow initially, then increases, and after roughly 0.4 ns decreases. This phenomenon is caused by the fact, that at low temperature only a few particles can react, only after a relatively long incubation time. Once a reaction takes place, a local increase of the temperature enhances the rate of the reactions and acts as a source of positive feedback. We do not observe this effect with endothermic reactions. The number of reactions with time (Fig. 16) shows that the speed of reaction decreases monotonically for all temperatures. Fig. 17 displays the temperature dependence of
Fig. 17. Number of reactions as a function of temperature for approximately a zero heat of reaction. Curves indicate subsequent instants.
number of reactions obtained in 10 runs with zero heat of the reaction. Runs have been performed with the same constant energy, as in the previous cases. Due to the potential chosen, the temperature remained roughly constant during the entire run, fluctuating only, with no drift. Two regions consisting of faster and slower growth of number of reactions varying with the temperature can be observed. There exists a threshold temperature dividing the two regions. It is roughly equal to 0.25 for short simulation time and moves down to T = 0.2 for t = 1.0 ns. This phenomenon can be explained by the scarcity of reactions at lower temperature due to the absence of any heat liberated by the reactions. Thus, the interfacial growth involving the C particles does not decelerate the global reaction rates.
5. Concluding remarks
Fig. 16. Number of endothermic reactions taken globally as a function of time for several temperatures. Five temperatures have been considered, T = 0.2, 0.4, 0.6, 0.8 and 1.0.
In this work, we have introduced for the first time a simple model for simulating chemical reactions, based on the formalism of classical MD. This technique has great potential of being utilized for investigating the collective nonlinear interaction involving millions to billions of reacting species being coupled to a larger-scale flow, thus allowing the possibilities for upscaling of chemical kinetics. Some of our prelim-
W. Alda et al. / Physica D 146 (2000) 261–274
inary results already offer some interesting insights, especially those concerning the enhancing effects of increasing the temperature, the marked heterogeneities produced at the reacting interface, the threshold nature of increasing the temperature and the fundamental issue regarding the concepts of local versus global rates of reaction. Our approach differs fundamentally from previous studies using discrete particles, by means of cellular automata or lattice-gas techniques, because we do not employ rules for updating the lattice dynamics but only two-body, LJ-like potentials with some chemical–physical meaning. This method is also different from other MD approaches, employing already assumed kinetic-rate equations [23]. Previous works in this area of MD with reactions have been concerned with equilibrium systems of rare gas solutions consisting of hard-spheres [25] and with a single LJ potential but relatively few particles, O(100), [26]. We are aware of the simplifications and drawbacks in this model, but our primary interest lies in building a computationally efficient scheme to study the collective interaction of millions to billions of reacting particles in the hydrodynamic mode and to understand some basic chemical–physical principles in many-body reactions. First, we do not account for quantum-mechanical effects, which would slow down considerably the computations and limit the number of participating atoms to at most several hundred. Second, this model lacks a three-body interaction term [11] for ensuring the formation of chemical bonds between the correct number of atoms. Third, we have not allowed for the possibilities of back reactions, involving C + A and C + B. Back reactions would certainly have a substantial nonlinear feedback effects and studies are currently under way to address this important issue. Lastly, there is the issue of longer timescales, necessary for achieving both hydrodynamic effects and more realistic reaction time rates. Extending these calculations to the temporal domain of microseconds would definitely require the use of massively parallel processors with large shared memory architecture of 10 GB or more. By elongating the simulation time, it is possible to observe the reactions along with their interactions with mixing
273
processes induced by larger-scale external flow fields, pressure or velocities applied to particles, or by the reaction itself as well. In this case, we can expect collective hydrodynamic-type motion of particles with a wide-range of length-scales. In the very near future, we will complete simulations of up to 107 reacting particles. The initial conditions we have employed is only an end-member. Other initial conditions involving, e.g., a prepared state of well-mixed reactants would also be of considerable interest in studies of chemical combustion. This method may be extended to DPD for studying at the mesoscale the interaction between mixing and chemical reactions in the novel cross-scale sense [24,27]. Many of the analytical tools developed for mixing in the macroscale (e.g. [28]), can be brought to bear on these newer cross-scale issues in the micro- to mesoscale. By sacrificing some of the details of fine-scale structures found in quantum-chemical reactions, we may now be able to connect, with this approach, chemical reacting flows across considerable spatial–temporal scales, which will be of direct relevance in the environmental and engineering disciplines [29,30], as well as in the chaos control of chemical reaction systems [31].
Acknowledgements We thank Witold Dzwinel, Rashid Zia and Julio Ottino for stimulating discussions. This research has been supported by the Geosciences Program of the Department of Energy (DOE). References [1] K. Binder, M. Müller, F. Schmid, Comput. Sci. Eng. 3 (1999) 10. [2] J. Mai, W. von Niessen, J. Chem. Phys. 98 (1993) 2032. [3] B. Chopard, M. Drozd, Europhys. Lett. (Switzerland) 15 (1991) 459. [4] J.R. Weimar, J.P. Boon, Phys. Rev. E 49 (1994) 1749. [5] J.P. Boon, D. Dab, R. Kapral, A. Lawniczak, Phys. Rep. 273 (1996) 556. [6] J.T. Wells, D.R. Janecky, B.J. Travis, Physica D 47 (1991) 115. [7] D. Rothman, S. Zaléski, Rev. Mod. Phys. 66 (1994) 1417. [8] M. Li, W.L. Johnson, W.A. Goddard III, Phys. Rev. B 54 (1996) 12067.
274
W. Alda et al. / Physica D 146 (2000) 261–274
[9] L.L. Lee, Y.S. Li, K.R. Wilson, J. Chem. Phys. 95 (1991) 2458. [10] J.P. Bergsma, J.T. Hynes, J.R. Reimers, K.R. Wilson, J. Chem. Phys. 85 (1986) 5625. [11] F.H. Stillinger, Th.A. Weber, J. Chem. Phys. 88 (1988) 5123. [12] P. Ortoleva, S. Yip, J. Chem. Phys. 65 (1976) 2045. [13] I. Stich, R. Car, M. Parinello, Phys. Rev. B 44 (1991) 4262. [14] J.R. Chelikowsky, N. Troullier, Y. Saad, Phys. Rev. Lett. 72 (1994) 1240. [15] K.C. Hass, W.F. Schneider, A. Curioni, W. Andreoni, Science 282 (1998) 265. [16] J.P. Boon, S. Yip, Molecular Hydrodynamics, McGraw-Hill, New York, 1980. [17] J. Mo´sci´nski, W. Alda, M. Bubak, W. Dzwinel, J. Kitowski, M. Pogoda, D.A. Yuen, in: D. Stauffer (Ed.), Annual Review of Computational Physics, Vol. V, World Scientific, Singapore, 1997, p. 97. [18] W. Alda, W. Dzwinel, J. Kitowski, J. Mo´sci´nski, M. Pogoda, D.A. Yuen, Comput. Phys. 12 (1998) 595.
[19] W. Dzwinel, W. Alda, M. Pogoda, D.A. Yuen, Physica D 137 (2000) 157. [20] P.J. Hoogerbrugge, J.M.V.A. Koelman, Europhys. Lett. 19 (1992) 155. [21] P. Español, Phys. Rev. E 57 (1998) 2930. [22] W. Dzwinel, D.A. Yuen, Mol. Sim. 22 (1999) 369. [23] K. Geisshirt, E. Praestgaard, S. Toxvaerd, J. Chem. Phys. 107 (1997) 9406. [24] W. Dzwinel, W. Alda, D.A. Yuen, Mol. Sim. 22 (1999) 397. [25] J. Gorecki, J. Gryko, J. Stat. Phys. 48 (1987) 329. [26] J. Gorecki, J. Gryko, Comput. Phys. Commun. 54 (1989) 245. [27] W. Dzwinel, D.A. Yuen, J. Colloid Interf. Sci. 225 (2000) 179. [28] J.M. Ottino, Scient. Am. 260 (1989) 40. [29] A.A. Ten, Y.Y. Podladchikov, D.A. Yuen, T.B. Larsen, A.V. Malevsky, Geophys. Res. Lett. 25 (1998) 3205. [30] J.M. Ottino, A. Souvaliotis, G. Metcalfe, Chaos, Solitons and Fractals 6 (1995) 425. [31] A.A. Hoff, H.H. Diebner, G. Baier, Z. Naturforsch 50a (1995) 1141.