李磊 劉驍 周述光 何正華 姬廣富 葛妮娜

本文采用基于緊束縛密度泛函理論的分子動(dòng)力學(xué)， 研究高溫下碳基燒蝕材料三種模型（無(wú)缺陷、原子缺陷以及孔洞缺陷）的氧化反應(yīng)機(jī)制. 研究發(fā)現(xiàn)高溫下的反應(yīng)產(chǎn)物主要是CO和CO2. CO的產(chǎn)生過(guò)程主要源于環(huán)氧基團(tuán)中C-C鍵的裂解，而CO2的形成則較為復(fù)雜，主要源于小分子團(tuán)簇（C2O2、C3O1、C4O1）的裂解. C-C鍵裂解是石墨氧化反應(yīng)的主要途徑，C-O鍵形成是CO和CO2生成速率的控制因素. 此外，體系的溫度、缺陷以及孔洞對(duì)石墨的氧化反應(yīng)機(jī)制有重要的影響. 通過(guò)分析氧化反應(yīng)速率，計(jì)算得到三種模型氧化反應(yīng)的活化能分別為7.56、2.4和1.6 kcal/mol. 缺陷越明顯活化能越低，則氧化反應(yīng)速率較高，無(wú)缺陷的模型由于活化能最高，其氧化反應(yīng)速率最低.

分子動(dòng)力學(xué); 氧化反應(yīng)機(jī)制; 反應(yīng)產(chǎn)物; 氧化反應(yīng)速率

TB35A2023.034003

收稿日期： 2022-10-24

基金項(xiàng)目： 環(huán)境友好能源材料能源材料國(guó)家重點(diǎn)實(shí)驗(yàn)室自主課題（21fksy27）； 國(guó)家數(shù)值風(fēng)洞項(xiàng)目

作者簡(jiǎn)介： 李磊（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 引 言再入式飛行器進(jìn)入大氣層時(shí)，其外表層會(huì)與大氣層中氧原子發(fā)生碰撞，通過(guò)侵蝕和氧化導(dǎo)致材料發(fā)生物理和化學(xué)變化，從而會(huì)破壞飛行器的外表層. 因此需要在外表層裝備一層燒蝕材料，防止熱量傳遞到飛行器內(nèi)部，從而達(dá)到保護(hù)飛行器的目的. 碳基燒蝕材料由于其重量輕、燒蝕速率低和耐高溫等優(yōu)點(diǎn)，通常被用于再入式飛行器的熱防護(hù)系統(tǒng)，是當(dāng)前較為理想的高溫?zé)g型防熱材料［1-11］. 從微觀角度看，高溫下碳基材料的燒蝕是由高能氧、氮原子碰撞材料引發(fā)的氧化還原反應(yīng)過(guò)程. 從微觀角度研究碳基燒蝕材料的燒蝕機(jī)制對(duì)深入認(rèn)識(shí)燒蝕機(jī)理和改善飛行器熱防護(hù)效能具有重要的參考價(jià)值.

目前，實(shí)驗(yàn)工作經(jīng)常使用高度有序的熱解石墨（HOPG）作為原材料研究高能碰撞反應(yīng). Ngo等人［12］根據(jù)他們對(duì)放置在亞特蘭蒂斯號(hào)航天飛機(jī)上HOPG樣品的測(cè)量結(jié)果，發(fā)現(xiàn)每8個(gè)撞擊氧原子就會(huì)蝕刻1個(gè)碳原子. Kinoshita等人［13］在室溫下將HOPG暴露于5 eV的原子氧后，對(duì)其表面進(jìn)行了表征. 從XPS光譜中發(fā)現(xiàn)，HOPG表面每平方厘米內(nèi)含有約4×1017個(gè)氧原子，并通過(guò)XPS分析揭示了C-O、C=O和O=C=O物質(zhì)的存在. 隨后，Paci等人［14］通過(guò)分子束表面散射實(shí)驗(yàn)發(fā)現(xiàn)氧原子與HOPG在高溫下碰撞時(shí)，HOPG表面會(huì)先被環(huán)氧化合物功能化，隨后與半醌發(fā)生碰撞導(dǎo)致CO的釋放，接著CO與表面上的氧原子反應(yīng)生成CO2. Kenneth等人［15］在298～423 K的溫度下觀察到許多嵌入反應(yīng)區(qū)域存在多層圓形凹坑. 接著，Nicholson等人［16］通過(guò)一系列實(shí)驗(yàn)，研究HOPG暴露于高溫氧原子后的形態(tài)演變，實(shí)驗(yàn)發(fā)現(xiàn)當(dāng)溫度為493 K時(shí)，侵蝕速率為每入射8個(gè)氧原子對(duì)應(yīng)1個(gè)碳原子的消失，表明侵蝕速率強(qiáng)烈依賴于HOPG的表面溫度. 上述實(shí)驗(yàn)數(shù)據(jù)為理論建模提供了非常重要的參考，在此基礎(chǔ)上，Paci等人［17］建立了24碳原子石墨模型進(jìn)一步研究了碳原子在高溫下的去除機(jī)制及碳原子的侵蝕速率. Paci等人的理論模型得到了很多團(tuán)隊(duì)［18-22］的認(rèn)可.

高溫下HOPG與氧原子反應(yīng)機(jī)理的理論研究還處于初步探索階段. 楊德軍等人［23］通過(guò)有限元法分析了熱化學(xué)燒蝕， 燒蝕表明了退縮及溫度場(chǎng)耦合作用下C/C復(fù)合材料的燒蝕性能變化規(guī)律，求解了進(jìn)入材料內(nèi)部的凈熱流和燒蝕率，表明通過(guò)耦合計(jì)算可以真實(shí)地反應(yīng)材料的燒蝕特性. 隨后，劉驍［24］、國(guó)義軍等人［25］針對(duì)碳化燒蝕熱防護(hù)材料的三維溫度場(chǎng)計(jì)算問(wèn)題進(jìn)行了研究，清晰地顯示了該模型的溫度場(chǎng)與密度場(chǎng)分布的多維效應(yīng). 李偉等人［26］通過(guò)有限體積法并引入固體相體積分?jǐn)?shù)與界面重構(gòu)技術(shù)建立了微觀燒蝕數(shù)值模型，通過(guò)該模型進(jìn)行數(shù)值模擬，模擬發(fā)現(xiàn)碳纖維的傾斜角對(duì)材料的微觀燒蝕行為產(chǎn)生了較大的影響. 在微觀方面，首先是基于ReaxFF的分子動(dòng)力學(xué)（MD）的研究. Sriram等人［21］研究了高能原子氧碰撞金剛石表面的蝕刻過(guò)程，結(jié)果表明在碰撞時(shí)金剛石的表面形成各種官能團(tuán)，與早期的實(shí)驗(yàn)和基于第一性原理的計(jì)算一致. Poovathingal等人［27］通過(guò)擴(kuò)大模擬體，研究了超熱O原子與HOPG的碰撞，模擬觀察到碳原子的去除發(fā)生在缺陷周?chē)?，隨著更多的碳原子暴露于氧原子，凹坑會(huì)變得更寬. 其次是基于第一性原理的探究.Paci等人［17］通過(guò)模擬原子氧與石墨的超熱碰撞提出了從石墨中去除碳的機(jī)制. 當(dāng)石墨暴露于高溫氧原子時(shí)會(huì)形成CO和CO2，并且還會(huì)形成新的空穴.這些空穴被羰基官能化，在隨后的O原子碰撞時(shí)再次形成CO/CO2. 此外，Sun等人［28］研究了氧化石墨烯表面的還原反應(yīng)，確定了石墨和氧化石墨烯表面還原的表面前驅(qū)體和詳細(xì)反應(yīng)機(jī)制，強(qiáng)調(diào)了多種富氧前驅(qū)體在產(chǎn)生分子CO/CO2氣體產(chǎn)物的過(guò)程中所起的關(guān)鍵作用.

綜上所述，實(shí)驗(yàn)上雖可以獲取試樣的燒蝕氧化形貌，然而這些實(shí)驗(yàn)涉及的時(shí)間尺度為秒，而大多數(shù)化學(xué)反應(yīng)的時(shí)間尺度是皮秒或者更小的時(shí)間尺度. 因此，實(shí)驗(yàn)無(wú)法解釋這些防熱材料樣品侵蝕和氧化過(guò)程中的反應(yīng)機(jī)制. 與實(shí)驗(yàn)相比，理論模擬具有不受實(shí)驗(yàn)條件限制的優(yōu)點(diǎn)，既可以獲得實(shí)驗(yàn)無(wú)法觀測(cè)到的時(shí)間和空間尺度，又可以模擬實(shí)際反應(yīng)過(guò)程中所能達(dá)到的環(huán)境條件，預(yù)測(cè)和再現(xiàn)反應(yīng)中的物理與化學(xué)現(xiàn)象. 目前微觀尺度可以闡述化學(xué)反應(yīng)機(jī)制，已有學(xué)者通過(guò)DFT/ReaxFF來(lái)研究石墨的氧化還原反應(yīng)，但這兩種方法都有各自的不足. DFT方法不需要任何經(jīng)驗(yàn)參數(shù)計(jì)算精度高，但模擬的體系小，僅為幾十個(gè)原子. ReaxFF可以模擬較大的模型（可以達(dá)到幾千個(gè)原子），但反應(yīng)力場(chǎng)方法帶有較多的經(jīng)驗(yàn)參數(shù)，對(duì)于極端條件下材料性能的預(yù)測(cè)不夠準(zhǔn)確. 相比前兩種方法，采用基于緊束縛密度泛函理論的DFTB方法，模擬的體系可以達(dá)到上百個(gè)原子，在精度方面又僅次于DFT.由于在描述能量隨電荷密度變化的過(guò)程中做了簡(jiǎn)化處理，計(jì)算速度相對(duì)較快. 因此，本文采用DFTB的第一性原理對(duì)反應(yīng)過(guò)程進(jìn)行求解驗(yàn)證，針對(duì)高溫下石墨的氧化還原反應(yīng)過(guò)程進(jìn)行探究，建立了石墨的三種微觀燒蝕模型（無(wú)缺陷結(jié)構(gòu)模型、有原子缺陷結(jié)構(gòu)模型以及有孔洞缺陷結(jié)構(gòu)模型），分析不同溫度下產(chǎn)物的變化以及反應(yīng)速率來(lái)闡釋反應(yīng)的機(jī)理，為碳基燒蝕材料高溫抗氧化性能的提升以及碳基材料燒蝕量的預(yù)測(cè)提供一定的指導(dǎo).

2 模擬方法

本文使用CP2K代碼中實(shí)現(xiàn)的基于量子力學(xué)的MD方法模擬了石墨的氧化還原反應(yīng). 自洽電荷緊束縛密度泛函理論（SCC-DFTB）方法用于計(jì)算原子間力，它是基于密度泛函理論（DFT）的總能量相對(duì)于電荷密度的二階展開(kāi)，能夠以自洽的方式描述能量隨電荷密度變化的情況. 為了準(zhǔn)確描述能量隨電荷密度的變化情況，SCC-DFTB所描述的電荷密度是體系能量對(duì)電荷密度變化泰勒級(jí)數(shù)展開(kāi)的二階導(dǎo)數(shù).

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是在參考密度下計(jì)算的有效的科恩-沙姆哈密頓量，其中n（r）為總電荷密度，n0（r）為參考密度的二階導(dǎo)數(shù). φi是科恩-沙姆軌道，Exc和Vxc分別為交換關(guān)系能和勢(shì)，Eii為核-核排斥能. 為了推導(dǎo)出SCC-DFTB方法的總能量，對(duì)方程（1）中的能量作進(jìn)一步近似，SCC-DFTB總能量最終可以表達(dá)為

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

與密度泛函理論（DFT）相比，唯一不同的是DFT是對(duì)空間各處電荷密度進(jìn)行計(jì)算，而DFTB所考慮的電荷密度僅僅是原子的partial charge，使得較小的模擬體系都可以獲取到很精確的物理參數(shù). 通過(guò)這種計(jì)算方法既可以過(guò)濾出有用的微觀信息，又可以提高計(jì)算效率，將計(jì)算成本降級(jí)幾個(gè)數(shù)量級(jí). 這對(duì)于研究分子結(jié)構(gòu)燒蝕材料復(fù)雜的物理化學(xué)過(guò)程等具有實(shí)際意義. 與反應(yīng)力場(chǎng)相比，該方法可以提供更合理可靠的化學(xué)特征描述，特別是對(duì)于高溫高壓極端條件的反應(yīng)，它已在有機(jī)系統(tǒng)上成功測(cè)試，并準(zhǔn)確描述了反應(yīng)能. 如Sun等人［28］在800 K時(shí)計(jì)算出在石墨表面的反應(yīng)具有相應(yīng)更高的活化能（>1.6 eV）. 此外，該方法還被用于研究氧化石墨烯中碳原子產(chǎn)生CO和CO2的反應(yīng). 如Sriram等人［21］通過(guò)三組獨(dú)立的模擬，其計(jì)算結(jié)果、實(shí)驗(yàn)結(jié)果以及DFT理論計(jì)算結(jié)果吻合較好，驗(yàn)證了DFTB方法的可靠性.

3 模型構(gòu)建

為了研究碳基材料的微觀燒蝕機(jī)理，本文以高溫條件下的石墨為研究對(duì)象，探究影響石墨氧化還原反應(yīng)的因素，基于美國(guó)西北大學(xué)課題組的計(jì)算模擬經(jīng)驗(yàn)［14，17］，構(gòu)建了三種石墨模型，如圖1所示. 圖1a為單層無(wú)缺陷石墨模型，此模型是根據(jù)Paci等人所建立的由8個(gè)氧原子功能化的24原子石墨模型基礎(chǔ)上進(jìn)行擴(kuò)胞建成的，共包含360個(gè)碳原子與120個(gè)氧原子. 圖1b為單層有原子缺陷石墨模型，此模型是根據(jù)Paci等人所建立的有一個(gè)原子缺陷的24原子石墨模型基礎(chǔ)上進(jìn)行擴(kuò)胞建成的，共包含345個(gè)碳原子和120個(gè)氧原子. 圖1c為單層有空位缺陷的石墨模型，根據(jù)Poovathingal等人［27］模擬的結(jié)果——石墨在高溫下反應(yīng)后會(huì)在表面出現(xiàn)孔洞，因此建立了在碳層的內(nèi)部有孔洞缺陷的模型，以探究初始孔洞是否對(duì)反應(yīng)造成影響.此模型共包含354個(gè)碳原子和120個(gè)氧原子.

針對(duì)上述三種模型進(jìn)行了結(jié)構(gòu)優(yōu)化，在幾何優(yōu)化的過(guò)程中，每一步最大應(yīng)力不超過(guò)2.3×10-2 eV/，幾何變化均方差不高于7.9×10-4 ，所受力的均方差小于1.5×10-2 eV/. 通過(guò)優(yōu)化我們得到了最穩(wěn)定的結(jié)構(gòu)模型，隨后采用NVT系綜，模擬溫度分別為1500、2000、2500和3000 K，時(shí)間步長(zhǎng)為0.1 fs，獲取的石墨發(fā)生氧化反應(yīng)的軌跡文件，最后使用由Fortran90編譯的MD后處理代碼逐步分析. 在后處理中，通過(guò)分析熱點(diǎn)的形成、重要反應(yīng)產(chǎn)物的形成過(guò)程以及化學(xué)反應(yīng)速率，小分子產(chǎn)物的數(shù)量隨溫度和時(shí)間的變化關(guān)系，闡述影響石墨氧化反應(yīng)的關(guān)鍵因素.

4 結(jié)果分析與討論

4.1 反應(yīng)機(jī)理的分析

圖2為三種模型在不同溫度下能量和溫度隨時(shí)間的變化關(guān)系. 在0～1 ps時(shí)，溫度和能量處于下降狀態(tài)，在這個(gè)過(guò)程中體系處于不平衡狀態(tài)；在1～4 ps時(shí)，整個(gè)反應(yīng)的能量和溫度逐漸處于平衡狀態(tài)，上下波動(dòng)的誤差不超過(guò)5%.從溫度的收斂性可以看出，三種模型在不同溫度下都處于熱平衡狀態(tài).

圖3和圖4描述了三種模型不同溫度下CO和CO2的產(chǎn)生過(guò)程. 對(duì)于單層石墨模型而言，在石墨的氧化還原反應(yīng)過(guò)程中，可以看到有多種團(tuán)簇和小分子產(chǎn)物產(chǎn)生. 在眾多的團(tuán)簇和小分子產(chǎn)物中，大家關(guān)心的主要對(duì)象是CO和CO2. Paci等人［17］與魏明剛等人［29］分別從理論和實(shí)驗(yàn)的角度出發(fā)指出反應(yīng)過(guò)程中的主要反應(yīng)和反應(yīng)產(chǎn)物，主要關(guān)注的是碳原子的去除機(jī)制. Paci所用的是DFT方法，該方法所支持的模型體系較小，故使用了24原子石墨烯結(jié)構(gòu)，因此在反應(yīng)過(guò)程中觀察到CO和CO2的過(guò)程并不明顯，并且在反應(yīng)過(guò)程中CO和CO2產(chǎn)生時(shí)間和產(chǎn)生過(guò)程無(wú)法具體研究. 而本文所使用的DFTB方法可以模擬更大的模型體系，故可以對(duì)CO和CO2產(chǎn)生的過(guò)程作進(jìn)一步具體分析.

從圖3中可以看到， 三種模型第一次產(chǎn)生CO的時(shí)間都在10～14 fs左右，并且隨著溫度的升高，第一次出現(xiàn)CO的數(shù)量也逐漸增多. 另外值得注意的是， CO出現(xiàn)的位置都是在邊緣位置，而有孔洞的模型除了在邊緣附近會(huì)產(chǎn)生CO，在孔洞附近也會(huì)優(yōu)先產(chǎn)生CO，說(shuō)明石墨的氧化反應(yīng)大部分是從邊緣開(kāi)始的. 針對(duì)中間有孔洞的石墨模型，在孔洞附近也會(huì)同時(shí)產(chǎn)生CO.該現(xiàn)象與Poovathingal等人［27］的理論結(jié)果一致，并且在他們的實(shí)驗(yàn)中也觀察到了類(lèi)似現(xiàn)象，證明了石墨在反應(yīng)的過(guò)程中CO主要是從邊緣和孔洞附近產(chǎn)生的，隨后CO會(huì)脫離表面，并在表面留下原型凹坑.

由圖4可知，三種模型CO2產(chǎn)生的位置與CO大致相同，都是在石墨的邊緣和孔洞附近產(chǎn)生，但形成CO2需要較長(zhǎng)的時(shí)間. 這是由于CO2的產(chǎn)生過(guò)程要較為復(fù)雜一些，故第一次所需要形成的時(shí)間要晚一些. 研究發(fā)現(xiàn)在2500 K時(shí)，反應(yīng)初期有一些小團(tuán)簇C2O2和C3O1生成，C2O2和C3O1等再次分解生成CO2. 在3000 K時(shí)，首先形成比較大的團(tuán)簇如C10O1，隨后大團(tuán)簇先分解成小團(tuán)簇如C4O1等，然后才將CO2釋放出來(lái). 此外，在分析團(tuán)簇形成的過(guò)程中，較大的團(tuán)簇出現(xiàn)的時(shí)間要比較小團(tuán)簇出現(xiàn)的時(shí)間要晚. 團(tuán)簇的形成導(dǎo)致在3000和2500 K的溫度區(qū)間產(chǎn)生的CO2出現(xiàn)的時(shí)間滯后.

圖5是CO和CO2產(chǎn)生過(guò)程示意圖. 如圖5a所示，一開(kāi)始由吸附在石墨表面的氧原子發(fā)生裂解，隨后通過(guò)振動(dòng)激發(fā)形成環(huán)氧基團(tuán). 與羰基相鄰的兩個(gè)環(huán)氧基移動(dòng)到與羰基相鄰的碳原子上，隨后兩個(gè)C-C鍵的連續(xù)斷裂釋放了CO并產(chǎn)生了兩個(gè)新的羰基，原子缺陷和孔洞缺陷的差異并沒(méi)有影響到第一次產(chǎn)生CO的過(guò)程. 與Paci等人［17］所做的小體系相比，CO的形成過(guò)程基本相似，由此可見(jiàn)模型體系的大小并沒(méi)有影響到CO的形成過(guò)程.

CO2的產(chǎn)生過(guò)程共有7種產(chǎn)生方式，其中有4種形成過(guò)程較為簡(jiǎn)單. 圖5b前4種形成方式都是由一開(kāi)始吸附在石墨表面的氧原子形成的，而不同的是形成的過(guò)程.第一步的主要反應(yīng)是C-O鍵和C-C鍵的裂解，吸附在石墨表面的氧原子裂解成單個(gè)原子隨后在石墨的表面發(fā)生遷移; 第二步的主要反應(yīng)是新的C-C鍵和C-O鍵的生成，遷移后的氧原子與碳原子形成一個(gè)CO2分子吸附在石墨表面; 最后就是C-C鍵裂解釋放出CO2. 形成過(guò)程較為復(fù)雜的有3種， 如圖5b后3種所示. 在該方式下CO2是由不同的團(tuán)簇（C2O2、C3O1和C4O1）裂解開(kāi)而形成的. 由此可見(jiàn)， CO2的主要來(lái)源之一就是團(tuán)簇. 這些團(tuán)簇的形成會(huì)直接影響CO2的形成.

通過(guò)軌跡文件的分析，得出了對(duì)應(yīng)于軌跡文件的成鍵圖，如圖5b所示.從圖中可以看到C-C鍵的斷裂較容易形成C2O2團(tuán)簇，隨后在該團(tuán)簇上的C-O鍵斷裂形成CO2. 其次，反應(yīng)過(guò)程中小分子團(tuán)簇C3O1吸附到石墨表面，隨后C-O鍵斷裂一個(gè)單獨(dú)的氧原子從該團(tuán)簇中釋放，最后C-C鍵斷裂也會(huì)形成CO2.該形成過(guò)程與Paci等人［17］和Sriram等人［21］闡述的CO2的形成過(guò)程基本一致. 它們都是由于O原子的碰撞導(dǎo)致形成了環(huán)氧化物，最后C-C鍵斷裂釋放出CO2. 也就是說(shuō)這些團(tuán)簇是形成CO2的主要來(lái)源之一.

C4O1的形成主要源于團(tuán)簇C10O3裂解. 高溫下團(tuán)簇C10O3分解產(chǎn)生兩個(gè)小的團(tuán)簇C6O2和C4O1，隨后兩個(gè)團(tuán)簇結(jié)合成一個(gè)新的團(tuán)簇C7O3釋放出一個(gè)碳原子，最后該碳原子與附近的兩個(gè)氧原子結(jié)合生成了CO2. 該過(guò)程與Paci等人［17］所觀察到CO2的形成過(guò)程基本一致，都是通過(guò)大的團(tuán)簇釋放出CO，然后CO與游離在附近的O原子結(jié)合形成了CO2.

綜上所述，氧原子和石墨相互碰撞形成CO和CO2的過(guò)程中，C原子吸附O原子形成的羰基和氧原子的遷移形成環(huán)氧化物都起著重要作用. 高溫氧原子會(huì)導(dǎo)致形成新的空穴.這些空穴被羰基官能化，從而在隨后與O原子碰撞時(shí)，又成為CO/CO2另一來(lái)源. 另一方面，高溫原子氧與石墨表面具有高度反應(yīng)性. 它傾向于在基面上形成環(huán)氧基團(tuán)，一旦形成，這些基團(tuán)可以遷移，在缺陷處，原子氧可以形成半醌、羰基、內(nèi)酯、羧酸和醚. 這些官能團(tuán)是通過(guò)實(shí)驗(yàn)提出的，這些實(shí)驗(yàn)表明，在高溫下氧原子在石墨表面發(fā)生反應(yīng)時(shí)，C-O、C=O或O-C-O和O-C=O會(huì)填充表面，正是這些基團(tuán)被認(rèn)為是實(shí)驗(yàn)中觀察到的CO和CO2的來(lái)源. 以上類(lèi)型的反應(yīng)也被認(rèn)為在燃燒化學(xué)中發(fā)揮著重要作用.

4.2 反應(yīng)速率的分析

本文通過(guò)分析反應(yīng)過(guò)程中的主要化學(xué)鍵隨時(shí)間的變化情況，進(jìn)一步研究了石墨的氧化動(dòng)力學(xué)性質(zhì)［30］. 圖6分別顯示了三種模型在氧化反應(yīng)過(guò)程中主要化學(xué)鍵的數(shù)量隨時(shí)間的變化情況. 隨著反應(yīng)的進(jìn)行，C-C鍵的數(shù)量逐漸減少，而C-O鍵則相應(yīng)的形成.這些化學(xué)鍵的數(shù)量隨時(shí)間的變化遵循指數(shù)形式，并且可以擬合以下指數(shù)函數(shù).

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

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

其中，C0和C∞分別代表C-C鍵和C-O鍵的數(shù)量，b為擬合系數(shù). 眾所周知，化學(xué)反應(yīng)速率密切依賴于相應(yīng)反應(yīng)物和產(chǎn)物的濃度. 在石墨的氧化反應(yīng)過(guò)程中，C-C鍵為主要裂解的化學(xué)鍵，C-O鍵為主要生成的鍵，而所涉及的化學(xué)鍵的裂變和形成的反應(yīng)速率一般可以用以下的公式表示.

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

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

其中，r是反應(yīng)速率，k是有效反應(yīng)速率常數(shù)，n是反應(yīng)級(jí)數(shù). 將公式（5）（6）轉(zhuǎn)化為對(duì)數(shù)形式（見(jiàn)公式（7）（8）），得到相應(yīng)的線性函數(shù)：ln（-dC（t）/dt）～ln（C（t））h和ln （dC（t）/dt）～ln（C∞- C（t））. 它們的斜率和截距是分別對(duì)應(yīng)n和ln（k）.將公式（3）（4）代入公式（7）（8），發(fā)現(xiàn)對(duì)于反應(yīng)物的數(shù)量（裂變反應(yīng)）和產(chǎn)物（生成反應(yīng)），化學(xué)鍵的裂變和形成是一級(jí)反應(yīng)（n=1），速率常數(shù)k就是擬合系數(shù)b，可以直接通過(guò)化學(xué)鍵的數(shù)量演化的指數(shù)擬合得到.

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

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

圖6a～6c顯示了三種模型在氧化反應(yīng)過(guò)程中C-C鍵的裂變和C-O鍵的生成的指數(shù)擬合. 考慮到在1500 K時(shí)反應(yīng)C-C鍵的裂變和C-O鍵的生成較為不明顯，因此只研究了2000、2500和3000 K時(shí)的數(shù)量變化. 我們可以看出擬合曲線與數(shù)據(jù)的一致性，說(shuō)明擬合合理可靠. 為了進(jìn)一步揭示動(dòng)力學(xué)特性，本文使用經(jīng)典的Arrhenius形式（k=A*exp（-Ea/RT））來(lái)擬合這些數(shù)據(jù)，其中k是有效的速率常數(shù)，A是頻率因子，Ea是活化能，R是氣體常數(shù)，T是溫度（K）. C-C鍵和C-O鍵的反應(yīng)速率常數(shù)的溫度依賴性由Arrhenius公式的對(duì)數(shù)形式ln（k）=ln（A）-Ea/RT表示，如圖7a和7b所示. 三種模型在C-C鍵的裂解上計(jì)算得到的反應(yīng)活化能分別為7.56、2.4和1.6 kcal/mol. 從圖7a中可以看出， 單層無(wú)缺陷模型具有最高的反應(yīng)活化能，而另外兩種模型由于一部分氧原子代替了碳原子和一部分碳原子的缺失，導(dǎo)致活化能降低. 這說(shuō)明無(wú)缺陷模型的C-C鍵不易裂解，因此具有最低的反應(yīng)速率. 三種模型在C-O鍵的生成上計(jì)算得到的反應(yīng)活化能分別為5.15、3.57和4.82 kcal/mol. 從圖7b中可以看出， 單層無(wú)缺陷模型也具有最高的活化能，而另外兩種模型由于碳原子的缺失導(dǎo)致活化能降低，使反應(yīng)速率加快，說(shuō)明C-O鍵的形成是CO和CO2生成速率的主要的控制因素. 在這樣的溫度條件下，石墨在分解過(guò)程中會(huì)留下許多C-C鍵. 這就是在反應(yīng)過(guò)程中產(chǎn)生許多碳原子團(tuán)簇的原因. 此外，值得注意的是，C-C鍵的裂解是氧化反應(yīng)的主要途徑，而C-O鍵的生成是CO和CO2生成的主要控制因素，將顯著影響碳的氧化形成最終的氣態(tài)氧化物. 這與前面在反應(yīng)機(jī)理中觀察到的現(xiàn)象一致，進(jìn)一步驗(yàn)證了C-C鍵的裂解是氧化反應(yīng)主要反應(yīng)途徑，C-O鍵是CO和CO2生成速率的主要控制因素.

4.3 小分子產(chǎn)物的分析

實(shí)驗(yàn)上顯示石墨在氧化還原反應(yīng)過(guò)程的主要產(chǎn)物是CO和CO2.圖8所示為單層無(wú)缺陷石墨模型CO和CO2數(shù)量隨時(shí)間的變化. 從圖8a中可以看出， CO的數(shù)量呈現(xiàn)一定的規(guī)律性. 當(dāng)溫度為3000 K時(shí)，CO數(shù)量最多; 但當(dāng)溫度為1500 K時(shí)，只有少量的CO產(chǎn)生. 從CO產(chǎn)生的數(shù)量上可以看出，在高溫下石墨的氧化還原反應(yīng)是以CO為主導(dǎo)的. 當(dāng)溫度越高，產(chǎn)生的CO數(shù)量越多，反之則越少. 該模型在反映中出現(xiàn)的CO和CO2的數(shù)量上的趨勢(shì)， 與國(guó)義軍等人［31， 32］在碳基燒蝕材料的快慢反應(yīng)理論和燒蝕機(jī)理中計(jì)算的趨勢(shì)相吻合. 即， 在1700 K以上，溫度越高生成物主要是CO; 1500 K以下，溫度越低生成物主要是CO2. 從圖8b中可以看出， 在低溫時(shí)只有在反應(yīng)后期才會(huì)產(chǎn)生1～2個(gè)CO2; 當(dāng)溫度升高，CO2會(huì)逐漸消失. 在高溫時(shí)，反應(yīng)初始階段會(huì)產(chǎn)生1～2個(gè)CO2，在反應(yīng)的后期階段就不再會(huì)產(chǎn)生CO2. 由此可以看出， CO2在低溫下會(huì)較早產(chǎn)生，隨著反應(yīng)的進(jìn)行則不會(huì)再產(chǎn)生; 而在高溫時(shí)CO產(chǎn)生的數(shù)量較多，占主導(dǎo)地位.

根據(jù)上述單層無(wú)缺陷的石墨模型的分析，我們還做了單層有原子缺陷（圖8b）和孔洞缺陷（圖8c）的石墨模型進(jìn)行比較. 從圖中可以看到，單層有缺陷（圖8b）的石墨模型和無(wú)缺陷（圖8a）的石墨模型在分解產(chǎn)物中存在一定的差異. 單層原子缺陷的石墨在氧化反應(yīng)過(guò)程中產(chǎn)生的CO增長(zhǎng)較快，由于石墨表面的一部分碳原子被氧原子所替代，加快了氧化還原反應(yīng)，故CO增長(zhǎng)的速度較快. 所以缺陷對(duì)石墨模型的影響是不能忽略的. 實(shí)際材料中缺陷分為原子缺陷和空位缺陷. 圖8c為單層有孔洞的石墨模型CO和CO2數(shù)量隨時(shí)間的變化情況. 與上述的單層有原子缺陷的石墨模型相比，空位缺陷所產(chǎn)生的CO數(shù)量更多，并且空位缺陷所產(chǎn)生CO數(shù)量的趨勢(shì)比原子缺陷所產(chǎn)生的較為穩(wěn)定. 在CO的產(chǎn)生上，兩者都是在低溫下產(chǎn)生的數(shù)量較為穩(wěn)定，在高溫下都呈現(xiàn)緩慢增長(zhǎng)的狀態(tài). 而CO2的產(chǎn)生兩種模型存在差異，即是在反應(yīng)的初期，有空位缺陷的模型會(huì)產(chǎn)生CO2，有原子缺陷的模型沒(méi)有產(chǎn)生CO2. 因此缺陷的種類(lèi)也會(huì)影響石墨氧化還原反應(yīng)產(chǎn)生的CO和CO2數(shù)量.

5 結(jié) 論

本文采用基于緊束縛密度泛函理論分子動(dòng)力學(xué)方法研究高溫下三種石墨模型（無(wú)缺陷、原子缺陷以及孔洞缺陷）的氧化還原反應(yīng). 研究發(fā)現(xiàn)三種石墨模型的共同點(diǎn)是， 在高溫下CO和CO2是石墨在氧化還原反應(yīng)中的主要產(chǎn)物.CO的主要來(lái)源是由環(huán)氧基團(tuán)上的C-C鍵斷裂后，碳原子與氧原子結(jié)合形成CO;CO2的形成較為復(fù)雜，主要來(lái)源是小分子團(tuán)簇（C2O2、C3O1和C4O1）的裂解. 通過(guò)分析氧化反應(yīng)速率，C-C鍵斷裂是石墨燒蝕反應(yīng)的主要途徑，C-O鍵的斷裂是CO和CO2生成速率的主要控制因素. 此外，三種模型氧化反應(yīng)差異是，在高溫下無(wú)缺陷的石墨模型產(chǎn)生CO和CO2的數(shù)量較多，而有缺陷和有孔洞的模型產(chǎn)生數(shù)量較少，表明缺陷和孔洞在高溫下會(huì)影響CO和CO2產(chǎn)生的數(shù)量. 在氧化反應(yīng)速率方面，三種模型計(jì)算得到的反應(yīng)活化能分別為7.56、2.4和1.6 kcal/mol，表明無(wú)缺陷的石墨模型活化能高，表現(xiàn)出最低的燒蝕反應(yīng)速率，而原子缺陷和孔洞缺陷則會(huì)加快燒蝕反應(yīng)速率. 本文的研究對(duì)燒蝕材料的性能改進(jìn)具有參考意義.

參考文獻(xiàn)：

［1］ 邢亞娟， 孫波， 高坤， 等. 航天飛行器熱防護(hù)系統(tǒng)及防熱材料研究現(xiàn)狀［J］. 力學(xué)學(xué)報(bào)， 2011， 43： 289.

［2］ 胡連成. 前蘇聯(lián)輕質(zhì)防熱材料和碳/酚醛復(fù)合材料［J］. 宇航材料工藝， 2018， 48： 9.

［3］ 靳亮， 王德， 王耀輝， 等. 熱防護(hù)材料燒蝕性能測(cè)試方法研究［J］. 固體火箭技術(shù)， 2018， 41： 245.

［4］ 吳曉宏， 陸小龍， 李濤， 等. 輕質(zhì)燒蝕材料研究綜述［J］. 航天器環(huán)境工程， 2011， 28： 313.

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

［6］ 張威， 曾明， 肖凌飛， 等. 碳-酚醛材料燒蝕熱解對(duì)再入流場(chǎng)特性影響的數(shù)值計(jì)算［J］. 國(guó)防科技大學(xué)學(xué)報(bào)， 2014， 36： 41.

［7］ 高俊杰， 俞繼軍， 韓海濤， 等. 樹(shù)脂基燒蝕材料細(xì)觀傳熱特性預(yù)測(cè)［J］. 航空學(xué)報(bào)， 2017， 38： 170.

［8］ 王子龍， 邢素麗， 尹昌平， 等. 高殘?zhí)磕蜔g樹(shù)脂基體研究進(jìn)展［J］. 工程塑料應(yīng)用， 2018， 46： 131.

［9］ 賀廣東， 鄒寧宇. 對(duì)高硅氧纖維耐燒蝕材料的幾點(diǎn)探討［J］. 玻璃纖維， 2015（6）： 40.

［10］ 霍肖旭， 曾曉梅， 劉紅林. 炭/炭復(fù)合材料及炭纖維增強(qiáng)塑料在固體火箭噴管上的應(yīng)用［J］. 炭素技術(shù)， 2001（3）： 23.

［11］ 李林杰， 方國(guó)東， 易法軍， 等. 高硅氧/酚醛復(fù)合材料熱變形實(shí)驗(yàn)測(cè)試及表面燒蝕形貌分析［J］. 固體火箭技術(shù)， 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］ 楊德軍. 防熱復(fù)合材料燒蝕行為的數(shù)值模擬［D］. 蘭州： 蘭州理工大學(xué)， 2013.

［24］ 劉驍. 燒蝕熱防護(hù)系統(tǒng)三維熱響應(yīng)有限元計(jì)算方法研究［D］. 綿陽(yáng)：中國(guó)空氣動(dòng)力研究與發(fā)展中心， 2016.

［25］ 國(guó)義軍， 石衛(wèi)波， 曾慕， 等．高超聲速飛行器燒蝕防熱理論與應(yīng)用［M］． 北京： 科學(xué)出版社， 2019.

［26］ 李偉， 方國(guó)東， 李瑋潔， 等. 碳纖維增強(qiáng)復(fù)合材料微觀燒蝕行為數(shù)值模擬［J］. 力學(xué)學(xué)報(bào)， 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鋼的脫碳行為演化及防護(hù)研究［J］. 四川大學(xué)學(xué)報(bào)： 自然科學(xué)版， 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］ 國(guó)義軍， 周述光， 曾磊， 等. 碳基材料氧化燒蝕的快/慢反應(yīng)動(dòng)力學(xué)問(wèn)題［J］. 實(shí)驗(yàn)流體力學(xué)， 2019， 33： 17.

［32］ 國(guó)義軍， 桂業(yè)偉， 童福林， 等. 碳/碳化鋯復(fù)合材料燒蝕機(jī)理和計(jì)算方法研究［J］. 空氣動(dòng)力學(xué)報(bào)， 2012， 31： 22.

引用本文格式：

中 文： 李磊， 劉驍， 周述光， 等. 碳基燒蝕材料高溫氧化反應(yīng)機(jī)理的分子動(dòng)力學(xué)研究［J］. 四川大學(xué)學(xué)報(bào)： 自然科學(xué)版， 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.