王 權(quán) 莊宿國 黃 丹 朱正興 劉秀波

（1. 中南林業(yè)科技大學(xué)材料表界面科學(xué)與技術(shù)湖南省重點實驗室 湖南長沙 410004；2. 西安航天動力研究所 陜西西安 710100）

隨著機械端面密封在航天航空等領(lǐng)域的廣泛應(yīng)用， 對其在特定工況下的性能提出了更高要求。 其中， 航天領(lǐng)域液氧煤油發(fā)動機的端面密封通常在高pv值（p： 壓力，v： 速度）、 寬轉(zhuǎn)速區(qū)間內(nèi)運行， 并且服役時間長， 容易導(dǎo)致端面密封失效、 介質(zhì)泄漏等問題， 存在一定的安全隱患[1-6］。

為了提高端面密封摩擦副材料的摩擦學(xué)性能及工作穩(wěn)定性， 國內(nèi)外研究者已進行了大量研究和試驗。目前對于端面密封摩擦副性能的優(yōu)化主要集中在其機械結(jié)構(gòu)方面， 改善摩擦副接觸條件或傳動過程， 以匹配各種工況需求[7-9］。 如ZHOU 和ZOU[10］運用有限元分析的手段， 探究了壓力、 轉(zhuǎn)速對端面密封性能的影響。 但是對于該部件微觀尺度下的運動過程及對應(yīng)機制研究較少。

S-07 不銹鋼作為一種航天發(fā)動機用端面密封材料，是一種馬氏體-奧氏體雙相鎳鉻不銹鋼， 具有高強高韌和耐蝕性優(yōu)良等特點， 其化學(xué)成分如表1 所示[11］。

表1 S-07 不銹鋼化學(xué)成分[11］Table 1 Chemical composition of S-07 stainless steel[11］

機械端面密封在服役時， 接觸面緊密接觸且接觸應(yīng)力較大， 此時， 可將微觀尺度下的接觸行為視為諸多粒子相互作用。 分子動力學(xué)（Molecular dynamics，MD） 模擬是一種在納米尺度上研究材料變化的工具， 它基于牛頓經(jīng)典力學(xué)的計算方法， 從統(tǒng)計力學(xué)的基本原理出發(fā)， 能夠在微觀的分子甚至是原子層面，在計算機中仿真得到能量、 溫度、 應(yīng)力等物理量， 進而推導(dǎo)出工程問題所需的摩擦學(xué)參數(shù)[12-14］。 MD 模擬因其能從原子運動角度對材料體系進行系統(tǒng)模擬和計算， 現(xiàn)已廣泛運用于研究醫(yī)藥、 化學(xué)、 材料表面工程等科學(xué)研究領(lǐng)域[15-19］。 張宏亮等[20］利用MD 模擬通過改變滑動速度、 距離及外加載荷， 研究了納米單晶銅的磨料磨損行為， 發(fā)現(xiàn)納米單晶銅內(nèi)部缺陷及表層單晶銅原子的磨料磨損行為有較大差異。 YANG等[21］通過MD 模擬軟件， 模擬了AlCoCrFe 高熵合金涂層的力學(xué)和摩擦學(xué)性能， 模擬測試計算的彈性模量、 納米硬度、 摩擦因數(shù)和磨損體積均與試驗所得數(shù)據(jù)保持一致。 LI 等[22］采用分子動力學(xué)方法研究了銅基高熵合金涂層的摩擦磨損行為， 運用位錯萃取分析（Dislocation extraction analysis， DXA） 和共近鄰分析（Common neighbor analysis， CNA） 分析晶格結(jié)構(gòu)和位錯， 發(fā)現(xiàn)高熵合金涂層可以有效地釋放應(yīng)力， 降低銅基材料的損傷。

目前已經(jīng)有很多學(xué)者采用MD 模擬對各種材料的摩擦學(xué)性能進行了研究， 但大多數(shù)模擬過程中， 都是采用金剛石剛體模型等作為對偶件， 研究摩擦過程中原子運動、 位錯及應(yīng)力等信息。 但是， 結(jié)合工件實際工作場景， 選擇常見的對偶件結(jié)構(gòu)模型進行模擬， 可以使模擬結(jié)果對于實際應(yīng)用更具有直觀參考價值。

本文作者將分子動力學(xué)模擬技術(shù)應(yīng)用于S-07 不銹鋼端面密封材料體系， 基于發(fā)動機工況特征， 選擇同種材料作為摩擦副， 通過改變摩擦過程中壓入深度、 滑動速度等參數(shù)， 模擬其磨損量、 摩擦因數(shù)等摩擦學(xué)性能， 探索其摩擦磨損機制， 并基于此提出合理有效控制及改進方案， 以期提高端面密封材料的可靠性。

1 模型和勢函數(shù)

在S-07 不銹鋼分子動力學(xué)模擬建模過程中， 從S-07 不銹鋼的組成元素Fe、 Cr 和Ni 入手， 建立密封材料S-07 不銹鋼的結(jié)構(gòu)模型。 具體操作過程如下：

（1） 建立S-07 不銹鋼的原胞模型， 按照質(zhì)量分數(shù)為78.24% （Fe）、 16.45% （Cr）、 5.31% （Ni）將Fe、 Cr、 Ni 這3 種原子隨機填充到點陣中；

（2） 對建立好的模型進行幾何結(jié)構(gòu)優(yōu)化以及能量最小化處理；

（3） 選擇優(yōu)化結(jié)束后能量最低的S-07 不銹鋼原胞模型擴胞建立S-07 不銹鋼模型。

為符合所建立結(jié)構(gòu)模型的部分設(shè)定， 消除邊界效應(yīng)， 在其X、Y軸方向上設(shè)定周期性邊界條件， 在Z軸方向上設(shè)定非周期性邊界條件， 模型尺寸大小為17.873 8 nm×8.419 6 nm×3.889 8 nm， 共有43 200 個原子， 其中包含F(xiàn)e 原子33 480 個， Ni 原子7 560 個，Cr 原子2 160 個， 對偶件模型尺寸大小為5.957 9 nm×2.806 5 nm×1.944 9 nm， 將其設(shè)置為剛體。

經(jīng)過上述過程， 建立的目標(biāo)S-07 不銹鋼的結(jié)構(gòu)模型， 如圖1 所示。

圖1 S-07 不銹鋼分子動力學(xué)模擬結(jié)構(gòu)模型Fig.1 S-07 stainless steel molecular dynamics simulation structural model

隨著MD 模擬技術(shù)的不斷發(fā)展， 已經(jīng)有較為成熟的勢函數(shù)來描述金屬原子間的相互作用， 如嵌入原子勢（Embedded-atom method， EAM） 等。 文中模擬計算時， 采用EAM 勢函數(shù)描述模型中Fe-Cr-Ni 原子之間的相互作用[23］。 其勢函數(shù)如下：

式中：Ei為系統(tǒng)的總勢能；F為原子的電子密度ρ的函數(shù)；ρi為除i外其余原子在i處產(chǎn)生的電子云密度的和；?為對勢項；rij為i原子和j原子之間的距離。

2 數(shù)值模擬

MD 模擬具體步驟如下： （1） 建立所需結(jié)構(gòu)模型； （2） 優(yōu)化結(jié)構(gòu)模型參數(shù)； （3） 給定初始條件；（4） 開始模擬， 計算模擬運動過程中的宏觀物理量。文中研究采用LAMMPS 模擬計算軟件、 OVITO 可視化軟件， 實現(xiàn)MD 模擬， 并進行可視化處理， 便于直觀研究模擬的進程。

為了有效利用計算資源， 選取時間步長為1 fs?？紤]到端面密封摩擦副材料的摩擦磨損過程， 選取相似于S-07 不銹鋼模型的另一模型作為對偶件， 以壓入深度和滑動速度為變量， 進行摩擦磨損試驗?zāi)M，分析S-07 不銹鋼的摩擦磨損性能。 在模擬加載之前，要對構(gòu)建的結(jié)構(gòu)模型進行優(yōu)化， 采用共軛梯度算法（Conjugate gradient algorithm） 最小化整個結(jié)構(gòu)模型的能量， 消除模型中不合理的結(jié)構(gòu)。 為了使模型達到平衡狀態(tài)， 采用正則系綜（Canonical ensemble， NVT）方法， 將模型初始溫度設(shè)為300 K 以用于后續(xù)計算。

對S-07 不銹鋼的摩擦磨損性能模擬分為兩個階段。 在第一階段， 對偶件以0.01 nm／ps 的速度進入模型， 達到預(yù)定壓入深度； 在第二階段， 以預(yù)定的滑動速度沿X正方向滑動。 在這兩個階段中， 模型整體溫度保持在300 K， 在NVT 系綜下進行模擬。 仿真結(jié)果使用OVITO 軟件對模型進行可視化， 模型及仿真環(huán)境參數(shù)如表2 所示。

表2 模型和仿真條件參數(shù)Table 2 Model and simulation environment parameters

3 實驗結(jié)果與分析

3.1 摩擦因數(shù)分析

對偶件受到的力是材料磨損性能的具體體現(xiàn)， 最重要的性能值就是摩擦因數(shù)， 而摩擦因數(shù)與對偶件的切向力和法向力密切相關(guān)， 因此分析摩擦因數(shù)是研究材料減摩性能的一個重要方向。 文中利用LAMMPS軟件中compute 功能， 將對偶件受到的摩擦力與法向力的比計算為摩擦因數(shù)， 如式（3） 所示， 其中μ為摩擦因數(shù)，F(xiàn)x為摩擦力，F(xiàn)z為切向力[24］。

繪制對偶件在不同參數(shù)條件下的摩擦因數(shù)曲線如圖2 所示。

圖2 不同壓入深度下對偶件的平均摩擦因數(shù)隨相對速度的變化Fig.2 Variation of average friction coefficient of the friction pair with relative velocity under different indentation depth

從圖2 中可知， 在0.3、 0.6、 0.9 nm 壓入深度下， 隨著相對滑動速度增加， 摩擦因數(shù)均呈現(xiàn)上升趨勢。 在S-07 不銹鋼與對偶件的相對滑動過程中， 對偶件先壓入表面， 然后與不銹鋼發(fā)生相對滑動； 期間， 對偶件與不銹鋼表面存在剪切作用， 不銹鋼表面發(fā)生彈性和塑性變形。 當(dāng)作用在不銹鋼表面的機械能足以破壞S-07 不銹鋼原子間的金屬鍵后， 不銹鋼原子就會隨對偶件一起滑動， 進而被磨損去除形成磨屑。 摩擦因數(shù)增加是因為在單位時間內(nèi)， 所需要破壞的金屬鍵更多， 從而表現(xiàn)為摩擦因數(shù)上升。

3.2 磨損形貌分析

不同參數(shù)下的模擬磨損表面形貌如圖3 所示， 其中圖3 （a）、 （b）、 （c） 分別代表壓入深度為0.3、0.6、 0.9 nm 時磨損表面形貌。 由圖3 （a）、 （b） 可知， 當(dāng)對偶件滑動速度為50 和100 m／s 時， 磨損表面較滑動速度為150 m／s 時更為粗糙。 由圖3 （c）可知， 當(dāng)對偶件壓入深度為0.9 nm 時， 滑動速度為150 m／s 時， 不銹鋼磨損表面形貌較滑動速度為50 和100 m／s 時更為粗糙。

圖3 壓入深度為0.3、 0.6、 0.9 nm 時不同滑動速度下表面磨損形貌Fig.3 Surface wear morphologies at different sliding velocity when the indentation depth is 0.3 nm （a）， 0.6 nm （b）， and 0.9 nm （c）

當(dāng)對偶件壓入S-07 不銹鋼然后發(fā)生相對滑動，在不銹鋼表面留下一道劃痕， 部分不銹鋼表面原子在力的作用下發(fā)生彈性變形和塑性變形， 被對偶件擠壓出不銹鋼表面。 一些不銹鋼原子在對偶件運動方向正前方和兩側(cè)堆積形成磨屑原子， 當(dāng)對偶件經(jīng)過某一位置后， 部分被擠壓的磨屑原子回到原來位置， 部分變形得到恢復(fù)。

為計算磨損體積損失， 導(dǎo)出S-07 不銹鋼在不同參數(shù)下摩擦磨損試驗后的模擬磨損軌跡， 使用OVITO 軟件的構(gòu)建表面網(wǎng)格（construct surface mesh） 對磨損軌跡進行表征， 選擇alpha-shape method 下的identify volumetric regions 選項工具對磨損區(qū)域進行識別和展示， 并計算磨損體積。

圖4 所示為滑動速度為50、 100、 150 m／s 時不同壓入深度下的磨損體積變化。 可知， 在滑動速度一定的情況下， 對偶件壓入深度越大， 不銹鋼的磨損量越大。 這是由于隨著對偶件壓入深度的增加， 在對偶件的作用下產(chǎn)生了更多的不銹鋼原子脫離原來的位置堆積形成磨屑， 伴隨對偶件一起滑動。 隨著滑動速度的增大， 體積磨損量大體上呈現(xiàn)出下降趨勢； 根據(jù)牛頓運動定律， 在相同的位移下， 靜止的原子獲得更大速度需要更大的驅(qū)動力， 隨著速度的增大， 部分原子未能獲得較大的速度， 無法與對偶件共同運動， 導(dǎo)致磨損量下降。

圖4 不同滑動速度和不同壓入深度下S-07 不銹鋼體積磨損量Fig.4 The wear volumeof S-07 stainless steel under different sliding velocity and indentation depth： （a） 50 m／s； （b） 100 m／s； （c） 150 m／s

3.3 位錯分析

當(dāng)S-07 不銹鋼與對偶件發(fā)生相對運動時， 在材料表面下方會發(fā)生亞表面損傷， 其損傷程度可以用統(tǒng)計分子動力學(xué)模擬中的位錯密度來表征， 位錯密度為位錯長度與材料體積的比。 研究中所使用的模型體積大小不變， 所以位錯長度所體現(xiàn)的變化趨勢與位錯密度變化趨勢一致。 采用OVITO 軟件中Dislocation Extraction Analysis （DXA） 模塊用于計算位錯線的長度。 圖5 所示為壓入深度為0.3、 0.6、 0.9 nm 時不同滑動速度下位錯線的長度變化。 可見摩擦滑動速度與位錯線長度大體上呈現(xiàn)負相關(guān)， 當(dāng)滑動速度較低時， 位錯線長度相對較長； 另外， 隨著相對滑動距離的增加， 大多數(shù)位錯線長度呈現(xiàn)上升趨勢， 預(yù)示隨著滑動距離的增加， 亞表面損傷程度增加， 內(nèi)部缺陷增多； 但壓入深度為0.6、 0.9 nm 和滑動速度為150 m／s時， 滑動距離為2～4 nm 時， 位錯線長度呈下降趨勢。

圖5 不同壓入深度和滑動速度下位錯線長度隨著相對滑動距離變化Fig.5 Variation of the length of dislocation line with relative sliding distance under different indentation depth and sliding velocity： （a） 0.3 nm； （b） 0.6 nm； （c） 0.9 nm

選壓入深度為0.6 nm、 滑動速度為150 m／s 作為例子進行了分析， 結(jié)果如圖6 所示。 從圖6 （a） 可知， 在滑動距離為2 nm 時， 存在5 條特殊的位錯線1／3＜1 0 0 ＞， 這種位錯線被稱為Hirth 位錯。 由于Hirth 位錯的滑移方向不在密排平面上， 很難滑動，因此它也被稱為Hirth 位錯鎖。 從圖6 （b）、 （c） 可知， 在滑動距離增加到4 nm 過程中， 位錯線長度下降。 分析圖7 可知， 在Hirth 位錯上方存在低應(yīng)力區(qū)，由于這個區(qū)域的存在， 降低了應(yīng)力， 因此不會通過產(chǎn)生大量新的位錯來釋放應(yīng)力； 原本存在的位錯逐漸向表面移動并消失， 導(dǎo)致位錯線長度下降； 在滑動距離為4 nm 之后， Hirth 位錯與高應(yīng)力區(qū)靠近， Hirth 位錯獲得能量并分解為Shockley 位錯。 位錯鎖消失后，應(yīng)力上升， 不銹鋼材料繼續(xù)產(chǎn)生位錯來釋放應(yīng)力， 因此位錯線長度又呈現(xiàn)上升趨勢。

圖6 不同滑動距離下位錯分布（h＝0.6 nm， v＝150 m／s）Fig.6 The dislocation distribution under different sliding distance（h＝0.6 nm， v＝150 m／s）： （a）2 nm； （b）3 nm； （c）4 nm

圖7 不同滑動距離下應(yīng)力分布（h＝0.6 nm， v＝150 m／s）Fig.7 The stress distribution under different sliding distance （h＝0.6 nm， v＝150 m／s）： （a） 2 nm； （b） 3 nm； （c） 4 nm

4 結(jié)論

通過分子動力學(xué)模擬研究S-07 不銹鋼滑動摩擦磨損過程， 探究在不同滑動速度、 壓入深度下S-07不銹鋼摩擦因數(shù)變化規(guī)律、 表面磨損形貌特征以及內(nèi)部位錯演變過程及相關(guān)機制。 主要結(jié)論如下：

（1） 隨著摩擦滑動速度的上升， 摩擦因數(shù)均增大， 這是由于在單位時間內(nèi)需要破壞的金屬鍵更多，導(dǎo)致摩擦因數(shù)變大。

（2） 磨損量大體上與滑動速度呈負相關(guān)， 原因在于隨滑動速度上升部分原子無法與對偶件一起協(xié)同滑動形成磨屑。

（3） 滑動摩擦速度與位錯線長度大體上呈負相關(guān)， 當(dāng)滑動速度較低時， 位錯線長度相對較長； 當(dāng)滑動距離為2～5 nm， 速度為150 m／s 時， 由于Hirth 位錯和低應(yīng)力區(qū)共同作用的原因， 導(dǎo)致位錯線長度出現(xiàn)先下降后上升。