基于分子动力学模拟的轮胎橡胶催化热解制氢机理

陶礼，杨启容，李昭莹，亓昊，王力伟，马欣如

（青岛大学机电工程学院，山东 青岛 266071）

近年来我国能源消费全面增长导致能源的需求量不断增加，传统的一次能源已经无法满足需求，同时传统能源带来的环境问题也不容忽视，因此寻求新型的清洁能源成为刻不容缓的任务。其中氢气作为一种多功能清洁能源载体，在能源安全和环境保护方面已引起全世界关注。氢具有显著的优势，如零碳、高效、能量互连介质、可以无损失地存储，此外在燃烧过程中氢仅释放水，安全环保。但是传统制氢方法的原料不可持续且碳排放量高，太阳能和电解水制氢存在能耗成本高的缺点。因此以固体废物为原材料进行热解制氢，既能缓解目前固体废物处理的紧迫性，又能使资源得到合理化利用。轮胎橡胶作为固体废物的重要类型，具有高碳氢比，是热解制氢的可选材料。

目前，研究者从热解终温、升温速率、热解压力和催化剂等几方面对废橡胶的热解进行了实验研究。为了使目标产物氢气含量增加，也有许多学者通过实验研究不同催化剂的效果。Li 等采用浸渍法制备了六种不同类型的催化剂（镍、铁和钴，分别由γ-AlO和活性炭负载）热解轮胎，实验结果表明对于相同的催化剂载体，Ni 催化剂的氢产率最高， 其次是Fe 催化剂和Co 催化剂。Arabiourrutia 等分析了HZSM-5和HY催化剂在流化床反应器中的轮胎热解，这些沸石催化剂增加了热解气体产量。Yu等采用初湿浸渍法制备了Cu/HZSM-5、Zn/HZSM-5、Cu/MCM-41和Zn/MCM-41催化剂对橡胶进行热解，由于过渡金属引入的路易斯酸性位点的高脱氢活性，锌和铜的结合增加了气体产量，并显著促进了H的形成。现阶段主要通过实验手段对轮胎橡胶进行热解研究，但是热解过程复杂，宏观试验难以解释热解机理，分子动力学模拟以经典牛顿力学为基础，能够分析分子间作用力、分子形态等微观性质的变化，可以作为宏观试验的补充。Liu 等利用反应分子动力学（RMD）和密度泛函理论（DFT）对甲基环己烷的热解进行了研究。发现2300～3100K的温度下，从RMD模拟中提取的表观活化能为263.60kJ/mol，这与实验结果相当一致。

因此本文利用Material Studio 软件，选取Ni、ZSM-5 以及Ni/ZSM-5 催化剂，在研究出的轮胎橡胶生成氢气路径的基础上对每条催化路径进行热力学与动力学计算，选取制氢较好的催化剂进行热解模拟，并且探究不同比例的催化热解过程及其产物。

1 模拟方法

1.1 MS模型建立及优化

轮胎橡胶主要由天然橡胶（NR）、丁苯橡胶（SBR）和顺丁橡胶（BR）组成。三者的比例可根据不同轮胎型号而定，为了使模型能够与实际情况结合，根据文献以及之前已做过的研究，选取三者橡胶比例为NR∶SBR∶BR=5∶4∶1，并根据此比例进行模型建立。利用Materials Studuio 软件中的Visualizer 模块分别建立三种橡胶单体模型，经过几何优化使模型的能量达到最低，具有合理化的结构，根据热力学原理，当体系为稳定平衡状态时，其热力学势具有最小值。然后将单体拓展成长链。三个橡胶的主要成分分子式、单链及长链模型见表1。

表1 橡胶主要成分分子式、单链及长链模型

在Amorphous Cell 模块中将三种橡胶长链按照5∶4∶1的比例在298K温度下聚合在一起，选用周期性边界条件，包含400个C原子，606个H原子，总大小为1006个原子，立方体模拟盒边长为20.56Å（1Å=0.1nm），密度为1.034g/cm。为了使模型结构更加合理化，利用Forcite等模块对聚合成功的模型进行几何优化、退火循环等操作。退火是将分子“加热”到较高的温度，让分子结构在较高温度下充分舒展和松弛，之后让其冷却下来，目的是试图通过利用每个退火周期的较高温度周期来克服能量势垒，从而找到不同的局部最小能量结构，这样就可以达到全局能量最低点获得三维最优结构，以达到没有虚频的最稳定状态。虚频是指能量二阶导数矩阵（Hessian矩阵）的负数本征值，在这个振动方向上，体系结构处于能量的最大值，因此有虚频的结构都是不稳定结构。选取300～500K之间的温度进行退火循环，采用NPT系综，最终获得稳定结构的模型，见图1。建立Ni、ZSM-5、Ni/ZSM-5（Ni与ZSM-5 的质量比为1%、5%、10%）路径模拟、热解轮胎橡胶模型时，首先从Materials Studio 软件中的模型库导入Ni、ZSM-5模型，见图2、图3。将轮胎橡胶模型与Ni催化剂模型按照1∶1、3∶1以及5∶1的质量比例混合后也通过Forcite等模块进行退火和几何优化，其模型见图4。

图1 优化结构后轮胎橡胶模型

图2 Ni催化剂模型

图3 ZSM-5催化剂模型

图4 加入1∶1 Ni催化剂轮胎橡胶模型

1.2 路径模拟

在研究出来的三种橡胶的反应路径基础上，选取生成氢气的路径对其进行催化模拟研究，判断该催化剂是否利于生成氢气。利用密度泛函理论进行路径模拟前，为了使生成物以及反应物的结构模型较稳定，对其进行几何优化，该过程利用MS软件中DMol3 模块，基于GGA、BP 密度泛函数对橡胶三种单体分别进行结构优化。将优化后的路径原子进行配对，然后通过LST/QST（线性同步变换/二次同步变换）工具进行过渡态搜索。过渡态是以量子力学对反应过程中的能量变化的研究为依据，认为从反应物到生成物之间形成了势能较高的活化络合物，活化络合物所处的状态叫过渡态。通过过渡态搜索能够得到系统的热力学、动力学参数。经过计算分析可以得出路径的能垒、焓值以及吉布斯自由能等参数。几何优化获得的结果给出了系统的总电子和离子能（）。热焓变见式(1)。

系统的内能（）来源于电子的振动、平移和旋转的贡献，见式(2)。

吉布斯自由能差见式(3)。

∆见式(4)。

反应进行的能垒∆见式(5)。

式中，分别表示过渡态的焓值、反应物的焓值、反应物的分子系数、气体常数、温度。一旦基态结构被确定了，并且使用统计力学获得固有模式的振动频率，利用DMol3模块计算上面的参数。能垒越大越不容易形成活化的中间产物，反应难以进行。吉布斯自由能的变化可作为恒温、恒压过程自发与平衡的判据。在等温、等压且不做非膨胀功的条件下：∆< 0，表示反应能自发进行；∆= 0，表示系统已处于平衡状态，反应不进行；∆> 0，表示反应不能自发进行。焓变反应，当∆< 0，反应放热；∆> 0，反应吸热。

1.3 热解模拟

通过MS 软件的CULP 模块，选取ReaxFF 力场，采用NVT系综分别对橡胶进行催化热解模拟，时间步长为0.1fs，总时间为100ps，时间极短。为了使短时间内模型能够接受足够的能量，设置热解温度会高于实际热解温度，但这并不会影响反应机理。温度从500K 开始，以200K 为间隔单位进行热解模拟，直至橡胶完全断裂成小分子产物，不存在大量的大分子产物。观察记录产物分布，分析不同比例催化剂对热解模拟的影响。ReaxFF 是一种经验反作用力，它允许从非键合相互作用到键合相互作用的平稳过渡。由电负性平衡法（EEM）计算的原子间距离和部分电荷计算的键级在每个时间步都会更新，从而允许在没有任何预先定义的反应路径的情况下动态描述化学反应。

2 结果与讨论

2.1 路径模拟

针对三种橡胶热解产物进行研究，发现热解大致分为低温和高温两个阶段。在低温热解阶段，系统当中存在的大量基团是天然橡胶产生的异戊二烯（CH）单体、丁苯橡胶产生的苯乙烯（CH）和1,3-丁二烯（CH）单体以及顺丁橡胶产生的1,3-丁二烯（CH）单体，进入高温热解阶段，其过程是多步自由基反应生成，伴随着自由基的生成、传递和消失，得到目标产物。在高温阶段初期，主要是各单体裂解成自由基，存在的小分子自由基主要有·CH、H·、CH==CH·。随着温度继续升高，自由基攻击单体生成小分子气体以及小分子结构，小分子气体主要有CH、H、CH、CH、CH。本文作者课题组针对三种橡胶产物找到了链引发和链传递的反应路径，见图5～图7，这些反应路径是单体生成自由基路径。为了研究催化剂对生成氢气的影响，选取H 基攻击单体的路径，对加入Ni、ZSM-5、Ni/ZSM-5催化剂进行路径模拟。

图5 异戊二烯产生自由基的反应路径

图6 1,3-丁二烯产生自由基的反应路径

路径1～路径12 经过动力学计算后的反应能垒已经标在图1～图7中，从图中发现，最容易释放H基的单体是1,3-丁二烯，苯乙烯单体所有释放H基的路径能垒都比较高，说明苯乙烯比较难释放自由基。比较路径1～路径12 所有的能垒，发现最低的是路径3异戊二烯释放甲基。自由基生成后，会继续攻击没有裂解的单体，与单体中的自由基结合生成小分子气体产物，主要是H。图8是H基攻击各单体的反应路径以及能垒。

图7 苯乙烯产生自由基的反应路径

图8 H基攻击三种单体的反应路径

路径13～路径24经过动力学计算后的反应能垒已经标在图8 中，从图中发现，H 基最容易攻击1,3-丁二烯单体中路径17位置的C—H键，计算结果与杜鸟锋等计算结果相似。同时H基攻击苯乙烯的路径也是较难发生的，与H基攻击苯乙烯的六条路径相比，发现H存在于乙烯基与苯基的结合处时，更容易被夺去。比较路径13～路径24所有的能垒，发现最低的是H·攻击异戊二烯单体生成甲烷。说明热解时气体中甲烷容易生成，徐宗平等、张会亮等、张冰等实验研究结果也表明橡胶热解时，热解气中甲烷含量比较高。

2.1.1 加入Ni催化剂反应路径

2.1.1.1 自由基的生成路径

加入Ni 催化剂对原有自由基生成的路径经过动力学计算后能垒以及催化前能垒见表2，经过计算加入催化剂后各路径（路径25～路径36）的吉布斯自由能变以及焓变随温度变化见图9、图10。

图9 催化后各路径吉布斯自由能变随温度变化

图10 催化后各路径焓变随温度变化

表2 加入Ni催化剂前后自由基生成能垒

加入催化剂后异戊二烯单体产生自由基的路径中，三条路径的能垒都有所降低，能垒降低顺序为∆< ∆< ∆，从表中发现未加入催化剂前路径3 是最容易发生的，表明异戊二烯最易生成·CH自由基，但是加入催化剂后路径2能垒明显降低很多，并接近于路径3的能垒，此路径是生成H基的路径之一，说明Ni 的加入促使异戊二烯催化裂解生成H基。对于1,3-丁二烯单体，能垒降低顺序为∆< ∆< ∆，未加入催化剂前路径5 最容易发生，加入催化剂后路径5的能垒降低最多，使其更容易发生，同时也降低路径4和路径6的能垒，说明Ni 催化剂对1,3-丁二烯生成H 基有明显促进效果。对于苯乙烯单体，能垒降低顺序为∆<∆< ∆< ∆< ∆< ∆，未加入催化剂前路径8最易发生，路径9最难发生，加入催化剂后路径9 的能垒明显降低很多，使原本不容易产H·的路径变得较为容易，同时比较苯环释放H基的路径9～路径11 发现，Ni 的加入会明显降低最不易释放H基的路径9的能垒，说明Ni的加入促使苯乙烯生成H基。

比较异戊二烯、1,3-丁二烯和苯乙烯单体反应的路径，生成H基的九条路径中，路径2能垒降低最为明显，说明Ni 催化剂对异戊二烯单体产生H基效果最为明显，同时从表2中不难发现，苯乙烯单体相较于其他两种单体更不容易释放H基，但是加入Ni 后苯乙烯释放H 基的路径能垒降低程度很大，促进了H 基的生成。未加入催化剂前生成·CH路径3 能垒最低，加入催化剂后生成H 基路径5能垒变得最低，这些说明Ni催化剂对单体生成H基有明显促进作用。图9各路径的吉布斯自由能变的大小顺序与反应能垒相近，且都大于零，说明反应不能自发进行。图10 显示焓变均大于零，说明各路径均为吸热反应。

2.1.1.2 H基攻击单体路径

加入Ni催化剂对原有H基攻击单体的路径经过动力学计算后能垒以及催化前能垒见表3。经过计算，加入催化剂后各路径（路径37～路径48）的吉布斯自由能变以及焓变随温度变化见图11、图12。

图11 催化后各路径吉布斯自由能变随温度变化

图12 催化后各路径焓变随温度变化

加入催化剂后H 基攻击异戊二烯单体生成H的两条路径能垒都降低，尤其是路径14，能垒降低最为明显，同时提高了路径15生成CH的能垒。从表3中发现未加入催化剂前路径15是最容易发生的，但是加入催化剂后路径14 变成最容易发生，说明Ni 的加入促使H生成的同时也抑制了CH的生成。对于1,3-丁二烯单体，加入催化剂前路径17 最容易发生，加入Ni 后路径17 的能垒降低最多，使其更容易发生，说明Ni 催化剂对1,3-丁二烯生成H有明显的促进效果。对于苯乙烯单体，加入催化剂各路径的能垒都有所降低，针对较难发生的路径21 效果最好，同时比较H 基夺取苯环上H 的三条路径，不同位置Ni 的效果不同，对路径21 效果最好，但三条路径能垒都会降低，说明Ni的加入促进苯乙烯生成H。

表3 加入Ni催化剂前后H基攻击单体能垒

比较自由基攻击各单体的反应路径，生成氢气的路径能垒都有所降低，其中路径14 和路径17 能垒明显降低，生成其他气体产物能垒变化不大甚至提高，同时使H基攻击苯乙烯的路径变得较为容易发生，未加入催化剂前路径15 最容易发生，加入催化剂后路径17 变得最为容易，可见Ni 对生成氢气有明显的正向促进作用。图11 各路径的吉布斯自由能变的大小顺序与反应能垒相近，且多数路径都小于零，说明在高温热解阶段，在相应温度下反应能自发进行。图12 显示多数路径焓变小于零，说明各路径均为放热反应。

2.1.2 加入ZSM-5、Ni/ZSM-5催化剂反应路径

2.1.2.1 自由基的生成路径

加入ZSM-5、Ni/ZSM-5 催化剂对原有自由基生成的路径经过动力学计算后能垒以及催化前能垒见表4。从表4 中发现，加入ZSM-5 催化剂后对三种单体释放H基的路径能垒降低并不明显，同时有的路径能垒略有提高。因此对于单体释放H 基来说，ZSM-5 催化剂没有Ni 催化剂效果明显。加入Ni/ZSM-5 催化剂，能垒虽然相较于ZSM-5 降低较多，但是与Ni 相比，还是有一定差距。因此对于单体释放H基来看，催化效果顺序为Ni>Ni/ZSM-5>ZSM-5。值得注意的是路径3，未加入催化剂时异戊二烯单体最容易释放甲基，但是ZSM-5 以及Ni/ZSM-5 的加入使该路径能垒提高很多，说明从抑制异戊二烯单体释放甲基的效果来看，ZSM-5作用比Ni明显。

表4 加入ZSM-5、Ni/ZSM-5催化剂前后自由基生成能垒

2.1.2.2 H基攻击单体路径

加入ZSM-5、Ni/ZSM-5 催化剂对原有H 基攻击单体的路径经过动力学计算后能垒以及催化前能垒见表5。从表5 中发现，加入ZSM-5、Ni/ZSM-5催化剂后，H 基攻击异戊二烯单体，能垒明显降低，但是与Ni 相比都有一定差距。尤其针对H 基攻击苯乙烯路径来看，ZSM-5 以及Ni/ZSM-5 催化制氢效果并不明显。从数据来看，1,3-丁二烯单体是最容易生成H 基以及被H 基攻击生成氢气的，ZSM-5 以及Ni/ZSM-5 的加入也会促进1,3-丁二烯被H基攻击的反应。说明这两种催化剂的加入也会促进制氢反应的进行。

表5 加入ZSM-5、Ni/ZSM-5催化剂前后H基攻击单体能垒

2.2 热解模拟

通过路径分析，已知Ni、ZSM-5、Ni/ZSM-5都有促进生成氢气的可能性，其中Ni 催化剂效果最好，其次是Ni/ZSM-5，最后是ZSM-5，因此选择Ni 催化剂对轮胎橡胶模型进行热解模拟，分析不同比例对热解模拟的影响，观察产物种类及其分布情况。为了探究不同催化比例对热解的影响，选取橡胶与Ni 质量比例分别为5∶1、3∶1、1∶1、1∶0 进行热解模拟，其不同比例的催化剂原子数目分别为18 个、30 个、92 个、0 个。根据各个温度下催化热解情况，热解大致分为两个阶段，低温阶段（500～1500K）和高温阶段（1500～2800K）。

2.2.1 催化热解低温阶段

在低温阶段，一开始温度较低，橡胶模型变化不大，但是随着温度升高，橡胶模型接受能量增多，见表6，1100K 时橡胶热解模型，原子键逐渐拉长，键角发生变化，直到长链中原子的能量升高到可以脱离原子间静电的束缚，发生键的断裂。从表中发现，未加催化剂时变化表现在多数原子键被拉长，加入催化剂后，明显热解较为快一些，部分长链断链，尤其是1∶1 比例，分离出较多单体，说明催化剂的加入对低温阶段有促进作用。当温度升高到1500K时，大量的键会发生断裂，产生多数橡胶单体，少量的H基和含C小分子结构，以及C以上的分子结构。如表6所示，天然橡胶生成的主要单体是异戊二烯，丁苯橡胶生成的主要单体是苯乙烯和1,3-丁二烯，顺丁橡胶生成的主要单体是1,3-丁二烯。1500K 时四个橡胶模型热解形态相似，但是加入催化剂后长链分解出更多单体以及小分子结构，且催化剂比例高者效果明显，比较1∶0和1∶1 两个比例，明显发现1500K 时后者基本不存在较长的长链。这是催化剂加入后的明显变化。

表6 低温阶段热解情况

2.2.2 催化热解高温阶段

随着温度继续升高，热裂解的分子基团发生二次反应，裂解断键形成小分子气体以及小分子自由基。如图13 所示，当温度升高到2800K 时热解基本完成，此时不存在多原子数的结构，热解产物主要以C以下分子结构为主，以及H 自由基、少量未热解的苯环和C以上的分子结构。主要的小分子气体为CH、H、CH、CH，少量存在的C为CH、CH。

图13 2800K未加催化剂热解模型

加入质量比为5∶1 的催化剂后，温度升到2500K时热解基本完成。与未加催化剂相比，热解速率加快，热解产物中小分子产物更多，最为明显的变化是H含量增多，CH以及碳氢化合物减少。这是由于Ni 的加入使生成氢气路径的反应能垒降低，使轮胎橡胶热解更易产生H。说明Ni 作为催化剂对轮胎橡胶模型施加了作用力，改变了橡胶的几何结构，造成了能垒的降低。

加入3∶1 的催化剂后，热解速率比5∶1 比例更快，在2300K时完全裂解。可以看出此时热解产物中小分子产物更多，而且与前两个模型相比H自由基增多，CH以及C～C含量进一步减少。由于催化剂的增多导致热解速率进一步加快，由之前路径模拟结果分析得知，Ni的加入使生成CH的路径能垒增加，因此随着Ni的含量增加，CH含量进一步减少。

在这个阶段，催化剂与轮胎橡胶比例为1∶1时，热解模型形态变化与前几个比例相似，但是热解速率最快，在2100K时基本完全裂解。明显可以看出此时热解产物中小分子气体H产物增多，比5∶1 以及3∶1 比例的更多，同时H 自由基也在增加。CH以及C～C含量进一步减少，虽然产物以C以下结构为主，但是与未加入催化剂相比，H急剧增加。对比四个催化比例热解情况可以发现，Ni添加比例越大，热解终止温度越低，也就说明反应越快，同时对高温阶段产物影响很大。

2.3 综合分析

通过对加入Ni 催化剂的产氢路径计算以及对橡胶催化热解模拟综合来看，Ni 的加入有利于提高氢的产量，且Williams 等使用镍基-镁铝催化剂在固定床进行热解轮胎实验，结果表明，在增加催化剂镍含量、高温度、高催化剂/轮胎比和高蒸汽流量的情况下，可以提高制氢产量。Zhang等采用Ni/AlO等催化剂对废旧轮胎进行热解催化重整实验，结果也表明相同载体情况下，加入Ni 的H产量最高。模拟结果与实验结果基本相符，说明Ni 对于提高热解轮胎产氢率是有帮助的。

3 结论

（1）通过模拟计算出催化热解路径数据，发现加入Ni 催化剂后，单体释放H 基所需能垒都有所下降，对异戊二烯单体释放H基效果最为明显，同时苯乙烯单体相较于其他两种单体更不容易释放H基，但是加入Ni 后苯乙烯释放H 基的路径能垒显著降低。

（2）通过比较H基攻击单体加入催化剂前后能垒数据发现，加入Ni 后，所有生成H路径的能垒都降低，同时提高了生成CH路径的能垒，说明Ni提高了轮胎橡胶热解氢气的产量。

（3）比较加入Ni、ZSM-5、Ni/ZSM-5 催化剂前后能垒发现，它们对生成氢气都有促进作用，其催化效果顺序为Ni>Ni/ZSM-5>ZSM-5。值得注意的是，未加入催化剂前异戊二烯单体较易释放甲基，ZSM-5 以及Ni/ZSM-5 的加入使该路径能垒提高很多，说明从抑制异戊二烯单体释放甲基效果来看，ZSM-5的作用比Ni明显。

（4）对轮胎橡胶加入质量比为5∶1、3∶1、1∶1 以及1∶0（橡胶∶催化剂）的Ni 催化热解过程进行模拟，发现热解分为两个阶段：低温阶段长链释放单体，单体主要是异戊二烯、苯乙烯以及1,3-丁二烯；高温阶段是自由基攻击单体生成小分子结构。加入催化剂后降低热解终止温度，催化剂比例高者终止温度最低。催化剂的加入在低温阶段主要表现为加快热解进程，增加低温阶段时单体数量；在高温阶段主要表现为改变了气体产物分布，Ni的加入使H比例增加，尤其1∶1比例H含量最高，同时降低了C～C比例。