张金利 何正华 韩 优,2,* 李 韡 武江洁星 甘中学 谷俊杰

(1天津大学化工学院,天津300072;2天津大学,天津市膜科学与海水淡化技术重点实验室,天津300072; 3新奥集团煤基低碳能源国家重点实验室,河北廊坊065001)

超临界水中碳酸钠团簇成核与生长分子动力学模拟

张金利1何正华1韩 优1,2,*李 韡1武江洁星1甘中学3,*谷俊杰3

(1天津大学化工学院,天津300072;2天津大学,天津市膜科学与海水淡化技术重点实验室,天津300072;3新奥集团煤基低碳能源国家重点实验室,河北廊坊065001)

应用分子动力学方法研究了碳酸钠颗粒在超临界水中的成核与生长过程.计算了温度为700-1100 K、压力在23-30 MPa下碳酸钠的团聚过程,计算时间为1 ns.对体系结合能与径向分布函数的分析表明,碳酸钠成核过程主要受静电作用的影响.在超临界态下,水分子与Na+和之间的静电作用降低,Na+与能够很容易碰撞形成Na2CO3小团簇.在Na2CO3整个成核过程中,单个离子的碰撞在前50 ps内完成,同时离子碰撞速率达到1030cm-3·s-1.另外,在成核阶段温度的影响比压力更加明显,温度越高,离子碰撞速率越快,形成的初始团簇越多.而压力对Na2CO3团簇的进一步生长影响较大.

超临界水;碳酸钠;结合能;碰撞速率;分子动力学

1 Introduction

Crystallization of salts is a phenomenon of great practical relevance.In fact,it is also one of the most important methods for industrial separation.In recent years,the nucleation of salts in supercritical water,as well as the properties of those aqueous systems under supercritical conditions,has received increasing interest because of its importance in hazardous organic waste water treatment,1natural geothermal processes,2design of the new materials,3,4and supercritical water oxidation.5-8

In the supercritical water(critical temperature:647 K and critical pressure:22.1 MPa),the salt solubility would change obviously.Therefore,supercritical water as an important solvent in salt nucleation research received extensive attention. Svishchev et al.9investigated the nucleation of strontium chloride nano-particles in supercritical water by molecular dynamic simulations.Their results showed that water molecules could reside both on the surface and in the interior of SrCl clusters, and the distribution of the clusters showed a very strong dependence on the density of the system.In their other work,10they examined the NaCl nano-particles generated by a rapid quench of supercritical aqueous solutions,and analyzed the spectroscopic signatures of NaCl clusters with different sizes and estimated their relative stabilities.Lümmen and Kvamme11,12studied the aggregation of ferrous chloride,the growth and properties of ferrous chloride with sodium chloride nano-particles in supercritical water by molecular dynamic simulations.They found that the growth rate of FeCl2particles was affected by the temperature and the density of the system.The growth rate of FeCl2particles was faster with lower temperature and higher density of the system.Finally,the disordered FeCl2crystalline was formed.

Recently,the supercritical water has been introduced in the reaction of coal catalytic gasification to enhance the efficiency.13,14The experimental results showed the Na2CO3catalyst had higher catalytic activity under supercritical water than non-supercritical water for the coal gasification reaction.But the effect mechanism of supercritical water on the catalyst as well as its contribution during the reaction process has been unclear yet.In fact,the nucleation and growth process of the Na2CO3particles in the supercritical water depend on their shape,size distribution,and loading status,which are crucial for the whole catalytic reaction process.But the nucleation process of Na2CO3particles in the supercritical water is difficult to be analyzed by using experimental method.However,with the developments of computer technology and molecular dynamics theory,theoretical studies on the micro-structure of solution system as well as the interaction mechanism of different moleculeswith moleculardynamicssimulation method have been carried out extensively.15-19Therefore,the nucleation and growth process of Na2CO3particles in supercritical water was studied in this work using molecular dynamics simulation method.The binding energy and radial distribution function (RDF)were introduced to analyze the interaction mechanism of water and Na2CO3.The effects of temperature and pressure on the distribution of Na2CO3particles as well as on ions collision rate were discussed.The important conclusions of Na2CO3nucleation in the supercritical water will shed light on the real production process.

2 Computational methods

2.1 Simulation details

The Forcite package20in the Materials Studio software21was applied here for the whole molecular dynamics simulations. The force filed of compass26,22which is suitable for the supercritical water system,23was chosen for the aqueous system.The O－H sp3hybrid orbital energy model of the o2*forcefield type,which provides a good prediction of water construction and properties,was assigned for the O atom in water molecules.The h1o forcefield type,which describes the bond between H and O,is suitable for the H atom in water molecules. The c3i forcefield type provides a good description of the C atom in the.The o2c forcefield type,which describes the O－X sp3hybrid orbital energy model in acid,is suitable for the O atom in the.The an+forcefield type was assigned for Na+.

The equations of motion were solved by Verlet leapfrog algorithm,24with the time step of 0.5 fs.Simulations have been performed in the NPT ensemble(Isobaric-Isothermal conditions). In order to determine a reasonable method to control the temperature and pressure during the dynamics simulation process, four temperature control methods(Andersen,25Berendsen,26Nose,27and Velocity28)combined with two pressure control methods(Andersen29and Berendsen25)were tested.When Andersen thermostatand Berendsen barostat were chosen,the temperature and pressure of the system can be quickly stabilized within 10 ps.Therefore,the temperature and the pressure were controlled by Andersen thermostat and Berendsen barostat methods during the whole simulation.Electrostatic interaction and van der Waals interaction were calculated using Ewald method.30The calculated density of supercritical water at 673 K and 28 MPa by the above method is 0.246 g·cm-3.It is approximate to the experimental value(0.25 g·cm-3),31indicating that the simulation method is reasonable.

A series of simulations ranging from 700 to 1100 K and from 23 to 30 MPa were explored.In each simulation,1024 water molecules,60 Na+,and 30(14.7%(w,mass fraction)Na2CO3)were randomly distributed in the simulation cell. And then,the geometry of this system was optimized with smart minimization method(which automatically combines appropriate calculation methods in a cascade.It starts with the steepest descent method,followed by the conjugate gradient method and ends with a Newton method)to prevent formation of the initial ion pairs and the unreasonably overlap with water molecules.After that,the dynamics simulations of the relaxed system at different state conditions in the supercritical region were carried out.The total simulation time was 1000 ps,and one sample could be recorded every 1000 simulation steps for follow-up analysis.

2.2 Binding energy

Based on the framework of the classical nucleation theory (CNT),32,33the whole(Gibbs)free energy of the nucleation system increased at first and then decreased as the cluster radius increased.Thus the nucleation is not a spontaneous process.It requires some impetus to overcome the energy barrier.Generally,the impetus comes from the spontaneous fluctuation of density or component in the metastable phase.34According to the CNT,a certain energy barrier is present during the nucleation process.Once the system has enough impetus,the energy barrier would be overcome and a new phase would be formed.At this moment,the thermodynamic energy is changed to the kinetic energy.As the temperature rising,the molecules move faster,which may lead single molecules to collide with each other,and then the total energy of system decreases with the new particles appearing.But the interaction of water molecules with Na2CO3is the main obstruction for the above process. Therefore,investigation of the binding energy of water and sodium carbonate is useful for better understanding of the mechanism of Na2CO3nucleation in supercritical water.

The binding energy is defined by the following equation,

where Ebindrepresents the binding energy,EH2OandENa2CO3are the energies of water and sodium carbonate in the Na2CO3solution,respectively,andEH2O+Na2CO3is the energy of the Na2CO3solution system.All the above energies are calculated in the same conditions.In the Forcite package,EH2O,ENa2CO3,Ekin,and EH2O+Na2CO3are contributed by their kinetic energy(Ekin)and potential energy(Epot),and the potential energy is defined as,

where Enon-bond,Evan,Eelect,EH-bond,Ecross,and Evalenceare the nonbond energy,van der Waals energy,electrostatic energy,hydrogen bond energy,cross term energy(which comes from the valence bond stretch,bend,and torsion),and valence interaction energy,respectively.Under the supercritical condition,the interaction of hydrogen bonds between water molecules is weak, so the hydrogen bond energy is very low,which could be ignored over the whole simulation ranges.Thus,the equation(3) was simplified as,

Therefore,

where ΔE represents the change of corresponding energies, which has been obtained by measuring the energy difference between the pure water,Na2CO3,and the Na2CO3solution at the same condition.

2.3 Cluster size

The size of Na2CO3clusters can be determined using Stillinger criteria.35That is,any two atoms belong to one cluster if their distance is less than 0.32 nm(which corresponds to the average value of the first minima peak of the Na+-Na+,Na+-,andradial distribution functions at melted state).If a continuous path which has been mentioned above is present between two atoms,they belong to the same cluster.

2.4 Collision rate

The nucleation rate is defined as the number of the clusters which is larger than the number of critical nucleus generated per unit time unit volume.The size of critical nucleus could estimate with the kinetic method reported by Yasuoka and Matsumoto.36Firstly,one threshold should be chosen.The threshold is a specific number corresponding to the number of atoms or ions in the cluster.If the atom numbers of a cluster exceed the threshold,the number of the clusters is plotted versus the simulation time,and a growth-decay evolution curve is produced.There will be a linear region at the beginning simulation time,and the slope of the linear region would reduce with the increasing threshold.When the threshold exceeds the critical value,the slope would stop decreasing.37The slope divided by the lattice volume is the nucleation rate.

3 Result and discussion

3.1 Interactions between water molecules and ions

The binding energy of water and sodium carbonate at each state point(temperature:700-1100 K,pressure:23-30 MPa) was calculated here and shown in Fig.1(a).It clearly showed that the impact of temperature on the binding energy was more obvious than that of pressure.When the pressure was constant, the binding energy decreased with temperature increasing.Especially,when the temperature was lower than 900 K,Ebindwas reduced quickly,whereas it declined slower when the temperature was higher than 900 K.The result implied that the interaction between water molecules and Na2CO3was reduced as temperature increasing.

In section 2,we introduced that the binding energy was contributed by the changes of kinetic energy,van der waals energy,electrostatic energy,the cross term energy,and the valence interaction energy(shown in formula(5)).During the simulation process,the energies of the cross terms and the valence interaction have no changes.Meanwhile,the changes of the kinetic energy and van der Waals energy are too small compared with the change of electrostatic energy.Thus,the binding energy of water and Na2CO3is mainly contributed by the change of electrostatic energy,which is shown in the Fig.1(b).That is, the electrostatic interaction between water and Na2CO3decreases at higher temperature.

Fig.1 Binding energy and the change of electrostatic energy at different temperatures and different pressures

In order to analyze the impact of interaction between water and Na2CO3on the Na2CO3nucleation process,ΔEpot,ΔEkin,and ΔEelectwere calculated during whole simulation process,their trends with simulation time at different temperatures and pressures are shown in Fig.2.Again,ΔEpotand ΔEelecthas no obvious dependence on the pressure(shown in Fig.2(a)),whereas they become lower at higher temperature(shown in Fig.2(b)) during the whole process.Besides,the change of kinetic energy almost keeps at 0 kJ·mol-1during the simulation process. The changes of potential energy and electrostatic energy decrease quickly from～800 to～150 kJ·mol-1at the initial 100 ps,indicating that the nucleation and particle initial collision of Na2CO3in supercritical water mainly occurres during this period.This period is shorter than that of FeCl2clusters,which is～200 ps.38

乔木林(纯林和混交林)碳储量3337416 t，占总碳储量的82.45%。按地类分：纯林碳储量3143808 t，混交林碳储量193608 t，纯林碳储量是混交林碳储量的16.24倍。从各龄组的碳储量分析：以中龄林、近熟林碳储量为主，其碳储量之和为2566391 t，占纯林碳储量的76.90%；中龄林的碳储量最大，其碳储量1606955 t，占纯林碳储量的48.15%。……