碳基燒蝕材料高溫氧化反應(yīng)機理的分子動力學(xué)研究

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

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

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

TB35A2023.034003

收稿日期： 2022-10-24

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

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

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

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

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

2 模擬方法

本文使用CP2K代碼中實現(xiàn)的基于量子力學(xué)的MD方法模擬了石墨的氧化還原反應(yīng). 自洽電荷緊束縛密度泛函理論（SCC-DFTB）方法用于計算原子間力，它是基于密度泛函理論（DFT）的總能量相對于電荷密度的二階展開，能夠以自洽的方式描述能量隨電荷密度變化的情況. 為了準確描述能量隨電荷密度的變化情況，SCC-DFTB所描述的電荷密度是體系能量對電荷密度變化泰勒級數(shù)展開的二階導(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是在參考密度下計算的有效的科恩-沙姆哈密頓量，其中n（r）為總電荷密度，n0（r）為參考密度的二階導(dǎo)數(shù). φi是科恩-沙姆軌道，Exc和Vxc分別為交換關(guān)系能和勢，Eii為核-核排斥能. 為了推導(dǎo)出SCC-DFTB方法的總能量，對方程（1）中的能量作進一步近似，SCC-DFTB總能量最終可以表達為

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

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

3 模型構(gòu)建

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

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

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

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

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

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

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

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

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

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

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

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

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

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

本文通過分析反應(yīng)過程中的主要化學(xué)鍵隨時間的變化情況，進一步研究了石墨的氧化動力學(xué)性質(zhì)［30］. 圖6分別顯示了三種模型在氧化反應(yīng)過程中主要化學(xué)鍵的數(shù)量隨時間的變化情況. 隨著反應(yīng)的進行，C-C鍵的數(shù)量逐漸減少，而C-O鍵則相應(yīng)的形成.這些化學(xué)鍵的數(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)過程中，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)級數(shù). 將公式（5）（6）轉(zhuǎn)化為對數(shù)形式（見公式（7）（8）），得到相應(yīng)的線性函數(shù)：ln（-dC（t）/dt）～ln（C（t））h和ln （dC（t）/dt）～ln（C∞- C（t））. 它們的斜率和截距是分別對應(yīng)n和ln（k）.將公式（3）（4）代入公式（7）（8），發(fā)現(xiàn)對于反應(yīng)物的數(shù)量（裂變反應(yīng)）和產(chǎn)物（生成反應(yīng)），化學(xué)鍵的裂變和形成是一級反應(yīng)（n=1），速率常數(shù)k就是擬合系數(shù)b，可以直接通過化學(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)過程中C-C鍵的裂變和C-O鍵的生成的指數(shù)擬合. 考慮到在1500 K時反應(yīng)C-C鍵的裂變和C-O鍵的生成較為不明顯，因此只研究了2000、2500和3000 K時的數(shù)量變化. 我們可以看出擬合曲線與數(shù)據(jù)的一致性，說明擬合合理可靠. 為了進一步揭示動力學(xué)特性，本文使用經(jīng)典的Arrhenius形式（k=A*exp（-Ea/RT））來擬合這些數(shù)據(jù)，其中k是有效的速率常數(shù)，A是頻率因子，Ea是活化能，R是氣體常數(shù)，T是溫度（K）. C-C鍵和C-O鍵的反應(yīng)速率常數(shù)的溫度依賴性由Arrhenius公式的對數(shù)形式ln（k）=ln（A）-Ea/RT表示，如圖7a和7b所示. 三種模型在C-C鍵的裂解上計算得到的反應(yīng)活化能分別為7.56、2.4和1.6 kcal/mol. 從圖7a中可以看出， 單層無缺陷模型具有最高的反應(yīng)活化能，而另外兩種模型由于一部分氧原子代替了碳原子和一部分碳原子的缺失，導(dǎo)致活化能降低. 這說明無缺陷模型的C-C鍵不易裂解，因此具有最低的反應(yīng)速率. 三種模型在C-O鍵的生成上計算得到的反應(yīng)活化能分別為5.15、3.57和4.82 kcal/mol. 從圖7b中可以看出， 單層無缺陷模型也具有最高的活化能，而另外兩種模型由于碳原子的缺失導(dǎo)致活化能降低，使反應(yīng)速率加快，說明C-O鍵的形成是CO和CO2生成速率的主要的控制因素. 在這樣的溫度條件下，石墨在分解過程中會留下許多C-C鍵. 這就是在反應(yīng)過程中產(chǎn)生許多碳原子團簇的原因. 此外，值得注意的是，C-C鍵的裂解是氧化反應(yīng)的主要途徑，而C-O鍵的生成是CO和CO2生成的主要控制因素，將顯著影響碳的氧化形成最終的氣態(tài)氧化物. 這與前面在反應(yīng)機理中觀察到的現(xiàn)象一致，進一步驗證了C-C鍵的裂解是氧化反應(yīng)主要反應(yīng)途徑，C-O鍵是CO和CO2生成速率的主要控制因素.

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

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

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

5 結(jié) 論

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

參考文獻：

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

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

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

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

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

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

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

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

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

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

［11］ 李林杰， 方國東， 易法軍， 等. 高硅氧/酚醛復(fù)合材料熱變形實驗測試及表面燒蝕形貌分析［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］ 劉驍. 燒蝕熱防護系統(tǒng)三維熱響應(yīng)有限元計算方法研究［D］. 綿陽：中國空氣動力研究與發(fā)展中心， 2016.

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

［26］ 李偉， 方國東， 李瑋潔， 等. 碳纖維增強復(fù)合材料微觀燒蝕行為數(shù)值模擬［J］. 力學(xué)學(xué)報， 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］. 四川大學(xué)學(xué)報： 自然科學(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］ 國義軍， 周述光， 曾磊， 等. 碳基材料氧化燒蝕的快/慢反應(yīng)動力學(xué)問題［J］. 實驗流體力學(xué)， 2019， 33： 17.

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

引用本文格式：

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