氧在熔融SiO2中擴散的分子動力學(xué)模擬

時曉婷，高希光，鐘 毅

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

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

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

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

1 分子動力學(xué)模型

1.1 模型構(gòu)建

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

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

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

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

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

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

式中：ε為勢阱深度；σ為原子作用直徑；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ù)需通過量子力學(xué)從頭算進(jìn)行擬合，根據(jù)Bedra 等［17?18］對高溫下氧在石英表面復(fù)合的研究，考慮到氧在SiO2表面可能發(fā)生解離或吸附，做出如下合理假設(shè)：

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

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

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

1.3 模擬方法

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

2 計算與分析

2.1 計算方法

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

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

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

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

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

2.2 模型驗證

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

圖2 高溫?zé)o定形SiO2結(jié)構(gòu)Fig.2 Structure of high temperature amorphous SiO2

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

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

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

表3 SiO2結(jié)構(gòu)特性計算結(jié)果Table 3 Calculation results of SiO2 structural feature

2.3 擴散系數(shù)分析

在此結(jié)構(gòu)基礎(chǔ)上，分別對950、1 100、1 200、1 300 及1 400 ℃下氧氣在SiO2中的擴散進(jìn)行了分子動力學(xué)模擬。以950 ℃情況為例，圖4 描述了不同時刻下的體系擴散情況，圖5 為系統(tǒng)在不同溫度下運行6 ns 時的擴散狀態(tài)對比。根據(jù)圖示，隨著模擬時間的加長，氧氣在SiO2中的分布逐漸均勻，在10 ns 時氧氣近似均勻分布在SiO2中。在其他溫度下，也表現(xiàn)出了相同的擴散規(guī)律。隨著溫度的升高，達(dá)到平衡狀態(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)抽樣時的氧氣在體系中的分布細(xì)節(jié)，可以看出在擴散達(dá)到穩(wěn)定后，SiO2內(nèi)部氧仍以分子形式存在，沒有出現(xiàn)晶格替換等其他擴散形式。

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

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

表4 分子動力學(xué)模擬結(jié)果Table 4 Results of MD simulation

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

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

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

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

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

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

圖8 Arrhenius 關(guān)系曲線Fig.8 Arrhenius curve

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

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

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

3 結(jié)論

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

（1）模擬高溫?zé)o定形SiO2結(jié)構(gòu)，將結(jié)構(gòu)參數(shù)與實驗數(shù)據(jù)進(jìn)行對比，驗證模型準(zhǔn)確性；

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

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