碳基烧蚀材料高温氧化反应机理的分子动力学研究

李磊 刘骁 周述光 何正华 姬广富 葛妮娜

本文采用基于紧束缚密度泛函理论的分子动力学， 研究高温下碳基烧蚀材料三种模型（无缺陷、原子缺陷以及孔洞缺陷）的氧化反应机制. 研究发现高温下的反应产物主要是CO和CO2. CO的产生过程主要源于环氧基团中C-C键的裂解，而CO2的形成则较为复杂，主要源于小分子团簇（C2O2、C3O1、C4O1）的裂解. C-C键裂解是石墨氧化反应的主要途径，C-O键形成是CO和CO2生成速率的控制因素. 此外，体系的温度、缺陷以及孔洞对石墨的氧化反应机制有重要的影响. 通过分析氧化反应速率，计算得到三种模型氧化反应的活化能分别为7.56、2.4和1.6 kcal/mol. 缺陷越明显活化能越低，则氧化反应速率较高，无缺陷的模型由于活化能最高，其氧化反应速率最低.

分子动力学; 氧化反应机制; 反应产物; 氧化反应速率

TB35A2023.034003

收稿日期： 2022-10-24

基金项目： 环境友好能源材料能源材料国家重点实验室自主课题（21fksy27）； 国家数值风洞项目

作者简介： 李磊（1996-）， 男， 江苏泰州人， 硕士研究生， 主要研究材料物性模拟.

通讯作者： 葛妮娜. E-mail： genina911@163.com

Molecular dynamics study on oxidation mechanism of carbon-based ablative materials under high temperature

LI Lei1， LIU Xiao2， ZHOU Shu-Guang2， HE Zheng-Hua3， JI Guang-Fu3， GE Ni-Na1

（1. State Key Laboratory of Environment-Friendly Energy Materials， Southwest University of Science and Technology， Mianyang 621010， China;

2. China Aerodynamics Research and Development Center， Mianyang 621010， China;

3. Institute of Fluid Physics， Chinese Academy of Engineering Physics， Mianyang 621010， China）

Molecular dynamics simulations based on the self-consistent charge density-functional tight-binding （SCC-DFTB） method is used to study the oxidation mechanism of three models （defect-free， atomic-defect and pore-defect） of carbon-based ablative materials under high temperature. It is found that the reaction products at high temperature are mainly CO and CO2. The production process of CO mainly comes from decomposition of C-C bond in epoxy group. The formation of CO2 is relatively complex， mainly from fragmentation of small molecular clusters （C2O2， C3O1， C4O1）. It is found that C-C bond is the main way of graphite oxidation reaction， and C-O bond is the controlling factor of CO and CO2 generation rate. In addition， the temperature， defects and holes of the system have important effects on the oxidation mechanism of graphite. By analyzing the oxidation reaction rate， the calculated activation energies of the three model oxidation reactions are 7.56， 2.4 and 1.6 kcal/mol， respectively. The model with atomic defects and holes corresponds to a low activation energy and a high oxidation reaction rate， while the case without defects associates with the lowest oxidation reaction rate due to the highest activation energy.

Molecular dynamics; Oxidation mechanism; Reaction products; Oxidation reaction rate

1 引 言再入式飞行器进入大气层时，其外表层会与大气层中氧原子发生碰撞，通过侵蚀和氧化导致材料发生物理和化学变化，从而会破坏飞行器的外表层. 因此需要在外表层装备一层烧蚀材料，防止热量传递到飞行器内部，从而达到保护飞行器的目的. 碳基烧蚀材料由于其重量轻、烧蚀速率低和耐高温等优点，通常被用于再入式飞行器的热防护系统，是当前较为理想的高温烧蚀型防热材料［1-11］. 从微观角度看，高温下碳基材料的烧蚀是由高能氧、氮原子碰撞材料引发的氧化还原反应过程. 从微观角度研究碳基烧蚀材料的烧蚀机制对深入认识烧蚀机理和改善飞行器热防护效能具有重要的参考价值.

目前，实验工作经常使用高度有序的热解石墨（HOPG）作为原材料研究高能碰撞反应. Ngo等人［12］根据他们对放置在亚特兰蒂斯号航天飞机上HOPG样品的测量结果，发现每8个撞击氧原子就会蚀刻1个碳原子. Kinoshita等人［13］在室温下将HOPG暴露于5 eV的原子氧后，对其表面进行了表征. 从XPS光谱中发现，HOPG表面每平方厘米内含有约4×1017个氧原子，并通过XPS分析揭示了C-O、C=O和O=C=O物质的存在. 随后，Paci等人［14］通过分子束表面散射实验发现氧原子与HOPG在高温下碰撞时，HOPG表面会先被环氧化合物功能化，随后与半醌发生碰撞导致CO的释放，接着CO与表面上的氧原子反应生成CO2. Kenneth等人［15］在298～423 K的温度下观察到许多嵌入反应区域存在多层圆形凹坑. 接着，Nicholson等人［16］通过一系列实验，研究HOPG暴露于高温氧原子后的形态演变，实验发现当温度为493 K时，侵蚀速率为每入射8个氧原子对应1个碳原子的消失，表明侵蚀速率强烈依赖于HOPG的表面温度. 上述实验数据为理论建模提供了非常重要的参考，在此基础上，Paci等人［17］建立了24碳原子石墨模型进一步研究了碳原子在高温下的去除机制及碳原子的侵蚀速率. Paci等人的理论模型得到了很多团队［18-22］的认可.

高温下HOPG与氧原子反应机理的理论研究还处于初步探索阶段. 杨德军等人［23］通过有限元法分析了热化学烧蚀， 烧蚀表明了退缩及温度场耦合作用下C/C复合材料的烧蚀性能变化规律，求解了进入材料内部的净热流和烧蚀率，表明通过耦合计算可以真实地反应材料的烧蚀特性. 随后，刘骁［24］、国义军等人［25］针对碳化烧蚀热防护材料的三维温度场计算问题进行了研究，清晰地显示了该模型的温度场与密度场分布的多维效应. 李伟等人［26］通过有限体积法并引入固体相体积分数与界面重构技术建立了微观烧蚀数值模型，通过该模型进行数值模拟，模拟发现碳纤维的倾斜角对材料的微观烧蚀行为产生了较大的影响. 在微观方面，首先是基于ReaxFF的分子动力学（MD）的研究. Sriram等人［21］研究了高能原子氧碰撞金刚石表面的蚀刻过程，结果表明在碰撞时金刚石的表面形成各种官能团，与早期的实验和基于第一性原理的计算一致. Poovathingal等人［27］通过扩大模拟体，研究了超热O原子与HOPG的碰撞，模拟观察到碳原子的去除发生在缺陷周围，随着更多的碳原子暴露于氧原子，凹坑会变得更宽. 其次是基于第一性原理的探究.Paci等人［17］通过模拟原子氧与石墨的超热碰撞提出了从石墨中去除碳的机制. 当石墨暴露于高温氧原子时会形成CO和CO2，并且还会形成新的空穴.这些空穴被羰基官能化，在随后的O原子碰撞时再次形成CO/CO2. 此外，Sun等人［28］研究了氧化石墨烯表面的还原反应，确定了石墨和氧化石墨烯表面还原的表面前驱体和详细反应机制，强调了多种富氧前驱体在产生分子CO/CO2气体产物的过程中所起的关键作用.

综上所述，实验上虽可以获取试样的烧蚀氧化形貌，然而这些实验涉及的时间尺度为秒，而大多数化学反应的时间尺度是皮秒或者更小的时间尺度. 因此，实验无法解释这些防热材料样品侵蚀和氧化过程中的反应机制. 与实验相比，理论模拟具有不受实验条件限制的优点，既可以获得实验无法观测到的时间和空间尺度，又可以模拟实际反应过程中所能达到的环境条件，预测和再现反应中的物理与化学现象. 目前微观尺度可以阐述化学反应机制，已有学者通过DFT/ReaxFF来研究石墨的氧化还原反应，但这两种方法都有各自的不足. DFT方法不需要任何经验参数计算精度高，但模拟的体系小，仅为几十个原子. ReaxFF可以模拟较大的模型（可以达到几千个原子），但反应力场方法带有较多的经验参数，对于极端条件下材料性能的预测不够准确. 相比前两种方法，采用基于紧束缚密度泛函理论的DFTB方法，模拟的体系可以达到上百个原子，在精度方面又仅次于DFT.由于在描述能量随电荷密度变化的过程中做了简化处理，计算速度相对较快. 因此，本文采用DFTB的第一性原理对反应过程进行求解验证，针对高温下石墨的氧化还原反应过程进行探究，建立了石墨的三种微观烧蚀模型（无缺陷结构模型、有原子缺陷结构模型以及有孔洞缺陷结构模型），分析不同温度下产物的变化以及反应速率来阐释反应的机理，为碳基烧蚀材料高温抗氧化性能的提升以及碳基材料烧蚀量的预测提供一定的指导.

2 模拟方法

本文使用CP2K代码中实现的基于量子力学的MD方法模拟了石墨的氧化还原反应. 自洽电荷紧束缚密度泛函理论（SCC-DFTB）方法用于计算原子间力，它是基于密度泛函理论（DFT）的总能量相对于电荷密度的二阶展开，能够以自洽的方式描述能量随电荷密度变化的情况. 为了准确描述能量随电荷密度的变化情况，SCC-DFTB所描述的电荷密度是体系能量对电荷密度变化泰勒级数展开的二阶导数.

Etot=∑OCCi〈φiH︿［n0］φi〉-

12∫n0（r）n0（r′）r-r′dr′dr+Exc［n0］-

∫Vxc［n0］（r）n0（r）dr+Eii+12∫［1r-r′-

δ2Excδn（r）δn（r′）n0］×δn（r）δn（r′）drdr′（1）

H︿0=H︿n0是在参考密度下计算的有效的科恩-沙姆哈密顿量，其中n（r）为总电荷密度，n0（r）为参考密度的二阶导数. φi是科恩-沙姆轨道，Exc和Vxc分别为交换关系能和势，Eii为核-核排斥能. 为了推导出SCC-DFTB方法的总能量，对方程（1）中的能量作进一步近似，SCC-DFTB总能量最终可以表达为

Etot=∑iμvCiμCivH0μv+12∑αβγαβΔqαΔqβ+Erep（2）

与密度泛函理论（DFT）相比，唯一不同的是DFT是对空间各处电荷密度进行计算，而DFTB所考虑的电荷密度仅仅是原子的partial charge，使得较小的模拟体系都可以获取到很精确的物理参数. 通过这种计算方法既可以过滤出有用的微观信息，又可以提高计算效率，将计算成本降级几个数量级. 这对于研究分子结构烧蚀材料复杂的物理化学过程等具有实际意义. 与反应力场相比，该方法可以提供更合理可靠的化学特征描述，特别是对于高温高压极端条件的反应，它已在有机系统上成功测试，并准确描述了反应能. 如Sun等人［28］在800 K时计算出在石墨表面的反应具有相应更高的活化能（>1.6 eV）. 此外，该方法还被用于研究氧化石墨烯中碳原子产生CO和CO2的反应. 如Sriram等人［21］通过三组独立的模拟，其计算结果、实验结果以及DFT理论计算结果吻合较好，验证了DFTB方法的可靠性.

3 模型构建

为了研究碳基材料的微观烧蚀机理，本文以高温条件下的石墨为研究对象，探究影响石墨氧化还原反应的因素，基于美国西北大学课题组的计算模拟经验［14，17］，构建了三种石墨模型，如图1所示. 图1a为单层无缺陷石墨模型，此模型是根据Paci等人所建立的由8个氧原子功能化的24原子石墨模型基础上进行扩胞建成的，共包含360个碳原子与120个氧原子. 图1b为单层有原子缺陷石墨模型，此模型是根据Paci等人所建立的有一个原子缺陷的24原子石墨模型基础上进行扩胞建成的，共包含345个碳原子和120个氧原子. 图1c为单层有空位缺陷的石墨模型，根据Poovathingal等人［27］模拟的结果——石墨在高温下反应后会在表面出现孔洞，因此建立了在碳层的内部有孔洞缺陷的模型，以探究初始孔洞是否对反应造成影响.此模型共包含354个碳原子和120个氧原子.

针对上述三种模型进行了结构优化，在几何优化的过程中，每一步最大应力不超过2.3×10-2 eV/，几何变化均方差不高于7.9×10-4 ，所受力的均方差小于1.5×10-2 eV/. 通过优化我们得到了最稳定的结构模型，随后采用NVT系综，模拟温度分别为1500、2000、2500和3000 K，时间步长为0.1 fs，获取的石墨发生氧化反应的轨迹文件，最后使用由Fortran90编译的MD后处理代码逐步分析. 在后处理中，通过分析热点的形成、重要反应产物的形成过程以及化学反应速率，小分子产物的数量随温度和时间的变化关系，阐述影响石墨氧化反应的关键因素.

4 结果分析与讨论

4.1 反应机理的分析

图2为三种模型在不同温度下能量和温度随时间的变化关系. 在0～1 ps时，温度和能量处于下降状态，在这个过程中体系处于不平衡状态；在1～4 ps时，整个反应的能量和温度逐渐处于平衡状态，上下波动的误差不超过5%.从温度的收敛性可以看出，三种模型在不同温度下都处于热平衡状态.

图3和图4描述了三种模型不同温度下CO和CO2的产生过程. 对于单层石墨模型而言，在石墨的氧化还原反应过程中，可以看到有多种团簇和小分子产物产生. 在众多的团簇和小分子产物中，大家关心的主要对象是CO和CO2. Paci等人［17］与魏明刚等人［29］分别从理论和实验的角度出发指出反应过程中的主要反应和反应产物，主要关注的是碳原子的去除机制. Paci所用的是DFT方法，该方法所支持的模型体系较小，故使用了24原子石墨烯结构，因此在反应过程中观察到CO和CO2的过程并不明显，并且在反应过程中CO和CO2产生时间和产生过程无法具体研究. 而本文所使用的DFTB方法可以模拟更大的模型体系，故可以对CO和CO2产生的过程作进一步具体分析.

从图3中可以看到， 三种模型第一次产生CO的时间都在10～14 fs左右，并且随着温度的升高，第一次出现CO的数量也逐渐增多. 另外值得注意的是， CO出现的位置都是在边缘位置，而有孔洞的模型除了在边缘附近会产生CO，在孔洞附近也会优先产生CO，说明石墨的氧化反应大部分是从边缘开始的. 针对中间有孔洞的石墨模型，在孔洞附近也会同时产生CO.该现象与Poovathingal等人［27］的理论结果一致，并且在他们的实验中也观察到了类似现象，证明了石墨在反应的过程中CO主要是从边缘和孔洞附近产生的，随后CO会脱离表面，并在表面留下原型凹坑.

由图4可知，三种模型CO2产生的位置与CO大致相同，都是在石墨的边缘和孔洞附近产生，但形成CO2需要较长的时间. 这是由于CO2的产生过程要较为复杂一些，故第一次所需要形成的时间要晚一些. 研究发现在2500 K时，反应初期有一些小团簇C2O2和C3O1生成，C2O2和C3O1等再次分解生成CO2. 在3000 K时，首先形成比较大的团簇如C10O1，随后大团簇先分解成小团簇如C4O1等，然后才将CO2释放出来. 此外，在分析团簇形成的过程中，较大的团簇出现的时间要比较小团簇出现的时间要晚. 团簇的形成导致在3000和2500 K的温度区间产生的CO2出现的时间滞后.

图5是CO和CO2产生过程示意图. 如图5a所示，一开始由吸附在石墨表面的氧原子发生裂解，随后通过振动激发形成环氧基团. 与羰基相邻的两个环氧基移动到与羰基相邻的碳原子上，随后两个C-C键的连续断裂释放了CO并产生了两个新的羰基，原子缺陷和孔洞缺陷的差异并没有影响到第一次产生CO的过程. 与Paci等人［17］所做的小体系相比，CO的形成过程基本相似，由此可见模型体系的大小并没有影响到CO的形成过程.

CO2的产生过程共有7种产生方式，其中有4种形成过程较为简单. 图5b前4种形成方式都是由一开始吸附在石墨表面的氧原子形成的，而不同的是形成的过程.第一步的主要反应是C-O键和C-C键的裂解，吸附在石墨表面的氧原子裂解成单个原子随后在石墨的表面发生迁移; 第二步的主要反应是新的C-C键和C-O键的生成，迁移后的氧原子与碳原子形成一个CO2分子吸附在石墨表面; 最后就是C-C键裂解释放出CO2. 形成过程较为复杂的有3种， 如图5b后3种所示. 在该方式下CO2是由不同的团簇（C2O2、C3O1和C4O1）裂解开而形成的. 由此可见， CO2的主要来源之一就是团簇. 这些团簇的形成会直接影响CO2的形成.

通过轨迹文件的分析，得出了对应于轨迹文件的成键图，如图5b所示.从图中可以看到C-C键的断裂较容易形成C2O2团簇，随后在该团簇上的C-O键断裂形成CO2. 其次，反应过程中小分子团簇C3O1吸附到石墨表面，随后C-O键断裂一个单独的氧原子从该团簇中释放，最后C-C键断裂也会形成CO2.该形成过程与Paci等人［17］和Sriram等人［21］阐述的CO2的形成过程基本一致. 它们都是由于O原子的碰撞导致形成了环氧化物，最后C-C键断裂释放出CO2. 也就是说这些团簇是形成CO2的主要来源之一.

C4O1的形成主要源于团簇C10O3裂解. 高温下团簇C10O3分解产生两个小的团簇C6O2和C4O1，随后两个团簇结合成一个新的团簇C7O3释放出一个碳原子，最后该碳原子与附近的两个氧原子结合生成了CO2. 该过程与Paci等人［17］所观察到CO2的形成过程基本一致，都是通过大的团簇释放出CO，然后CO与游离在附近的O原子结合形成了CO2.

综上所述，氧原子和石墨相互碰撞形成CO和CO2的过程中，C原子吸附O原子形成的羰基和氧原子的迁移形成环氧化物都起着重要作用. 高温氧原子会导致形成新的空穴.这些空穴被羰基官能化，从而在随后与O原子碰撞时，又成为CO/CO2另一来源. 另一方面，高温原子氧与石墨表面具有高度反应性. 它倾向于在基面上形成环氧基团，一旦形成，这些基团可以迁移，在缺陷处，原子氧可以形成半醌、羰基、内酯、羧酸和醚. 这些官能团是通过实验提出的，这些实验表明，在高温下氧原子在石墨表面发生反应时，C-O、C=O或O-C-O和O-C=O会填充表面，正是这些基团被认为是实验中观察到的CO和CO2的来源. 以上类型的反应也被认为在燃烧化学中发挥着重要作用.

4.2 反应速率的分析

本文通过分析反应过程中的主要化学键随时间的变化情况，进一步研究了石墨的氧化动力学性质［30］. 图6分别显示了三种模型在氧化反应过程中主要化学键的数量随时间的变化情况. 随着反应的进行，C-C键的数量逐渐减少，而C-O键则相应的形成.这些化学键的数量随时间的变化遵循指数形式，并且可以拟合以下指数函数.

C（t）=C0exp（-b·t）（3）

C（t）=C∞［1-exp（-b·t）］（4）

其中，C0和C∞分别代表C-C键和C-O键的数量，b为拟合系数. 众所周知，化学反应速率密切依赖于相应反应物和产物的浓度. 在石墨的氧化反应过程中，C-C键为主要裂解的化学键，C-O键为主要生成的键，而所涉及的化学键的裂变和形成的反应速率一般可以用以下的公式表示.

rfission=k·C（t）n=-dC（t）dt（5）

rformation=k·C∞-C（t）n=dC（t）dt（6）

其中，r是反应速率，k是有效反应速率常数，n是反应级数. 将公式（5）（6）转化为对数形式（见公式（7）（8）），得到相应的线性函数：ln（-dC（t）/dt）～ln（C（t））h和ln （dC（t）/dt）～ln（C∞- C（t））. 它们的斜率和截距是分别对应n和ln（k）.将公式（3）（4）代入公式（7）（8），发现对于反应物的数量（裂变反应）和产物（生成反应），化学键的裂变和形成是一级反应（n=1），速率常数k就是拟合系数b，可以直接通过化学键的数量演化的指数拟合得到.

ln-dC（t）dt=lnk+n·ln （C（t））（7）

lndC（t）dt=ln（k）+n·ln （C∞-C（t））（8）

图6a～6c显示了三种模型在氧化反应过程中C-C键的裂变和C-O键的生成的指数拟合. 考虑到在1500 K时反应C-C键的裂变和C-O键的生成较为不明显，因此只研究了2000、2500和3000 K时的数量变化. 我们可以看出拟合曲线与数据的一致性，说明拟合合理可靠. 为了进一步揭示动力学特性，本文使用经典的Arrhenius形式（k=A*exp（-Ea/RT））来拟合这些数据，其中k是有效的速率常数，A是频率因子，Ea是活化能，R是气体常数，T是温度（K）. C-C键和C-O键的反应速率常数的温度依赖性由Arrhenius公式的对数形式ln（k）=ln（A）-Ea/RT表示，如图7a和7b所示. 三种模型在C-C键的裂解上计算得到的反应活化能分别为7.56、2.4和1.6 kcal/mol. 从图7a中可以看出， 单层无缺陷模型具有最高的反应活化能，而另外两种模型由于一部分氧原子代替了碳原子和一部分碳原子的缺失，导致活化能降低. 这说明无缺陷模型的C-C键不易裂解，因此具有最低的反应速率. 三种模型在C-O键的生成上计算得到的反应活化能分别为5.15、3.57和4.82 kcal/mol. 从图7b中可以看出， 单层无缺陷模型也具有最高的活化能，而另外两种模型由于碳原子的缺失导致活化能降低，使反应速率加快，说明C-O键的形成是CO和CO2生成速率的主要的控制因素. 在这样的温度条件下，石墨在分解过程中会留下许多C-C键. 这就是在反应过程中产生许多碳原子团簇的原因. 此外，值得注意的是，C-C键的裂解是氧化反应的主要途径，而C-O键的生成是CO和CO2生成的主要控制因素，将显著影响碳的氧化形成最终的气态氧化物. 这与前面在反应机理中观察到的现象一致，进一步验证了C-C键的裂解是氧化反应主要反应途径，C-O键是CO和CO2生成速率的主要控制因素.

4.3 小分子产物的分析

实验上显示石墨在氧化还原反应过程的主要产物是CO和CO2.图8所示为单层无缺陷石墨模型CO和CO2数量随时间的变化. 从图8a中可以看出， CO的数量呈现一定的规律性. 当温度为3000 K时，CO数量最多; 但当温度为1500 K时，只有少量的CO产生. 从CO产生的数量上可以看出，在高温下石墨的氧化还原反应是以CO为主导的. 当温度越高，产生的CO数量越多，反之则越少. 该模型在反映中出现的CO和CO2的数量上的趋势， 与国义军等人［31， 32］在碳基烧蚀材料的快慢反应理论和烧蚀机理中计算的趋势相吻合. 即， 在1700 K以上，温度越高生成物主要是CO; 1500 K以下，温度越低生成物主要是CO2. 从图8b中可以看出， 在低温时只有在反应后期才会产生1～2个CO2; 当温度升高，CO2会逐渐消失. 在高温时，反应初始阶段会产生1～2个CO2，在反应的后期阶段就不再会产生CO2. 由此可以看出， CO2在低温下会较早产生，随着反应的进行则不会再产生; 而在高温时CO产生的数量较多，占主导地位.

根据上述单层无缺陷的石墨模型的分析，我们还做了单层有原子缺陷（图8b）和孔洞缺陷（图8c）的石墨模型进行比较. 从图中可以看到，单层有缺陷（图8b）的石墨模型和无缺陷（图8a）的石墨模型在分解产物中存在一定的差异. 单层原子缺陷的石墨在氧化反应过程中产生的CO增长较快，由于石墨表面的一部分碳原子被氧原子所替代，加快了氧化还原反应，故CO增长的速度较快. 所以缺陷对石墨模型的影响是不能忽略的. 实际材料中缺陷分为原子缺陷和空位缺陷. 图8c为单层有孔洞的石墨模型CO和CO2数量随时间的变化情况. 与上述的单层有原子缺陷的石墨模型相比，空位缺陷所产生的CO数量更多，并且空位缺陷所产生CO数量的趋势比原子缺陷所产生的较为稳定. 在CO的产生上，两者都是在低温下产生的数量较为稳定，在高温下都呈现缓慢增长的状态. 而CO2的产生两种模型存在差异，即是在反应的初期，有空位缺陷的模型会产生CO2，有原子缺陷的模型没有产生CO2. 因此缺陷的种类也会影响石墨氧化还原反应产生的CO和CO2数量.

5 结 论

本文采用基于紧束缚密度泛函理论分子动力学方法研究高温下三种石墨模型（无缺陷、原子缺陷以及孔洞缺陷）的氧化还原反应. 研究发现三种石墨模型的共同点是， 在高温下CO和CO2是石墨在氧化还原反应中的主要产物.CO的主要来源是由环氧基团上的C-C键断裂后，碳原子与氧原子结合形成CO;CO2的形成较为复杂，主要来源是小分子团簇（C2O2、C3O1和C4O1）的裂解. 通过分析氧化反应速率，C-C键断裂是石墨烧蚀反应的主要途径，C-O键的断裂是CO和CO2生成速率的主要控制因素. 此外，三种模型氧化反应差异是，在高温下无缺陷的石墨模型产生CO和CO2的数量较多，而有缺陷和有孔洞的模型产生数量较少，表明缺陷和孔洞在高温下会影响CO和CO2产生的数量. 在氧化反应速率方面，三种模型计算得到的反应活化能分别为7.56、2.4和1.6 kcal/mol，表明无缺陷的石墨模型活化能高，表现出最低的烧蚀反应速率，而原子缺陷和孔洞缺陷则会加快烧蚀反应速率. 本文的研究对烧蚀材料的性能改进具有参考意义.

参考文献：

［1］ 邢亚娟， 孙波， 高坤， 等. 航天飞行器热防护系统及防热材料研究现状［J］. 力学学报， 2011， 43： 289.

［2］ 胡连成. 前苏联轻质防热材料和碳/酚醛复合材料［J］. 宇航材料工艺， 2018， 48： 9.

［3］ 靳亮， 王德， 王耀辉， 等. 热防护材料烧蚀性能测试方法研究［J］. 固体火箭技术， 2018， 41： 245.

［4］ 吴晓宏， 陆小龙， 李涛， 等. 轻质烧蚀材料研究综述［J］. 航天器环境工程， 2011， 28： 313.

［5］ 高守臻， 魏化震， 李大勇， 等. 烧蚀材料综述［J］. 化工新型材料， 2009， 37： 19.

［6］ 张威， 曾明， 肖凌飞， 等. 碳-酚醛材料烧蚀热解对再入流场特性影响的数值计算［J］. 国防科技大学学报， 2014， 36： 41.

［7］ 高俊杰， 俞继军， 韩海涛， 等. 树脂基烧蚀材料细观传热特性预测［J］. 航空学报， 2017， 38： 170.

［8］ 王子龙， 邢素丽， 尹昌平， 等. 高残炭耐烧蚀树脂基体研究进展［J］. 工程塑料应用， 2018， 46： 131.

［9］ 贺广东， 邹宁宇. 对高硅氧纤维耐烧蚀材料的几点探讨［J］. 玻璃纤维， 2015（6）： 40.

［10］ 霍肖旭， 曾晓梅， 刘红林. 炭/炭复合材料及炭纤维增强塑料在固体火箭喷管上的应用［J］. 炭素技术， 2001（3）： 23.

［11］ 李林杰， 方国东， 易法军， 等. 高硅氧/酚醛复合材料热变形实验测试及表面烧蚀形貌分析［J］. 固体火箭技术， 2015， 38： 445.

［12］ Ngo T， Snyder E J， Tong W M， et al. O atom etching of graphite in low earth orbit ［J］. Surf Sci， 1994， 314： 817.

［13］ Kinoshit H， Umeno M， Tagawa M， et al. Hyperthermal atomic oxygen beam-induced etching of HOPG（0001） studied by X-ray photoelectron spectroscopy and scanning tunneling microscopy ［J］. Surf Sci， 1999， 440： 49.

［14］ Paci J T， Upadhyaya H P， Zhang J M， et al. Theoretical and experimental studies of the reactions between hyperthermal O（3P） and graphite： graphene-based direct dynamics and beam-surface scattering approaches ［J］. J Phys Chem A， 2009， 113： 4677.

［15］ Kenneth T N. Nucleation andgrowth of nanoscale to microscale cylindrical pits in highly-ordered pyrolytic graphite upon hyperthermal atomic oxygen exposure ［J］. High Perform Polym， 2004， 16： 197.

［16］ Nicholson K T， Minton T K， Sibener S J. Spatially anisotropic etching of graphite by hyperthermal atomic oxygen ［J］. J Phys Chem B， 2005， 109： 8476.

［17］ Paci J T， Minton T K， Schatz G C. Hyperthermal oxidation of graphite and diamond ［J］. Acc Chem Res， 2012， 45： 1973.

［18］ Hariharan S， Majumder M， Edel R， et al. Exploratory direct dynamics simulations of O-3（2） reaction with graphene at high temperatures ［J］. J Phys Chem C， 2019， 122： 51.

［19］ Yi M， Shen Z G， Zhao X H， et al. Exploring few-layer graphene and graphene oxide as fillers to enhance the oxygen-atom corrosion resistance of composites ［J］. Phys Chem Chem Phys， 2014， 16： 23.

［20］ Meng D Z， Yan G， Yue W， et al. Thermal damage mechanisms of Si-coated diamond powder based polycrystalline diamond ［J］. J Eur Ceram Soc， 2018， 38： 13.

［21］ Sriram G S， Adri C T. Direction dependent etching of diamond surfaces by hyperthermal atomic oxygen： a ReaxFF based molecular dynamics study ［J］. Carbon， 2015， 82： 314.

［22］ Srinivasan S G， Adri C T. Molecular-dynamics-based study of the collisions of hyperthermal atomic oxygen with graphene using the ReaxFF reactive force field ［J］. J Phys Chem A， 2011， 115： 13269.

［23］ 杨德军. 防热复合材料烧蚀行为的数值模拟［D］. 兰州： 兰州理工大学， 2013.

［24］ 刘骁. 烧蚀热防护系统三维热响应有限元计算方法研究［D］. 绵阳：中国空气动力研究与发展中心， 2016.

［25］ 国义军， 石卫波， 曾慕， 等．高超声速飞行器烧蚀防热理论与应用［M］． 北京： 科学出版社， 2019.

［26］ 李伟， 方国东， 李玮洁， 等. 碳纤维增强复合材料微观烧蚀行为数值模拟［J］. 力学学报， 2019， 51： 835.

［27］ Poovathingal S， Schwartzentruber T E， Srinivasan S G， et al. Large scale computational chemistry modeling of the oxidation of highly oriented pyrolytic graphite ［J］. J Phys Chem A， 2013， 117： 2692.

［28］ Sun T， Stefano B， Stefano F. Surface precursors and reaction mechanisms for the thermal reduction of graphene basal surfaces oxidized by atomic oxygen ［J］. J Phys Chem C， 2011， 115： 4730.

［29］ 魏明刚， 罗恒军， 张海成， 等. 300M钢的脱碳行为演化及防护研究［J］. 四川大学学报： 自然科学版， 2021， 58： 143.

［30］ He Z H， Yu Y， Huang YY， et al. Reaction kinetic properties of 1， 3， 5-triamino-2， 4， 6-trinitrobenzene： a DFTB study of thermal decomposition ［J］. New J Chem， 2019， 43： 18027.

［31］ 国义军， 周述光， 曾磊， 等. 碳基材料氧化烧蚀的快/慢反应动力学问题［J］. 实验流体力学， 2019， 33： 17.

［32］ 国义军， 桂业伟， 童福林， 等. 碳/碳化锆复合材料烧蚀机理和计算方法研究［J］. 空气动力学报， 2012， 31： 22.

引用本文格式：

中 文： 李磊， 刘骁， 周述光， 等. 碳基烧蚀材料高温氧化反应机理的分子动力学研究［J］. 四川大学学报： 自然科学版， 2023， 60： 034003.

英 文： Li L， Liu X， Zhou S G， et al. Molecular dynamics study on oxidation mechanism of carbon-based ablative materials under high temperature ［J］. J Sichuan Univ： Nat Sci Ed， 2023， 60： 034003.