時曉婷，高希光，鐘 毅

（1.南京航空航天大學能源與動力學院，航空發(fā)動機熱環(huán)境與熱結構工業(yè)和信息化部重點實驗室，南京 210016；2.南京航空航天大學能源與動力學院，江蘇省航空動力系統(tǒng)重點實驗室，南京 210016；3.中國海誠工程科技股份有限公司，上海 200031）

碳化硅（SiC）陶瓷材料在1 000 ℃以上具有耐磨損、耐腐蝕、機械性能優(yōu)等特點，其制備的陶瓷基復合材料（Ceramic matrix composites，CMCs）輕質(zhì)、耐高溫、高強韌，現(xiàn)已廣泛應用于航空航天領域，成為先進發(fā)動機熱端部件最具有潛力的替代材料。但SiC 在900～1 700 ℃的高溫足氧環(huán)境下將與氧氣發(fā)生被動反應［1］，表面形成SiO2氧化物，影響SiC 的高溫力學性能。研究表明，SiO2在SiC 上的生長由氧氣通過氧化物的擴散控制［2］。因此，研究氧在高溫熔融SiO2中的擴散規(guī)律，對于判斷氧氣到達未氧化界面的速率、反映氧化膜的保護作用都具有重要的理論意義。同時，確定擴散系數(shù)的經(jīng)驗公式，可用于求解結構件中的氧化程度，為SiC?CMCs 在航空發(fā)動機上的工程應用提供計算所需的初始參數(shù)。

擴散系數(shù)作為描述擴散過程的一個重要參數(shù)，其測定主要有實驗、經(jīng)驗公式和理論計算3 種方法。Norton［3］采用由球壁分隔的滲透室、電容、質(zhì)譜儀設備等測定了950～1 080 ℃范圍內(nèi)氧氣在玻璃狀石英內(nèi)的擴散系數(shù)。為保證實現(xiàn)真正的擴散而非脫氣，需實時觀測氣體壓強并保證穩(wěn)態(tài)流。試驗結果較好地擬合了Arrhenius 圖。Sucov［4］根據(jù)氣體交換技術，利用18O 同位素示蹤法及質(zhì)譜儀，測定了925～1 225 ℃溫度范圍內(nèi)氧離子在玻璃狀石英中的擴散率，并用Arrhenius 方程描述了實驗測量結果。Kalen 等［5］利用氣體交換技術研究了氧在玻璃石英中的擴散機制，實驗結果與核反應分析及二次離子質(zhì)譜技術等進行對比，認為擴散系數(shù)結果的分散性可定性解釋為存在兩種擴散機制。上述研究均以實驗方式對氧氣在SiO2中的擴散系數(shù)開展了測定，給出了相關的實驗現(xiàn)象，涉及各類先進儀器設備，實驗成本較高，且受制于實驗條件，測定的溫度范圍有限；另外，氧氣在高溫SiO2中主要以分子形式擴散，僅當溫度超過1 400 ℃時才以離子形式擴散。因此，在前人實驗研究的基礎上，有必要發(fā)展相關的仿真計算，用以擴展實驗區(qū)域以外的預測，同時相互支撐，進一步完善反應擴散機制的分析。

近幾十年來，隨著計算機技術的發(fā)展，分子動力學（Molecular dynamics，MD）模擬提供了一種新的計算分析方法，與理論、實驗研究相互促進、補充。分子動力學模擬跟蹤體系隨時間的動態(tài)演化，描述微小時間尺度內(nèi)粒子的動態(tài)擴散行為，可以從微觀上獲得粒子擴散的動力學信息，進行大量數(shù)據(jù)的反復統(tǒng)計與分析處理。通過模擬，不僅可以得到給定模型的擴散系數(shù)，以檢驗近似理論解析解的有效性，還可以突破實驗條件限制，計算高溫高壓等極端條件下的擴散系數(shù)。對于其他物質(zhì)在不同介質(zhì)中的擴散過程研究，國內(nèi)外已開展了相關工作并進行了計算模擬［6?9］。因此，本文基于上述學者的研究方法，利用Lammps 分子動力學開源模擬軟件，對于氧在950～1 400 ℃下熔融SiO2中的擴散進行了一系列的模擬，計算高溫下的擴散系數(shù)，具體工作包括：（1）高溫無定形SiO2模型的建立；（2）氧與SiO2的相互作用；（3）氧擴散系數(shù)的計算；（4）溫度對擴散系數(shù)的影響作用。

1 分子動力學模型

1.1 模型構建

以β?方石英為起始晶胞，空間群為FD?3M，晶胞參數(shù)為a=b=c=7.160 ?，α=β=γ=90°，進行周期性擴展生成6×6×5 的超晶胞，在表面生成60 個氧分子，共得到3 696 個原子作為初始模型，如圖1 所示。模擬盒子采用周期性邊界條件，選用Velocity Verlet 算法對牛頓運動方程進行求解；靜電力計算使用PPPM 算法，精度為10-4；截斷半徑為9 ?；時間步長為1 fs，調(diào)溫和調(diào)壓使用Nose?Hoover 方法。

1.2 勢函數(shù)及參數(shù)

分子動力學模擬結果的準確性主要取決于是否選取了能夠準確描述分子微觀結構的勢能模型。在本體系中，系統(tǒng)總勢能包括3 部分：SiO2勢能USiO2，O2間勢能UO2，SiO2與O2間相互作用勢能USiO2?O2。

對于SiO2，Si—O 鍵既有離子性成分，又有共價性成分，根據(jù)丁元法等［10］對石英玻璃高溫勢函數(shù)的研究及Govers 等［11］、孫義程［12］在相關模擬中的應用，在本模擬中，選用典型對勢模型Morse 勢疊加庫侖勢描述SiO2間的勢能作用，其函數(shù)形式為

式中：右邊第1 項表示庫侖作用；第2 項為原子對之間的相互作用；rij表示兩原子之間的距離；q表示原子的電荷，由于Si 和O 只有部分離子性，根據(jù)C’ag?in 等［13］的擬合研究，取qSi=+1.3e，qo=-0.65e；D0、γ、ρ為相應的勢函數(shù)參數(shù)。各參數(shù)值列于表1，其中ρ=rij/R0。

對于O2，氧分子由2 個原子通過鍵連接構成。氧原子間的勢能采用Lennard?Jones（12?6）勢能函數(shù)，其表達式為

式中：ε為勢阱深度；σ為原子作用直徑；rij為兩原子之間的距離；rc為原子作用截斷半徑。在截斷半徑處兩原子間勢能為0。

氧分子在研究中一般被視為一種具有彈性的分子，氧分子中原子間鍵的作用常采用一種簡諧勢函數(shù)表示

式中：K為彈性常數(shù)；xc為平衡鍵長。各勢函數(shù)參數(shù)列于表2。

表2 O2勢函數(shù)參數(shù)[15?16]Table 2 Potential function parameter of O2[15?16]

對于SiO2?O2，其精確的相互作用勢及參數(shù)需通過量子力學從頭算進行擬合，根據(jù)Bedra 等［17?18］對高溫下氧在石英表面復合的研究，考慮到氧在SiO2表面可能發(fā)生解離或吸附，做出如下合理假設：

（1）氧氣中的O 與SiO2中Si 的相互作用采用SiO2中Si—O 間的作用勢及參數(shù)；

（2）氧氣中的O 與SiO2中O 的相互作用采用氣相中O—O 間的作用勢及參數(shù)。

以上設定提供了氧進入SiO2晶格、發(fā)生晶格擴散的可能性，能更全面對擴散規(guī)律進行模擬。

1.3 模擬方法

通過初速度的形式將系統(tǒng)的初始溫度（T0=300 K）賦予每個粒子，初速度遵循統(tǒng)計原理并服從Maxwell?Boltzman 分布。首先將氧分子固定，在NPT 系綜下將初始結構升溫至5 000 K，運行20 ps，以去除SiO2初始構型的影響，隨后以100 K/20 ps的速率使系統(tǒng)溫度下降到模擬溫度，得到SiO2高溫無定形結構。以該模型作為起點，釋放氧氣分子，保持溫度恒定對系統(tǒng)進行充分弛豫。在各次的降溫、調(diào)壓過程中，根據(jù)輸出的能量和溫度、壓力值進行統(tǒng)計分析，表明系統(tǒng)均已達到平衡后進行采樣，采樣間隔為100 步（即0.1 ps）。

2 計算與分析

2.1 計算方法

根據(jù)統(tǒng)計力學中的漲落?耗散理論，分子動力學模擬中可用兩個等價的公式確定擴散系數(shù)，一種是利用均方位移（Mean square displacement，MSD）計算的Einstein 關系式

另一種是利用速度自相關函數(shù)（Velocity auto?correlation function，VACF）計算的Green?Kubo公式

式中：D為擴散系數(shù)；ri(t)、ri(0)和vi(t)、vi(0)分別對應粒子t時刻和0 時刻的位移矢量和速度矢量； 表示系綜平均。

Green?Kubo 公式需要對速度自相關函數(shù)進行大量的積分統(tǒng)計才能獲得擴散系數(shù)；而Einstein 公式直接與模擬結果的均方位移具有對應關系，更加簡便快捷。為節(jié)省數(shù)據(jù)處理的時間，在保證計算結果可靠性的前提下，本模擬中通過MSD 計算高溫下氧氣在SiO2氧化膜中的擴散系數(shù)。

為了避免氧氣系綜的壓力對模擬盒子大小的非物理影響，整個系統(tǒng)的動力學計算在NVT 系綜下進行。經(jīng)過一段初始的充分弛豫，氧氣分子的當前坐標被存儲起來。進一步的模擬伴隨著氧氣分子相對于這些坐標的均方位移的周期性計算。對于抽樣時間內(nèi)采集到的MSD 數(shù)據(jù)，繪制MSD 隨時間變化的曲線圖，當曲線呈線性增長時，可根據(jù)Einstein 公式計算氧氣的擴散系數(shù)。

2.2 模型驗證

SiO2結構的靜態(tài)性質(zhì)是后續(xù)動態(tài)研究的基礎，為驗證模型的有效性，首先對SiO2結構進行分析。圖2 為分子動力學模擬高溫無定形SiO2體系的結構，模型中原子呈現(xiàn)長程無序結構，與圖1 相比沒有了晶體結構中的原子三維周期性分布。

圖2 高溫無定形SiO2結構Fig.2 Structure of high temperature amorphous SiO2

為進一步考察高溫無定形SiO2的結構特性，計算平衡時的Si—Si、O—O、Si—O 對關聯(lián)函數(shù)如圖3 所示。對關聯(lián)函數(shù)表示兩個粒子屬性之間的空間相關性，使用快速傅里葉變換算法（FFT）計算卷積，然后在互反空間和實空間中計算徑向平均。通過對關聯(lián)函數(shù)分析，可以了解粒子間作用的相對強弱，第一峰位置對應為原子對間平衡距離。

圖3 分別為950 ℃與1 400 ℃下的對關聯(lián)函數(shù)圖，曲線除微小波動近似一致，由圖線可以清楚看出體系呈現(xiàn)近程有序、遠程無序結構，第一峰峰位計算結果與文獻［19］中的試驗結果對比情況列于表3，結果表明結構符合SiO2高溫無定形特性。

圖3 SiO2對關聯(lián)函數(shù)Fig.3 Correlation functions of SiO2 pair

表3 SiO2結構特性計算結果Table 3 Calculation results of SiO2 structural feature

2.3 擴散系數(shù)分析

在此結構基礎上，分別對950、1 100、1 200、1 300 及1 400 ℃下氧氣在SiO2中的擴散進行了分子動力學模擬。以950 ℃情況為例，圖4 描述了不同時刻下的體系擴散情況，圖5 為系統(tǒng)在不同溫度下運行6 ns 時的擴散狀態(tài)對比。根據(jù)圖示，隨著模擬時間的加長，氧氣在SiO2中的分布逐漸均勻，在10 ns 時氧氣近似均勻分布在SiO2中。在其他溫度下，也表現(xiàn)出了相同的擴散規(guī)律。隨著溫度的升高，達到平衡狀態(tài)的時間越短，氣體擴散越劇烈。

圖4 950 ℃下系統(tǒng)的擴散過程Fig.4 Diffusion of the system at 950 ℃

圖5 不同溫度下運行6 ns 的擴散狀態(tài)Fig.5 Diffusion state at different temperatures of 6 ns

圖6 為終態(tài)抽樣時的氧氣在體系中的分布細節(jié)，可以看出在擴散達到穩(wěn)定后，SiO2內(nèi)部氧仍以分子形式存在，沒有出現(xiàn)晶格替換等其他擴散形式。

圖6 氧氣分布細節(jié)圖Fig.6 Distribution detail of oxygen in the system

對10 ns 內(nèi)的模擬過程進行采集，計算該過程中氧氣的均方位移并輸出至文本，對計算結果進行統(tǒng)計平均作為計算結果。圖7 繪制了氧氣在不同溫度下均方位移隨時間的變化關系。從圖像上可以看出，MSD 變化趨勢比較明顯，在前1 ns 內(nèi)呈拋物線增長，后呈近似線性增長。根據(jù)Einstein 公式計算擴散系數(shù)最終在一常數(shù)附近小幅漲落，對平衡后的結果取平均，計算結果列于表4。

表4 分子動力學模擬結果Table 4 Results of MD simulation

圖7 氧氣均方位移曲線Fig.7 MSD curves of oxygen

計算結果顯示，在不考慮更高溫度下產(chǎn)生結晶的情況下，隨著溫度不斷升高，氧在SiO2中的擴散系數(shù)逐漸增大，預示著氧氣穿過保護層的速率越快，能夠達到未氧化區(qū)域的氧氣也越多。因此，SiO2的保護作用可能會被削弱。

為了驗證計算結果的可靠性，對氧氣在SiO2中的擴散系數(shù)與其他文獻［4?5，20］中的數(shù)據(jù)進行了比較。本文計算得到擴散系數(shù)為10-7cm2/s，文獻［20］總結了前人采用同位素交換等各類方式測定的試驗值，由于存在考慮的溫度范圍及擴散機制差異等因素，測量值分散性較大，在10-6～10-8cm2/s范圍內(nèi)。而分子動力學方法經(jīng)多次模擬得到了較穩(wěn)定的計算值，且在試驗測量值范圍內(nèi)。另外，本文模擬的擴散氧氣以分子形式進行，計算結果可用于與實驗的相互補充，通過擴散能的比較，分析擴散機制。

根據(jù)Arrhenius 定律，擴散系數(shù)與溫度存在如下關系式

式中：D為擴散系數(shù)，D0為頻率因子，Q為擴散激活能，R為氣體常數(shù)，T為絕對溫度。

對式（6）取對數(shù)可轉(zhuǎn)化為線性方程，對數(shù)據(jù)點進行線性分析繪制圖線如圖8 所示，擬合結果為

圖8 Arrhenius 關系曲線Fig.8 Arrhenius curve

式中：D0=2.6×10-4cm2/s，Q=74.1 kJ/mol。

為驗證模擬結果的可靠性，與各文獻［5，21］中的計算或試驗結果進行比較。擴散激活能Q表示克服能壘實現(xiàn)原子從一個平衡位置躍遷到另一個平衡位置的基本躍遷的額外能量，直接考察了擴散能力，因此對比Q情況列于表5。Newsome 等［21］采用反應力場模擬了SiC 氧化及氧氣擴散過程，其計算結果與本文計算結果誤差較小。Muehlenbachs等［22］、Norton［3］利用示蹤原子對氧氣的擴散進行了試驗測定，其計算結果高于本文計算結果。這樣的誤差，一方面是由于示蹤原子的摩爾質(zhì)量高于實際氧氣中氧原子的摩爾質(zhì)量，另一方面，氧氣在SiO2中的擴散存在表面吸附過程，即存在解離能，而分子動力學模擬中沒有出現(xiàn)該情況，因此激活能偏小。

表5 Arrhenius 參數(shù)擬合結果對比Table 5 Comparison of Arrhenius parameter fitting re?sults

3 結論

本文利用Lammps 軟件對高溫下熔融SiO2中的氧擴散進行了分子動力學模擬，計算了不同溫度下的氧氣擴散系數(shù)。主要工作及結論如下：

（1）模擬高溫無定形SiO2結構，將結構參數(shù)與實驗數(shù)據(jù)進行對比，驗證模型準確性；

（2）模擬氧在SiO2中的擴散過程，分析擴散規(guī)律，計算了不同溫度下的擴散系數(shù)，其溫度依賴性可用所擬合的Arrhenius 擴散速率方程進行描述；

（3）上述計算結果可與試驗結果相互支撐、補充，為擴散控制的SiC 基及其復合材料氧化行為研究提供理論參考。