Concentration fluctuation caused by reaction-diffusion coupling near catalytic active sites

Haolei Zhang,Mingcan Zhao,Yanping Li,Chengxiang Li,3,*,Wei Ge,4,*

1 State Key Laboratory of Multiphase Complex Systems,Institute of Process Engineering,Chinese Academy of Sciences,Beijing 100190,China

2 School of Chemical Engineering,University of the Chinese Academy of Sciences,Beijing 100049,China

3 Dalian National Laboratory for Clean Energy,Chinese Academy of Sciences,Dalian,Liaoning 116023,China

4 Innovation Academy for Green Manufacture,Chinese Academy of Sciences,Beijing 100049,China

Keywords:Molecular dynamics simulation Concentration fluctuation Reaction-diffusion coupling Catalyst

ABSTRACT The coupling of reaction and diffusion between neighboring active sites in the catalyst pore leads to the spatiotemporal fluctuation in component concentration,which is very important to catalyst performance and hence its optimal design.Molecular dynamics simulation with hard-sphere and pseudo-particle modeling has previously revealed the non-stochastic concentration fluctuation of the reactant/product near isolated active site due to such coupling,using a simple model reaction of A →B in 2D pores.The topic is further developed in this work by studying the concentration fluctuation due to such coupling between neighboring active sites in 3D pores.Two 3D pore models containing an isolated active site and two adjacent active sites were constructed,respectively.For the isolated site,the concentration fluctuation intensifies for larger pores,but the product yield decreases,and for a given pore size,the product yield reaches a peak at a certain reactant concentration.For two neighboring sites,their distance (d) is found to have little effect on the reaction,but significant to the diffusion.For the same reaction competing at both sites,larger d leads to more efficient diffusion and better overall performance.However,for sequential reactions at the two sites,higher overall performance presents at a smaller d.The results should be helpful to the catalyst design and reaction control in the relevant processes.

1.Introduction

Heterogeneous catalyst pellets often have complex internal catalytic surface,active sites distribution and pore structures.With the coupling of reaction and diffusion,they can significantly affect the catalytic performance[1]in terms of selectivity[2],coke deposition [3-6],and so on.Understanding their coupling mechanism and then designing and modulating the reaction-diffusion coupling in these confined structures [7,8] are,therefore,important to optimal catalytic performance.

Many studies on the diffusion [9-11] or reaction [12] process alone at the molecular level have been reported widely in literature,and some have addressed their coupling,which roughly belongs to two categories:One for coupling with intrinsic competing mechanisms in reaction kinetics due to changes in reaction components or reaction sites,which generally leads to periodic oscillations in chemical reactions and can be observed by in-situ electron microscopy [13-17],for example,the chemical and morphological change of catalyst surface[16]is attributed to the competing of simultaneously presenting oxidizing and reducing agents.Oscillating spatio-temporal patterns formed by reaction fronts spreading during the ongoing H2oxidation on a Rh crystal have also been observed in experiments[17]as another example in this category.

The other category of reaction-diffusion coupling involves the steric effect,which leads to complex chemical dynamics as the reactants diffuse to the active sites in the pores and then the products diffuse out.For example,the competition between the diffusion of reactants into the catalyst pores and the outward diffusion of products leads to the decrease of the overall performance of the catalyst [18].It is generally considered that such steric effect is caused by confined structures and boundary factors such as micro-pores [19].In fact,even without the influence of confined space,the interaction between many reactant/product molecules may result in steric effect also,and the reaction-diffusion coupling will lead to complex kinetic process.However,such research is very few.In our earlier simulation work [20],the local concentration fluctuation due to steric hindrance effect between molecules,studied using a simple model reaction in a 2D pore structure containing an isolated active site,is found to be nonstochastic and its magnitude depends on the coupling of diffusion and reaction.It indicates that only static pictures of the active sites are insufficient to describe the concentration distribution inside the catalyst pellet,catalytic dynamics and mechanism.

In an actual reaction system,the reaction-diffusion coupling between neighboring active sites due to steric effect will be more complex.Recently,reactions at neighboring isolated sites were found to be interrelated and affect the formation of a synergetic reaction network [21].And for bifunctional catalysts,the distance between active sites were found to affect the products selectivity and even the reaction mechanism [22-25].An in-depth understanding of the reaction-diffusion coupling is important for establishing theoretic model to study catalytic reaction behavior.However,it is difficult to observe the dynamics process,such as the concentration fluctuation,at atomic-scale in experiments.Molecular dynamics (MD) simulation is very suitable for studying such microscopic details because it can provide information on the spontaneous temporal and spatial variation of concentration [26-30].Therefore,a 3D pore model with reaction resembling 2-butene isomerization into 1-butene over mesoporous silica catalyst was constructed based on our previous work [20].Then the reaction-diffusion coupling in two pore structures,containing an isolated active site and two adjacent active sites respectively,were investigated.The characteristics of concentration fluctuation and its relationship with the distance between the two adjacent sites were analyzed.

2.Computational Model and Method

2.1.Molecular simulation method

Hard-sphere/pseudo-particle modeling(HS-PPM)is a simplified MD approach with much higher efficiency while providing reasonable results,especially for dilute gas [26,31,32].The molecules are constructed by hard spheres undergoing free flights and instantaneous elastic collisions with their velocities updated by.

where,m,P,Vare mass,position and velocity of the colliding spheres denoted by subscripts 1 and 2,and the superscript (′)means ‘‘after collision”.

The interaction between two hard spheres is described by step potential.This potential only needs one parameter,the effective diameter σHS,which can be calculated according to the parameters σLJand εLJfrom the most commonly used Lennard-Jones (L-J)potential [33]:

whereT*=kBT/εLJ,andkBandTare the Boltzmann constant and the absolute temperature,respectively.In this work,the parameters of L-J potential for two butenes are σLJ=7.142 × 10-10m and εLJ/kB=209 K,respectively [34].

2.2.Reaction model

A simplified reaction model was established based on the isomerization of 2-butene to 1-butene over mesoporous silica catalyst.The 2-butene is adsorbed on the active site,then further converted to 1-butene through positional isomerization[35].Here,the formation of inert or immobile molecules are not considered.Though the typical gas-solid reaction includes a series of processes,such as diffusion,adsorption,reaction and desorption,for simplicity,several surface processes are treated as a lumped reaction process.Such simplification is often used as a model in infrared radiation (IR) microscopy[36] and Fourier transform infrared (FTIR) spectroscopy [37]microimaging experiments.In this way,the chemical reaction was simplified as a first order irreversible isothermal reaction process from reactant A (2-butene) to product B (1-butene)(A →B),then only the lumped reaction and diffusion processes were considered in our model.

In the reaction model,the physical parameters of molecules,such as massmand effective diameter σHSremain unchanged during the reaction,but a flag is used to distinguish the reactant (A) and product (B).Whether the reaction occurs or not is determined based on the chemical reaction collision theory.When a reactant molecule collides with the active site,the reaction will occur only when its translational energy corresponding to the component of their relative velocity along the line of their centers reaches a threshold value,/2.Here,the critical velocity vccan be expressed in terms of the thermal velocity v0which depends on the temperature (T).That is,when a collision occurs in simulation,a reactant molecule (A)will be labeled as product molecule (B) if its velocity vtsatisfieswherep0is the reaction coefficient.Obviously,the critical energy is related to the activation energy [38] and can be altered by changingp0.

2.3.Characterization of diffusion,reaction performance and fluctuation

As defined in our previous work [20],the diffusion efficiency(Deff),or the fraction of the molecules colliding with the active site to all those entering the pore;the reaction efficiency (Reff),or the fraction of effective reactional collisions to all collisions on the active site;and the yield (Y),or the number of product molecules leaving the pore in unit time,are used to characterize the whole process quantitatively.

To quantitatively characterize the local concentration fluctuation,a small region near the active site,Region S shown in Fig.1,was used for the data statistics.The number of reactants and products entering (or leaving) this region in unit time were counted and denoted asNrandNp,respectively.At the same time,the total number of molecules was denoted asNa=Nr+Np.Thus,Nr,NpandNacan be regarded as stochastic variables.After the system reaches steady state,taking reactant as an example,the standard deviation (σr) ofNrcan be calculated as:

wheretnis the number of time steps fromt0totandis the timeaveraged value ofNrin this time interval.Then,the coefficient of variation (μr) was used to represent the magnitude of fluctuation ofNrand calculated by:

Fig.1.Schematic diagram of the 3D pore model (the molecules are shown in larger size to improve the overall visibility).

In addition,the completely stochastic concentration fluctuation in a similar system should be quantified for comparison.As explained in our previous work[20],it corresponds to the concentration fluctuation of total molecules (Na) without reaction,which is characterized by the square root of the average number of molecules.So,a parameter φr?awas defined to compare the relative magnitudes of the fluctuations and obtained directly from μr:

Obviously,the fluctuation is non-stochastic if φr?ais larger than 1.0 and it can also quantify the relative magnitudes of fluctuation obtained from different simulation cases.

2.4.Simulation setup

The 3D pore model is shown in Fig.1,the parameters for mass,length and time were nondimensionalized by the massm(9.32 × 10-26kg),effective diameter σHS(6.87 × 10-10m),and characteristic time σHS/v0(1.94 × 10-12s),respectively.For the pore containing an isolated active site,as shown in Fig.1(a),the pore was cubic with side lengthL(which can be approximated as pore size) ranging from 40 to 140,which corresponds to mesoscale pores.The two ends along thexdirection were open.By symmetrizing the whole structure along thexdirection,the pore containing two adjacent active sites was constructed,as shown in Fig.1(b).Periodic boundary conditions (PBCs) were applied along thexdirection.That is,once a reactant or product molecule leaves one end of the pore,it will enter the other end as a new reactant,so the inverse reaction is not considered.The no-slip boundary condition was applied at the walls inyandzdirections.The red area in Fig.1 with side lengthls=20 on the upper wall represents the active site(s).As mentioned above,the reference Region S(the dotted lines in Fig.1)for investigating the fluctuation has a side length ofl=25.For the structure with two sites,only the concentration fluctuation near one site was investigated because of the symmetry of the structure.

At the initial stage of the simulation,all molecules were labeled as reactants and randomly distributed in the pore,and their random initial velocities were assigned according to the given temperature.For a given number density of all molecules(n),the effective packing fraction,η=can be obtained.Then,for a given pore sizeL,the diffusion efficiencyDeffand the reaction efficiencyReffcan be changed independently by changing η and the reaction efficientp0,respectively.

The simulation conditions and parameters refer to the actual reaction conditions of 2-butene isomerization into 1-butene over silica catalyst.The molar volumeVof butene was obtained by SRK (Soave-Redlich-Kwong) equation of state [39]:

where α（T）=,TandPare the temperature and pressure of the system,PcandTcare the critical pressure and critical temperature of the pure substances,ωis the acentric factor,andTr=T/Tcis the reduced temperature.For 1-butene[40],Tc=419.5 K,Pc=4.02 MPa,and ω=0.187,respectively,which are also used for 2-butene.Under the typical conditions of 423-523 K,4-7 MPa [41,42],as listed in Table 1,the range of η is from 0.0335 to 0.1551,which means that the interaction between molecules plays a dominant role rather than that between molecules and pore surface.The reaction efficientp0was assumed to be from 0.25 to 4.The time step was 1.0 and the results were stored at every 10th step.The number of reactant and/or product was sampled every 200 steps.The total simulation time was different for each simulation case,but it ensured that sufficient statistical results were obtained after the system reached steady state.

Table 1 The effective packing fraction η under typical T and P conditions

3.Results and Discussion

3.1.Reaction-diffusion coupling near an isolated active site

3.1.1.Characteristics of dynamic fluctuation in the concentration

In the simulations,when the statistical average of the products molecule number from the active site was equaled that diffused out of the pore in unit time,the system was considered to have reached the steady state.Then the averages over 100000 time steps shown in Fig.2 indicating thatReffincreases with the increase ofp0,but is almost independent of η.However,for a given pore sizeL,Deffdepends only on η and is hardly affected byp0.Therefore,the reaction and diffusion performance can be regulated independently,which provides a way to regulate the reaction-diffusion coupling.

To characterize the fluctuation,the whole simulation domain was meshed into many grids with side length of 3-5 (=depending on η) and the number of molecules (NrandNa) in each grid was counted.Then as mentioned in Section 2.3,to characterize how far the fluctuation is from complete randomness,the relative number density of the reactant and all molecules in Region S are introduced as.

respectively.With different reaction coefficientp0and effective packing fraction η,Fig.3 shows that the fluctuation ofNris much greater than that ofNafor almost all values ofp0and η.For largerp0and η,the fluctuation ofNris stronger,while that ofNais weaker and greatly affected by η,but hardly affected byp0.The results indicates that the fluctuation of the reactant concentration is different from the complete stochastic fluctuation due to the influence of reaction-diffusion coupling.

However,Fig.4 shows that whenDeffis large andReffis small,φr?ais minimized to about 1.0,indicating almost stochastic concentration fluctuation of reactants.With the decrease ofDeff,φr?aincreases up to 1.732,displaying non-stochastic concentration fluctuation.On the other hand,Reffhas little influence on φr?a.These results are qualitatively consistent with the predictions of our 2D theoretical model developed previously [20].However,all φr?aobtained in this work are smaller than the 2D model results with the samep0and η.This is because for the same η,the average distance between particles in 3D systemis larger than that in 2D systemwhich will reduce the diffusion resistance,and then weaken the concentration fluctuation.For example,whenp0=4,η=0.1551 andL=140,theDeffis 0.0018,about 1.2 times that in 2D model,and φr?ais 1.732,only about 35 % of that in 2D model.

To understand the reaction-diffusion coupling more intuitively,taking the simulation withp0=4,η=0.1551 andL=80 as an example,all reactants/products near the active site at a certain time (e.g.t=15,000) are tracked,and several projections of the positions of the tracked particles in thex-zplane at different times are given in Fig.5,in which the reactants,old and new products are distinguished.It is found that when a number of reactants reach the active site and react to form products,these products will prevent new reactants from approaching the active site,as shown in Fig.5(a).Only when these products diffuse away can new reactants reach the active site and react (as shown in Fig.5(b and c)).In addition,the changes of the centroid positions of the tracked reactants and products during this process are calculated respectively,in which the reactants and products that change in the reaction are not distinguished.The trajectory projections of the centroid positions in thex-zplane over a period of time are shown in Fig.6.It can be seen intuitively that the centroid position of the tracked reactants can reach the active site only after the centroid position of the tracked products is far away from the active site.The phenomenon that the products prevent the reactants from reaching the active site can also be observed when the conditions are changed.This steric hindrance effect results in the coupling of reaction and diffusion,which leads to the non-stochastic concentration fluctuation and complex dynamics near the active site.

Correspondingly,Fig.7 also suggests that φr?acan be correlated toalmost linearly.Based on the analysis above,for a fast reaction in a given pore structure,the local fluctuation is closely related to the gas diffusion resistance,which is mainly related to the number of molecules besides the influence of the pore structure.Therefore,φr?acan be roughly estimated bywithout detailed simulations.

Fig.2.(a) Variations of Reff versus η under different p0 and (b) Variations of Deff versus p0 under different η.

Fig.3.The variation of andunder different p0 and η.

3.1.2.Effect of the reaction-diffusion coupling on the overall performance

As mentioned above,the concentration fluctuation near an active site is not stochastic,but its magnitude depends on the reaction-diffusion coupling,and will further affect the overall performance (Y) of the pore structure.

Fig.4.The variation of φr?a versus Deff and Reff (the blue dots represent simulation data).

As shown in Fig.8(a),the typical pore sizeL=80 are used to study the effect ofDeffandReffon the yieldY,whenDeffremains unchanged,Yincreases with the increase ofReff.This is because with the increase ofReff,the reactants are easier to react at the active site and more products are produced,thus obtaining higherY.WhenReffremains constant,higher molecular number density(i.e.,higher η) leads to higher diffusion resistance and also the increase of reactant/product molecules,soYwill increase with the decrease ofDeffto some extent.AsReffincreases andDeffdecreases,the competition between the two factors eventually leads to the increase of both φr?aandY.It can be predicted that ifDeffdecreases further,that is,the molecular number density continues to increase,Ywill not increase further but reach a peak value,as shown by the dotted lines in Fig.8(b).However,the gas number density then will be higher than reasonable reaction conditions.Therefore,in the actual process of 2-butene isomerization to 1-butene,higherYwill be obtained under lower temperature and higher pressure.

3.1.3.Effect of pore size on the reaction-diffusion coupling

In the simulations,the pore size ranges from 40 to 140.Based on the results above,the typical parameters (η=0.1551,p0=4)leading to large φr?aandYare chosen to study the effect of pore size on the reaction-diffusion coupling.

Fig.5.The reaction-diffusion coupling caused by steric hindrance effect near the active site (p0=4,η=0.1551 and L=80,particles in red,green and magenta represent reactants,‘‘old” and ‘‘new” products,respectively).

Fig.6.The trajectory projections of the centroid positions of the tracked reactants and products near the active site during t=15000-15300 (p0=4,η=0.1551 and L=80,without considering the changes of reactants and products due to the reaction).

Fig.7.The relationship between φr?a and N-p/N-a.

As shown in Fig.9,with the increase of pore size (L),φr?aincreases whileYdecreases.To further understand it,the effect ofLonDeffandReffwas analyzed,as shown in Fig.10.It is shown thatReffis basically independent ofL,whileDeffis greatly affected byLfor a given η.By fitting the simulation results,it is found that the relationship betweenDeffandLsatisfiedDeff=104/L2.256.AsLincreases,theDeffdecreases rapidly.According to the definition ofDeff,for largerL,the proportion of molecules in the pore can effectively diffuse to the active site becomes small,which results in smallDeff.Therefore,with the increase ofL,theDeffdecreases,which leads to the increase of φr?a.SinceReffremains unchanged,the competition betweenDeffandReffeventually leads to the decrease ofY.The variation ofYversusLwas fitted toY=629.03/L3.24,as shown in Fig.9(b).From this point of view,to obtain higher yield,porous media with more relatively smaller pores should be used.It is noted that keeping η constant and increasingLleads to the decrease ofDeff,so thatYdecreases with the decrease ofDeff.This is different from the result that Y increases with the decrease ofDeffin Section 3.1.2,where keepingLconstant and increasing η leads to the decrease ofDeff.Of course,the diffusion performance of the whole catalyst particle should be further considered when the catalyst structure is optimized.

3.2.Reaction-diffusion coupling between two adjacent active sites

The discussion above shows that the reaction-diffusion coupling leads to the non-stochastic local concentration fluctuation near an isolated active site in the catalyst pore.However,in the actual catalytic system,the active sites will be non-uniformly distributed on the pore surface.There is also competition and compromise between the reaction and diffusion near adjacent active sites.

Similar to the pore structure with a single isolated active site,the length of the pore in thexdirection was 160,the cross section was square with side length 80.As shown in Fig.1(b),the two active sites were located on the upper surface of the pore and were symmetrically distributed about the center in thexdirection.The distance between the two sites wasd.Simulation conditions here were the same as those with the isolated site.In this section,the distancedvaried from 30 to 130,the typical parameters(η=0.1551,p0=4) were used and the effect ofdon the diffusion-reaction coupling was studied.

Fig.8.The variation of the yield Y versus Deff and Reff (The dashed lines in (b) represent the predicted values, L=80).

Fig.9.The variation of φr?a (a) and Y (b) versus L (η=0.1551, p0=4).

3.2.1.Effect of d on the reaction-diffusion coupling

Fig.11 shows thatReffis independent of the distanced,whileDeffincreases with the increase ofd,indicating that the proportion of molecules that can effectively diffuse to each active site becomes smaller.So,stronger concentration fluctuation is found for smallerdin Fig.12,indicating diffusion-reaction coupling between two adjacent active sites.With the increase ofd,their coupling gradually weakened,and the local concentration fluctuation near each site was closer to that of the isolated active site shown in Fig.9.

3.2.2.Effect of d on the overall performance

The coupling of reaction and diffusion between adjacent active sites will lead complex interactions among the components nearby.As shown in Fig.13,the yieldYincreases with the increase ofd.Compared with the results of the isolated active site with the same pore size,the reactions at the two adjacent active sites compete with each other,which leads to the decrease of reactant concentration near each site,thus making the yield of each active site slightly lower,but the total yield of the two sites is higher than that of the isolated sites.Therefore,for the simple reaction A →B,it is found that for the same pore size,the high yield of products will be obtained when the distance between adjacent sites is large.

However,the reactions on adjacent active sites are usually complex and involves multiple components and reaction steps.If sequential reactions occur on two adjacent sites respectively[43,44],for example,assuming that the reaction A →B occurs on one site and B →C on its neighbor site,then the variation ofYof product C versus the distancedis shown in Fig.14.It is found thatYincreases with the decrease ofd,which is just opposite to the reaction A →B discussed above.This is because the reduction ofdwill make the intermediates (B) easier to diffuse from one site to its neighbor site and the sequential reaction will occur,which is beneficial to the formation of the products and increases its yield.

Fig.10.The variation of Reff(line with square)and Deff(circle)versus the pore size L(η=0.1551, p0=4).

Fig.11.Variation of Reff (line with square) and Deff (line with circle) versus d.

Fig.12.The variation of φr?a versus d (η=0.1551, p0=4).

Based on the results of the two different reaction processes,the reactions transforming the same reactants at adjacent active sites compete with each other,and the larger distance between the two sites is beneficial to the overall performance.But the sequential reactions transforming different reactants at adjacent active sites are coordinated with each other,and the smaller distance is beneficial to the overall performance.Therefore,there will be an optimal value of the number density of active sites in a given pore structure.However,the present work is just based on very simple pore structure and reaction model.For more complex reaction pathways and pore structures,more work will be needed to obtain optimal active site distribution in the future.

Fig.13.The variation of Y versus d for reactions A →B on both sites.

Fig.14.The variation of Y versus d for sequential reaction A →B on the first site and B →C on its neighbor site.

At present,the concentration fluctuation due to the reactiondiffusion coupling near an isolated active site or even between two neighboring sites is still difficult to be observed directly in the experiment.It is conceivable that if the system contains a large number of regularly arranged sites,more complex coupling between these sites will occur.After careful regulation under different conditions,the ‘‘superposition” of concentration fluctuation between components may result in a great difference in the concentration or yield of specific reaction products at the macroscale.Then,it is expected to be observed and verified by experiments.

4.Conclusions

Molecular dynamics simulation combining hard-sphere and pseudo-particle modeling (HS-PPM) were carried out to study the reaction-diffusion coupling near active sites,and the steric effect between molecules was also discussed.Based on the non-stochastic concentration fluctuation of the reactant/product found near an active site with a simple reaction of A →B in a 2D pore model,two 3D pore models,containing an isolated active site and two adjacent active sites,respectively,were further studied in the present work.For a single isolated active site,it is found that the increase of the pore size leads to the increase of concentration fluctuation but the decrease of the product yield.When the pore size remains constant,a maximum product yield can be obtained by changing the reactant concentration.For two neighboring sites,the reaction performance is basically independent of the distance (d) between the two sites,but the diffusion performance is greatly affected byd.For the same reactions on the two neighboring sites,due to their competition,the largerdleads to the larger diffusion performance and overall performance.However,for the sequential reactions at the two sites,due to their coordination,the larger overall performance is obtained at smallerd.So,for a given pore structure,the optimal value of the active site distribution and high product yield can be obtained through optimization.Although the pore structure and reaction model in this work are relatively simple,the discussion of the nonstochastic concentration fluctuation near neighboring active sties may be helpful for optimizing the catalyst design and getting higher performance.Future work is needed to investigate the diffusion and reaction coupling in physical experiments and in the actual complex reaction pathways.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (92034302 and 22178347),the Dalian National Laboratory for Clean Energy (DNL) Cooperation Fund,the Chinese Academy of Sciences (DNL201905),the ‘‘Transformational Technologies for Clean Energy and Demonstration”,Strategic Priority Research Program of the Chinese Academy of Sciences(XDA21030700),the National Science and Technology Major Project (2017-I-0004-0005).

Nomenclature

Deffdiffusion efficiency

ddistance between two adjacent active sites

Llength of simulated pore (which can be approximated as pore size)

lslength of active site

llength of Region S for statistics

m1,m2mass of hard spheres

Natotal number of molecules entering the certain region in unit time

Npnumber of products entering the certain region in unit time

Nrnumber of reactants entering the certain region in unit time

Nrtime-averaged value ofNr

relative number density ofNr

Pccritical pressure,MPa

p0reaction coefficient

Reffreaction efficiency

Tabsolute temperature,K

Tccritical temperature,K

tone moment of evolution time

t0the moment after which the system reaches steady state

tnthe number of time steps fromt0tot

v0thermal velocity

vcthe critical velocity

vtvelocity at collision

Ythe yield of product

εLJdepth of the potential well in Lennard-Jones potential,J

η effective packing fraction

μafluctuation ofNa

μrfluctuation ofNr

σHSeffective diameter of hard sphere

σLJfinite distance at which the Lennard-Jones potential is

zero,10-10m

φr?arelative fluctuation ofNr