基于Gaussian和LAMMPS模拟的氧化钙脱除氯化氢机理

吴 畏 王统伟 汪德成 金保昇

(东南大学能源与环境学院，南京 210096)(东南大学能源热转换及其过程测控教育部重点实验室，南京 210096)

生活垃圾中由于含有大量的塑料制品，所以在焚烧过程中会产生大量有毒有害的氯化氢气体，目前通常的做法是采用氧化钙来脱除氯化氢[1-3].工程中常采用的方法有3种：①在垃圾焚烧的过程中往焚烧炉内喷洒氧化钙粉末，通常认为，当焚烧温度控制在850 ℃以上、焚烧时间超过2 s时，有害物质基本上都会分解掉[4].②在尾气排入大气前与脱硫脱硝等一起进行处理，通常安排在低温段(<100 ℃)，主要使用钙基吸附剂.其原理是酸性气体与吸附剂中活性组分发生酸碱中和反应，生成稳定的金属化合物而被固定下来[5-6].③在中高温时将烟气中的氯化氢脱除，以免当烟气的温度降到200～400 ℃时，烟气中的氯化氢在某些过渡金属的催化作用下生成二噁英[7-8].尾气中的二噁英脱除比较困难而且不可控，最好的方法是将生成二噁英的核心元素Cl提前脱除.氧化钙脱除氯化氢的机理通常被认为是酸碱中和反应，但是对其反应中的具体细节目前还缺乏深入的理论研究，尤其是在脱除效率方面.

采用氧化钙低温脱除氯化氢，其依据是中和反应理论，并有大量试验效果得到验证，但目前还没从原子、电子层面进行分析.焚烧炉中喷洒氧化钙粉末也是一种基于中和反应理论的经验做法，但还缺乏深入研究.中高温阶段深度脱除氯化氢时，总体认为钙基吸附剂在中高温段对氯化氢的脱除有一定的效果，但是大量的试验结果表明当温度超过某个限值时，脱氯效率会大大降低.Daoudi等[9-10]通过研究发现当温度大于650 ℃时，氯化氢的转化率会迅速减小；王蕊[11]通过试验研究发现600 ℃是反应的最佳温度；王婷等[12]的研究结果也表明在650 ℃时脱氯效率最大；刘金生等[13]通过离子色谱试验，得出脱氯的最佳反应温度为550～650 ℃；Shemwell等[14]则认为在固定床反应器中，500～600 ℃之间的脱氯效果最好.这与通常情况下化学反应速率随着温度的升高而升高的理论不符.因此，通过量子化学计算对氧化钙脱除氯化氢的机理进行深入研究具有重要意义.

Gaussian是一款基于量子化学的计算化学软件，它研究的对象是核外电子及其相互作用，可以计算过渡态能量和结构、化学键的生成和断裂、热力学性质、反应路径等，非常适合研究化学反应的机理.Zhao等[15]、Bian等[16]、Xu等[17]以及Hu等[18]利用Gaussian研究了所研究对象的反应机理、能垒、反应速率常数、指前因子以及产物的大气寿命等.LAMMPS是一款常用的分子模拟软件，重点考虑原子、分子结构以及相互作用，必要时通过反应力场来考虑化学反应的影响，能够很好地模拟分子结构的变化与特性.Zhang等[19]利用LAMMPS研究了二氧化碳的存在对氧化钙高温烧结的影响；Wang等[20]利用LAMMPS研究了水蒸气的存在对纳米氧化钙颗粒碳酸化过程的增强作用.Han等[21]利用LAMMPS对金属有机骨架的抗水稳定性进行了研究；Lu等[22]利用LAMMPS对聚酰亚胺的热解机理进行了研究，并研究了温度对其产物分布的影响.

本文首先采用Gaussian软件对氧化钙脱除氯化氢的多种可能反应路径进行探究，并且基于主要反应路径上的决速步分析了反应平衡常数以及反应速率常数；然后采用LAMMPS研究了2个氧化钙颗粒在不同温度下微观结构的改变.通过已有的试验结果和模拟结果相互验证，从而比较准确地揭示氧化钙脱除氯化氢的机理.

1 氧化钙和氯化氢的反应机理

1.1 模拟

本文中所有量子化学计算都采用Gaussian 09程序[23].使用MP2方法，在6-31G(d,p)水平上，对所有反应路径上的驻点，包括反应物、过渡态(TS)、中间体(IM)以及产物的几何构型进行了优化.在相同计算水平下，对所有驻点上的结构进行了振动频率计算，由频率分析表明，所有反应物、中间体以及产物的振动频率都是正值，所有过渡态有且仅有一个虚频，从而证实了结构的真实性.同时对每一个过渡态进行了内禀反应坐标(IRC)计算[24]，以确认每个过渡态所连接的反应物与产物.为了获得更精确的能量，在CCSD(full)/cc-pVTZ水平下计算了所有驻点的单点能(SPE).此外，为了描述反应过程中键的演变，还利用Multiwfn程序[25]计算了电子定域化函数(ELF)，并进行了分子内原子理论(AIM)分析.ELF是三维实空间函数，是一个用来衡量电子在局部区域定域化程度的指标，数值范围在0～1之间.ELF值较大的区域意味着该处电子定域性较强，电子不容易随意运动，这表明可能存在共价键、孤对电子或原子内壳层[26]；而ELF值较小的区域则表明电子在该区域很容易离域到其他区域.除了ELF外，通过分析AIM理论中键临界点(BCP)的性质也可以用来讨论成键特征[27].BCP处的电子密度(ρ)及其拉普拉斯值(▽2ρ)是键的拓扑性质，如果BCP处的▽2ρ符号为正，就认为该两原子间相互作用是闭壳层相互作用，如离子键、氢键、卤键之类；如果▽2ρ符号为负，则认为是共价作用，原因是由于共价相互作用的原子间会共享电子对，从而造成电子密度在成键区域聚集.此外，还计算了不同温度下的反应平衡常数，并且利用过渡态(TST)理论计算了不同温度下的反应速率常数，并考虑了一维隧穿效应的魏格纳校正[28-29].

1.2 氧化钙和氯化氢的反应路径

氧化钙脱氯反应在MP2/6-31G(d,p)水平上优化后的各个驻点的几何结构如图1所示.分步反应中各个驻点结构的单点能、总能量以及相对能量的数据如表1和表2所示，总能量的计算考虑了零点能(ZPE)校正，取相距无穷远的反应物的总能量作为参考零点.所有分步反应的势能剖面示意图(PES)如图2所示.

图1 反应物、中间体、过渡态和产物的几何结构图(单位：nm)

表1 第1步反应中各驻点的总能量与相对能量表kJ/mol

驻点结构单点能零点能校正总能量相对能量CaO-1 980 143.395.70-1 980 137.690HCl-1 212 131.2718.71-1 212 112.560P1-3 192 287.5627.15-3 192 260.41-10.15TS1-3 192 073.4615.40-3 192 058.06192.20IM1-3 192 355.6019.11-3 192 336.49-86.24TS2-3 192 204.6817.40-3 192 187.2862.98P2-3 192 441.4234.35-3 192 407.07-156.81TS3-3 192 406.0032.39-3 192 373.61-123.35CaClOH-3 192 841.7636.40-3 192 805.36-555.10

表2 第2步反应中各驻点的总能量与相对能量表

kJ/mol

对于CaO与HCl的反应机理研究，本文考虑CaO与1个HCl分子先反应，其反应产物随后再与第2个HCl分子反应.如图2(a)所示，CaO与HCl的第1步反应有2种反应路径：路径1，(CaO+HCl)→P1→TS1→IM1→TS2→CaClOH；路径2，(CaO+HCl)→P2→TS3→CaClOH.其中，路径1包含2个反应步骤：①相互独立的CaO的1Ca原子与HCl的3Cl原子不断靠近，在弱相互作用下直接形成稳定的反应前驱体P1，这一步不存在过渡态，所以该步为无能垒反应. ②随着1Ca-3Cl键键长的逐渐减小，形成了过渡态TS1，这个过程需要克服202.35 kJ/mol的能垒；TS1结构不稳定，4H原子会逐渐被1Ca原子所吸引，然后形成中间体IM1；随着2O-1Ca-4H键角的减小，IM1变成了过渡态TS2，在这个过程中需要克服149.22 kJ/mol的能垒；随着2O-1Ca-4H键角的继续减小，4H原

(a) 第1步

(b) 第2步

子又被2O原子所吸引，造成4H-1Ca键的断裂，同时生成新的4H-2O键，最终形成中间产物CaClOH.路径2只需要一个反应步骤：首先，相互独立的CaO的2O原子与HCl的4H原子在弱相互作用下形成反应前驱体P2，该步不存在过渡态，同样为无能垒反应；然后4H-3Cl键断裂，3Cl原子逐渐向1Ca原子转移，形成过渡态TS3，这个过程只需要克服33.46 kJ/mol的能垒；随着3Cl原子继续向1Ca原子转移，生成了新的1Ca-3Cl键，同样也形成了中间产物CaClOH.从能量角度看，路径1中最大的能垒为202.35 kJ/mol，远大于路径2中的33.46 kJ/mol，如此大的能垒也表明该反应步骤是不太可能发生的，所以在CaO与HCl的第1步反应的2种反应路径中，路径2更可能发生.

如图2(b)所示，CaO与HCl的第2步反应为路径3：(CaClOH+HCl)→P3→TS4→IM2→(CaCl2+H2O).路径3中只包含一个反应步骤：首先，第1步反应形成的中间产物CaClOH与第2个HCl分子在无能垒反应下形成反应前驱体P3；然后6H原子往垂直纸面向里的方向旋转，形成过渡态TS4，这个过程需要克服26.13 kJ/mol的能垒；随着6H原子的继续转移，最终与2O原子结合，生成新的2O-6H键，此时结构变为中间体IM2，也称为反应后驱体；随后1Ca-2O键断裂，不经历过渡态直接形成最终产物CaCl2和H2O，这一步同样是无能垒反应.

从以上分析可知，氧化钙脱氯机理的主要反应路径是首先通过路径2生成中间产物CaClOH，然后通过路径3生成CaCl2·H2O，最终分解为CaCl2和H2O.关于中间产物CaClOH，在试验中也被证实是真实存在的[11]，因此量子化学计算可以用来寻找在试验中因为存在时间很短而很难被观测到的反应中间产物.上述反应路径中的决速步是P2→TS3→CaClOH，该步反应计算得到的活化能垒为33.46 kJ/mol，与试验测得的28.1[30]、31.1[31]和45.0 kJ/mol[32]的均值34.7 kJ/mol仅相差3.54%.因此，氧化钙与氯化氢的反应可以写成如下形式：CaO+HCl→CaClOH，CaClOH+HCl→CaCl2·H2O，CaCl2·H2O→CaCl2+H2O，该反应形式与Gullett等[30]通过试验得出的结论完全一致，因此本文模拟分析得到的主要反应路径就是试验中最可能发生的反应路径.

1.3 键的演变分析

为了更深入地了解反应机理，本节对氧化钙脱氯反应中主要反应路径上的驻点结构进行了ELF和AIM分析.氧化钙脱氯反应中主要反应路径上，对应驻点结构在BCP处的AIM参数如表3所示，驻点结构的ELF填色图见图3.

由图3可知，CaO中1Ca-2O键的ELF值很小，意味着1Ca-2O键主要是离子键.从表3的AIM分析中可知，1Ca-2O键在BCP处的电子密度为0.091，并且▽2ρ是一个很小的正值，由此得出的结论与ELF分析的一致.值得注意的是，在HCl的AIM分析中，在4H-3Cl键之间没有找到BCP，这主要是由于4H-3Cl键是极性共价键，使得BCP落到H原子核内，所以表现为没有找到BCP，但是从ELF填色图可以看出，4H-3Cl键之间的ELF值接近于1，证明4H-3Cl键确实是共价键.在P2中，4H-3Cl键的▽2ρ为0.053 au，意味着4H-3Cl之间的共价键断裂，成为离子键；而2O-4H键的▽2ρ为-2.333 au，说明随着2O原子和4H原子的不断靠近，两者之间形成了共价键.在TS3中，4H-3Cl键的BCP消失并且两者之间的ELF值也变得很小，同时3Cl-2O键之间出现新的BCP且▽2ρ等于0.279 au，意味着在这过程中3Cl原子发生了转移.在CaClOH中，1Ca-3Cl键在BCP处的电子密度由0.020变为0.042，▽2ρ由0.072 au变为0.177 au，说明在这过程中1Ca-3Cl键的强度增强了.在CaO与HCl第2步反应的前驱体P3中，1Ca-5Cl键的电子密度为0.014，▽2ρ为0.064 au，说明CaClOH与HCl之间主要是离子键相互作用.随着6H原子的转移，6H-3Cl和1Ca-6H之间的离子键发生断裂，表现为在TS4中，6H-3Cl和1Ca-6H键之间的BCP消失.随着6H原子的继续转移，在IM2中的2O-6H键之间出现新的BCP，该处的电子密度为0.355，▽2ρ为-2.795 au，表明6H原子与2O原子之间形成了共价键，而该键之间的ELF值也接近于1，同样说明了该键主要是共价键.最终反应后驱体IM2在无能垒反应下，发生1Ca-2O键的断裂，生成最终产物CaCl2和H2O，从ELF值和AIM参数可知，CaCl2中主要为离子键，H2O中主要为共价键.

表3 CaO+2HCl主要反应路径上驻点结构在BCP处电荷密度的拓扑特性

(a) CaO

(b) HCl

(c) P2

(d) TS3

(e) CaClOH

(f) P3

(g) TS4

(h) IM2

(i) CaCl2

(j) H2O

图3 CaO+2HCl主要反应路径上驻点结构的ELF填色图

2 反应平衡常数和反应速率常数分析

2.1 反应平衡常数

由于氧化钙脱氯反应的决速步是在第1步反应中，所以本节仅对决速步反应的平衡常数进行计算，结果如表4所示.

平衡常数的计算公式如下：

lnK=-ΔG/(RT)

(1)

式中，ΔG为吉布斯自由能变，kJ/mol；R为理想气体常数，J/(mol·K)；T为开尔文温度，K.由表4可知，ΔG值随着温度的升高而缓慢增大，lnK随温度的升高不断减小，说明当温度升高时，虽然脱氯反应依然会发生，但是向正反应方向进行的程度会受到抑制，反应逐渐向逆反应方向移动，该结论与万旦[33]的研究结果一致.

表4 不同温度下决速步的反应平衡常数

2.2 反应速率常数

本节通过反应速率常数的计算来分析温度对正向与逆向反应的影响.

采用传统的TST理论[28]计算速率常数，公式如下：

k(T)=κ(T)kTST

(2)

(3)

式中，kTST为由TST理论计算得到的近似速率常数，s-1；κ(T)为隧穿校正系数；kB为玻尔兹曼常数，J/K；h为普朗克常数，J·s；QTS(T)为过渡态的配分函数，上标TS表示过渡态；QR(T)为反应物的配分函数，上标R表示反应物；ΔE为活化能垒高度，kJ/mol.

κ(T)采用魏格纳隧穿校正进行估计[29]，公式如下：

(4)

式中，Vm为过渡态的虚频，cm-1；c为光速，m/s.

QTS(T)和QR(T)的计算公式如下：

Qj=QtransQrotQvibQelectj=TS,R

(5)

式中，Qtrans、Qrot、Qvib和Qelect分别为平动配分函数、转动配分函数、振动配分函数和电子配分函数.这些参数可以在Gaussian的输出文件中直接获取，不需要额外的计算.

温度298～1 100 K之间的决速步正向与逆向的反应速率常数如表5所示，不同温度下正向与逆向反应速率常数对数值的变化曲线如图4所示.由图可见，随着温度的升高，氧化钙脱氯反应的正向反应速率常数与逆向反应速率常数都呈增长趋势，其中逆向反应速率常数的增长速度要大于正向的增长速度，说明温度的升高会抑制正向反应的进行，这与反应平衡常数分析得到的结论相一致.但是从表5可以看出，虽然温度的升高确实会抑制正向反应的进行，但是正向反应速率常数始终远大于逆向反应速率常数，所以正向反应速率常数在整体反应速率常数中占支配地位.因此从热力学的角度来看，造成温度升高脱氯效果降低的主要原因并不是由于逆反应的增强导致的整体反应速率的下降，同样的规律也可以在万旦[33]的研究中发现.

表5 不同温度下决速步正向与逆向的反应速率常数

图4 不同温度下决速步正向与逆向的反应速率常数变化曲线

3 不同温度下氧化钙的微观结构

3.1 模拟

量子化学计算的结果表明，氧化钙脱氯反应的逆反应不是造成温度升高脱氯效果降低的主要原因.为了探究其背后的真正原因，本文采用LAMMPS软件对氧化钙与氯化氢进行了分子动力学(MD)模拟.为了模拟不同温度下氧化钙微观结构的变化，采用了反应力场(ReaxFF)进行模拟，关于ReaxFF势函数中每项能量的详细描述见文献[34-35]，Ca/O/H/Cl的参数来自于Psofogiannakis等[36]对ReaxFF的开发.

MD模拟的具体步骤如下：①利用氧化钙的单胞制备出2个球形氧化钙颗粒，颗粒半径取1.2 nm，沿垂直z方向放置，如图5所示，球形颗粒间的距离取0.9 nm；②将2个颗粒放置于包含200个随机生成的氯化氢分子的模拟盒子中，盒子尺寸为5 nm×5 nm×8 nm；③将初始系统能量最小化以消除所有可能的重叠和紧密接触；④采用NVT系综在目标温度下将平衡系统运行60 ps，步长取0.1 fs，边界条件为周期性边界；⑤每隔0.2 ps将模拟轨迹输出，分析系统的特性(能量、温度、压力、原子坐标、速度等).上述参数的选取来自于Zhang等[19]的研究，已经被证实可以很好地再现氧化钙的烧结过程.

图5 包含2个氧化钙颗粒与200个氯化氢分子的模拟盒子示意图(单位：nm)

3.2 结果与讨论

模拟盒子在不同温度下运行60 ps后的微观结构变化如图6所示.由图可以看出，当温度在298～700 K之间时，随着温度的升高，氧化钙的微观结构由最初的有序晶体结构逐渐变为无序的非晶体结构，并且随着温度的升高，氧化钙颗粒的体积也变得比初始结构的体积大.王蕊等[37]在试验中也观察到，高温下的氧化钙表面会变得更加疏松，颗粒表面呈现松碎的形态.正是由于这种松碎的形态，使得颗粒的表面积变得更大，从而与氯化氢气体的接触面积也越大，同时疏松的结构会使得氯化氢气体向内部扩散的阻力变小，使得反应在高温下更容易进行，Duo等[38]的研究也得出了类似的结论.而且随着温度的升高，2个氧化钙颗粒之间的距离在不断减小，但是在298～700 K之间，氧化钙颗粒尚未发生烧结，此时氧化钙的脱氯效率主要由化学反应主导，因此表现为脱氯效率随着温度的升高而增加；当温度在700～900 K之间时，可以看到2个颗粒已经开始发生初步烧结，虽然烧结会导致传质阻力的增加，但温度升高的同时传质阻力会相应减小，因此当温度在700～900 K之间时，脱氯效率整体还是呈上升趋势；但是当温度大于900 K时，颗粒间的颈部区域已经完全被反应物与生成产物填满，意味着温度在900～1 100 K时氧化钙已经完全烧结.在Weinell等[39]的研究中也发现当温度超过900 K时，氧化钙与氯化氢的结合能力受到了限制.随着温度的继续升高，烧结现象也越来越严重.由于颗粒发生烧结，造成反应的表面积大大减小，导致传质阻力增加，此时温度升高对传质阻力的减小已经不能弥补由于烧结带来的传质阻力的增加，因此氧化钙的脱氯效率在温度大于900 K时出现下降的趋势.Daoudi等[9-10]和王恺等[40-41]通过研究同样认为氧化钙脱氯效果下降的主要原因是由于高温下材料结构烧结导致传质阻力的增加.

(a) 298 K

(b) 400 K

(c) 500 K

(d) 600 K

(e) 700

(f) 800 K

(g) 900 K

(h) 1 000 K

(i) 1 100 K

图6 2个氧化钙颗粒(相距0.9 nm)在不同温度下煅烧60 ps后微观结构变化示意图

4 结论

1) 氧化钙脱除氯化氢的主要反应路径是CaO+HCl→CaClOH，CaClOH+HCl→CaCl2·H2O，CaCl2·H2O→CaCl2+H2O.

2) 反应平衡常数随着温度的升高而减小.当温度升高时，向正反应方向进行的程度会受到抑制，反应逐渐向逆反应方向移动.

3) 反应速率常数随着温度的升高而增加，但是正向反应速率常数始终远大于逆向反应速率常数，因此造成温度升高脱氯效果降低的主要原因并不是由于逆反应的增强导致的整体反应速率的下降.

4) 当温度在298～700 K之间时，化学反应占主导地位，脱氯效率随温度的升高而增加；当温度在700～900 K之间，氧化钙颗粒开始烧结，但此时温度的升高可以克服由于初步烧结导致增大的传质阻力，整体的脱氯效率依然呈上升趋势；当温度大于900 K时，颗粒完全烧结，此时温度的升高无法弥补由烧结产生的传质阻力，从而导致氧化钙的脱氯效率随着温度的升高而下降.