端面密封材料S-07 不銹鋼滑動摩擦學行為的分子動力學模擬*

王 權 莊宿國 黃 丹 朱正興 劉秀波

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

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

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

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

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

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

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

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

1 模型和勢函數

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

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

（2） 對建立好的模型進行幾何結構優化以及能量最小化處理；

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

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

經過上述過程， 建立的目標S-07 不銹鋼的結構模型， 如圖1 所示。

圖1 S-07 不銹鋼分子動力學模擬結構模型Fig.1 S-07 stainless steel molecular dynamics simulation structural model

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

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

2 數值模擬

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

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

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

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

3 實驗結果與分析

3.1 摩擦因數分析

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

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

圖2 不同壓入深度下對偶件的平均摩擦因數隨相對速度的變化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 壓入深度下， 隨著相對滑動速度增加， 摩擦因數均呈現上升趨勢。 在S-07 不銹鋼與對偶件的相對滑動過程中， 對偶件先壓入表面， 然后與不銹鋼發生相對滑動； 期間， 對偶件與不銹鋼表面存在剪切作用， 不銹鋼表面發生彈性和塑性變形。 當作用在不銹鋼表面的機械能足以破壞S-07 不銹鋼原子間的金屬鍵后， 不銹鋼原子就會隨對偶件一起滑動， 進而被磨損去除形成磨屑。 摩擦因數增加是因為在單位時間內， 所需要破壞的金屬鍵更多， 從而表現為摩擦因數上升。

3.2 磨損形貌分析

不同參數下的模擬磨損表面形貌如圖3 所示， 其中圖3 （a）、 （b）、 （c） 分別代表壓入深度為0.3、0.6、 0.9 nm 時磨損表面形貌。 由圖3 （a）、 （b） 可知， 當對偶件滑動速度為50 和100 m／s 時， 磨損表面較滑動速度為150 m／s 時更為粗糙。 由圖3 （c）可知， 當對偶件壓入深度為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）

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

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

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

圖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 位錯分析

當S-07 不銹鋼與對偶件發生相對運動時， 在材料表面下方會發生亞表面損傷， 其損傷程度可以用統計分子動力學模擬中的位錯密度來表征， 位錯密度為位錯長度與材料體積的比。 研究中所使用的模型體積大小不變， 所以位錯長度所體現的變化趨勢與位錯密度變化趨勢一致。 采用OVITO 軟件中Dislocation Extraction Analysis （DXA） 模塊用于計算位錯線的長度。 圖5 所示為壓入深度為0.3、 0.6、 0.9 nm 時不同滑動速度下位錯線的長度變化。 可見摩擦滑動速度與位錯線長度大體上呈現負相關， 當滑動速度較低時， 位錯線長度相對較長； 另外， 隨著相對滑動距離的增加， 大多數位錯線長度呈現上升趨勢， 預示隨著滑動距離的增加， 亞表面損傷程度增加， 內部缺陷增多； 但壓入深度為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 作為例子進行了分析， 結果如圖6 所示。 從圖6 （a） 可知， 在滑動距離為2 nm 時， 存在5 條特殊的位錯線1／3＜1 0 0 ＞， 這種位錯線被稱為Hirth 位錯。 由于Hirth 位錯的滑移方向不在密排平面上， 很難滑動，因此它也被稱為Hirth 位錯鎖。 從圖6 （b）、 （c） 可知， 在滑動距離增加到4 nm 過程中， 位錯線長度下降。 分析圖7 可知， 在Hirth 位錯上方存在低應力區，由于這個區域的存在， 降低了應力， 因此不會通過產生大量新的位錯來釋放應力； 原本存在的位錯逐漸向表面移動并消失， 導致位錯線長度下降； 在滑動距離為4 nm 之后， Hirth 位錯與高應力區靠近， Hirth 位錯獲得能量并分解為Shockley 位錯。 位錯鎖消失后，應力上升， 不銹鋼材料繼續產生位錯來釋放應力， 因此位錯線長度又呈現上升趨勢。

圖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 不同滑動距離下應力分布（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 結論

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

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

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

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