陳 志，蔡 垚，顧燦鴻

(四川大學 化學工程學院，四川 成都 610065)

機械密封為端面密封，其依賴兩個固體表面間的接觸和相對運轉，以防止流體泄漏。常規(guī)設計中通常按照接觸端面的名義面積來計算其端面比壓。然而，機械加工表面在本質上都是粗糙的，動、靜環(huán)端面的實際接觸僅發(fā)生在部分微凸體上，導致實際接觸面積僅占名義接觸面積很小的一部分[1]。干摩擦是機械密封最惡劣的潤滑狀態(tài)，其接觸狀態(tài)直接決定了密封的性能與壽命[2]。張永振[3]指出在干摩擦狀態(tài)下，摩擦熱量與速度、接觸壓力存在冪函數的關系，外載荷和速度對端面的接觸、摩擦磨損行為的影響尤為突出。所以，研究外載荷和速度對機械密封干摩擦時端面接觸特性的影響對了解其密封機理和端面磨損具有重要意義。

目前，對相對運動摩擦端面的動態(tài)過程的研究仍然缺乏直接觀察和檢測評估的手段，因此，為了準確描述摩擦端面的接觸特性，需要建立合理的接觸模型，前人對此進行了大量研究。Majumdar等[4–5]首先基于2維W–M分形函數建立了彈塑性接觸分形模型，探究了總載荷與真實接觸面積的關系式；Wang等[6–7]在M–B模型的基礎上進行改進，得到了M–B分形接觸修正模型；孫見君[8]、魏龍[9]、丁雪興[10]等參照M–B模型建立了機械密封摩擦端面的分形模型，通過2維理論計算對密封端面泄漏、真實接觸面積、溫度分布進行了研究。基于分形理論建立的接觸模型在各個領域都得到了廣泛運用[11–16]。可見，利用分形理論可以實現對粗糙面的唯一表征，不受尺度的影響。

干摩擦的宏觀表現行為是眾多接觸摩擦學行為耦合作用的結果[3]。粗糙表面上的接觸壓力、熱應力、溫度場是相互影響的，因此，考慮微凸體間的相互作用及熱–力耦合對分析摩擦端面的接觸特性尤為重要。孫星星等[17]考慮密封動靜環(huán)的熱力變形，數值求解了密封環(huán)的溫度分布，但忽略了粗糙面微凸體的影響。房桂芳等[18]依據分形理論提出了機械密封端面摩擦熱的耦合計算方法。Li等[19]提出了考慮磨損和熱變形的接觸機械密封分形泄漏模型。黃健萌[20–21]、丁雪興[22]等建立了粗糙面/平面的1維滑移接觸模型，通過有限元模擬求解了滑動摩擦過程中微凸體變形及熱–力耦合對粗糙面應力、溫度變化的影響?？偟膩碚f，對于機械密封摩擦端面接觸的研究主要在理論計算方面，模擬計算也僅采用了直線滑動摩擦的方式；對于相互接觸的表面而言，直線運動與回轉運動的摩擦學特性是不同的[3]，因此采用3維回轉接觸模型進行仿真將更符合實際。

針對碳石墨與碳化鎢配對的機械密封，首先，作者基于W–M分形函數建立了模擬其摩擦副端面的3維粗糙實體與理想光滑剛體的轉動摩擦接觸模型；然后，在考慮熱–力耦合及微凸體間相互作用的基礎上，通過ABAQUS有限元分析軟件求解了其瞬時干摩擦時密封端面間的接觸特性，并討論不同外載荷、滑動速度對其特性的影響，并針對API682中串聯式機械密封的第2級接觸式干摩擦機械密封（抑制機械密封）進行了分析，以期為探索機械密封在干摩擦時的發(fā)熱、摩擦行為和密封設計提供借鑒。

1 計算模型的建立

1.1 粗糙表面及轉動接觸模型的建立

表征粗糙接觸界面的形貌是摩擦學的熱點，Yan等[23]的研究給出了能夠描述粗糙表面微觀特征的3維W–M函數：式中：Z(x，y)為3維粗糙表面形貌的高度；D和G分別為分形維數和分形粗糙度；L為取樣長度；γ為表面輪廓頻率密度的特征參數，通常取1.5；M為曲面褶皺的重疊數，值為10；φm,n為隨機相位；nmax為最大頻率，其值為：

式中，Ls為截止長度（近似取材料的原子間距）。

對于同一粗糙表面，Z(x，y)可由與尺度無關的分形參數D與G唯一確定，因此式（1）具有在不同長度尺度上描述粗糙表面的能力[23]，并且隨著D增大，G減小，表面粗糙度降低，趨于光滑[24–25]。本文中，D取2.35，G取2.45×10–9m，此時表面粗糙度Ra的大小為0.8 μm。

由于接觸問題屬于典型的邊界非線性問題，且粗糙表面的特征尺度在微米級，考慮到計算成本，研究中將微元體的模型尺寸也控制在微米級，其中，靜環(huán)模型：外徑DAi=260 μm，內徑DAo=380 μm，厚度hA=60 μm；動環(huán)模型：外徑DBi=250 μm，內徑DBo=390 μm，厚度hB=60 μm。通過增大轉速的數量級可以實現模型微元體的真實線速度V，根據文獻[26]粗糙表面輪廓曲線的局部與整體的自相似性可知，微元體的接觸行為可以反映實際端面之間的接觸行為。

機械密封中，動、靜環(huán)材料常采用硬軟配對，摩擦副為硬質合金（WC）對碳石墨（M106K），兩者的彈性模量相差1～2個數量級[22]，見表1。因此，研究中將靜環(huán)（實體A）與動環(huán)（實體B）的接觸簡化為分形粗糙表面與剛性理想光滑平面的接觸[9]；同時，由于密封環(huán)模型軸對稱，為了減少計算時間，靜環(huán)采用1/8周期模型，如圖1所示。其中：實體A中，A2為頂面，A1為底面（分形粗糙面），A3、A4為側面，A5為內圓柱面，A6為外圓柱面；實體B中，B1為頂面（接觸面），B2為底面，B3為內圓柱面，B4為外圓柱面；RP為參考點；徑向X、周向Y、軸向Z（實體A、B厚度方向）為坐標系。

表1 摩擦副材料參數[22]Tab. 1 Material parameters of friction pair[22]

圖1 分形接觸模型Fig. 1 Fractalcontact model

1.2 位移邊界條件

計算時設置兩個時間段：在第1個時間段內，將均布載荷由0線性增加到P，加載到實體A的A2面上，同時對其4個側面施加沿X、Y方向的位移約束。對于實體B，約束其Z方向的位移，同時將實體B用參考點RP（0，0，0）進行剛性約束（設置實體B為剛性體，同時可以通過RP的運動控制實體B的運動），通過約束RP的6個自由度使B保持靜止，讓兩實體建立穩(wěn)定的接觸關系。在第2個時間段內，保持A2面上載荷P不變，同時給定實體B的參考點RP繞Z軸的轉動速度ω，使兩實體間發(fā)生轉動摩擦。

1.3 模型假設

由于計算條件的限制，本文在計算中做出如下假設：

1）滑動摩擦過程遵循庫侖定律，在計算時間范圍內認為摩擦系數保持不變，不受外載荷、滑動速度變化的影響。

2）構成實體A、B的材料各向同性，且在瞬態(tài)計算中，由于溫升不大，假設材料的熱物性參數不會改變。

3）忽略材料磨損的影響，認為摩擦消耗的功全部轉化為熱能，并被摩擦副吸收，則局部熱流密度可表示為：

式中，q為熱流密度，f為摩擦系數，Pb為粗糙表面微元體的接觸壓力，V為相對滑動速度，t為時間，r為半徑。

4）考慮摩擦熱在兩實體間的分配，認為接觸區(qū)域為理想熱傳導，即接觸點處兩實體的瞬時溫度相等，這樣摩擦熱流密度就能根據材料的熱物性參數及散熱環(huán)境在兩實體間自由分配[20]，同時忽略熱輻射的影響。

1.4 熱邊界條件

實體A、B的熱傳導方程可表示為：

式中，ρ為實體A、B材料的密度，c為比熱容，k為熱導率，T為溫度，t為時間，r為半徑。

初始條件：TA=TB=298 K，t=0。

耦合條件（假設4））：TA=TB，qA+qB=q（實體A、B接觸點）。

接觸面（A1、B1）的熱邊界條件為：忽略熱輻射，非接觸區(qū)域的傳熱為對流換熱；接觸區(qū)域為摩擦熱流的輸入[21]：

式中，接觸點處g(m)=1，非接觸點處g(m)=0。

由于本文研究的是瞬態(tài)問題，在很短的時間內，熱量只在兩實體接觸面的表層及亞表層傳遞，因此可以認為密封環(huán)非接觸面（A2、B2）是絕熱的；同樣地，干摩擦條件下可以認為實體A及實體B的側面A3～A6、B3～B4是絕熱的。

1.5 求解方法

采用ABAQUS/Explicit求解計算模型，對于復雜的接觸問題，顯示方法求解效率較高且結果可靠。其采用中心差分的方法求解運動方程，用一個時間段的運動學條件計算下一個時間段的運動學條件，通過自動劃分足夠小的時間增量段來確保模型狀態(tài)穩(wěn)定，從而獲得高精度的解[27]。求解前需要設置接觸對，將剛性光滑面B1設置為主面，分形粗糙面A1設為從面，采用罰函數法作為求解接觸問題的算法，兩表面間的相對滑動通過有限滑移公式描述。

考慮熱–力耦合計算，使用C3D8T六面體熱–力耦合單元對整個計算模型進行網格劃分，并對實體A的粗糙面（A1）的表層及亞表層進行網格細化，對光滑實體B使用較粗的網格劃分方式。為了證明求解結果與網格大小無關，對模型進行網格無關性驗證分析，由于主要研究對象是粗糙表面，所以采用調整粗糙表面網格大小的方式驗證接觸面積的變化。定義無量綱接觸面積為Ar/Aa，其中，Ar為真實接觸面積，Aa為粗糙面A1的名義接觸面積。ABAQUS軟件中可設置接觸面積的監(jiān)控，計算完成后輸出Ar隨時間的變化，則無量綱接觸面積隨網格尺寸的變化如圖2（a）所示?？梢园l(fā)現，隨著網格尺寸減小，無量綱接觸面積趨于穩(wěn)定，實際粗糙面網格尺寸為1.736 μm，可以較好地保證計算結果的準確性。網格劃分完成后，光滑實體單元數目為10 400，粗糙實體單元數為20 475，如圖2（b）、（c）所示。

圖2 網格劃分Fig. 2 Mesh generation

2 算例分析

2.1 計算參數選擇

動環(huán)、靜環(huán)的材料如表1所示，兩實體的干摩擦系數選擇為f=0.25[28]，靜環(huán)材料的屈服強度為σy=200 MPa[9]。在初期計算中，外載荷P過小時，數據變化不明顯，故P選取變化范圍為0.2～40.0 MPa；滑動速度V為光滑實體B外徑處的線速度，變化范圍為20～50 m/s。

在模擬計算過程中不考慮材料非線性問題，只考慮幾何非線性與接觸邊界非線性，由于計算過程難以收斂，同時粗糙表面實際處于不斷磨損的狀態(tài)中，表面形貌也不斷更新[24]，因此，只求解轉動發(fā)生后60 μs以內的瞬態(tài)過程。為了使兩實體發(fā)生平穩(wěn)接觸，將轉動前的計算時間設置為40 μs，外載荷以線性增加的方式施加在A2面上直到載荷等于P為止，之后保持載荷不變使兩實體發(fā)生相對轉動，速度為V；密封環(huán)初始溫度與環(huán)境溫度設置為298 K。每次計算改變外載荷P、滑動速度V兩個參數中的一個值，計算得到密封環(huán)不同工況下的7組摩擦數據。

2.2 結果分析與討論

2.2.1 真實接觸面積

定義無量綱載荷為P/E′，其中：

得到本文分形模型與魏龍分形模型[9]、M–B分形模型[5]理論值的對比如圖3所示。

圖3 本文分形模型與其他模型理論值的對比Fig. 3 Comparison of the theoretical values of fractal model in this paper with other models

從圖3中可以發(fā)現，隨著外載荷增大，真實接觸面積近似呈線性增大，外載荷增加3倍時，接觸面積增大了2.33倍，這是因為外載荷的增加導致粗糙表面受力增加，接觸面主峰微凸體發(fā)生彈塑性變形，次峰微凸體逐漸參與接觸，增加了接觸面積。同時，可以看出：隨著載荷的增大，模擬分析的結果與理論值之間具有相同的變化趨勢，在載荷較低時，它們間差異較小（串聯式機械密封第二級接觸式干運轉機械密封恰恰是在低載荷下工作）；隨著載荷的進一步增大，它們間的誤差也逐漸增大（無量綱載荷為0.001 5時，M–B模型計算得到的無量綱接觸面積是本文模擬計算的1.364倍，是魏龍模型的1.186倍）。產生這一差異的原因在于：魏龍分形模型、M–B分形模型的建立是基于2維粗糙輪廓的，而不是真實的3維粗糙表面，它們未考慮粗糙表面相鄰微凸體之間的相互作用。文獻[29]的研究表明，魏龍分形模型、M–B分形模型計算得到的真實接觸面積與實驗值相比較大。本文所采用的3維有限元模擬考慮了粗糙表面相鄰微凸體之間的相互作用，所以其結果更接近實際，從而驗證了模擬分析的正確性。

圖4為P=30 MPa時，不同滑動速度下，滑動0.3 μs及60.0 μs時粗糙面的無量綱接觸面積。可以發(fā)現，不同滑動時間下，接觸面積隨速度增加的斜率相同，速度越大，接觸面積越大。但在滑動速度增加1.5倍的情況下，接觸面積僅增大為原來的1.002倍，因此，粗糙面的真實接觸面積的大小主要取決于外載荷，而滑動速度影響較小。

圖4 滑動階段無量綱接觸面積隨滑動速度的變化Fig. 4 Non-dimensional contact area vs sliding velocity when ring sliding

2.2.2 接觸壓力

摩擦磨損是表面行為，因此，研究粗糙表面的真實接觸面積、接觸壓力、溫度的分布極其重要。圖5（a）～（d）所示為V=30 m/s時，P=10、20、30、40 MPa下，滑動60 μs時刻分形粗糙面的接觸壓力分布云圖。由圖5可以看出：隨著外載荷增大，真實接觸面積Ar上較高的接觸壓力節(jié)點在增加，粗糙表面上最大接觸壓力值位于節(jié)點13538；當外載荷為10、20、30、40 MPa時，節(jié)點13538的接觸壓力分別為415.88、416.64、423.25和431.59 MPa。對比不同載荷下各個接觸區(qū)域上接觸壓力的變化，可以發(fā)現外載荷增大3倍后，節(jié)點13538上接觸壓力值增加的幅度并不是很大，這說明相對運動的摩擦副的接觸壓力峰值達到一定值后變化緩慢，該值對外載荷的變化不敏感，應與摩擦副的材料力學特性有關。

粗糙表面節(jié)點13538的接觸壓力隨載荷增加變化緩慢，其原因是接觸壓力在達到一定值后，材料發(fā)生彈塑性變形，次峰微凸體進入接觸并承擔了額外的載荷，使得節(jié)點13538所在微凸體受力得到緩解。表現為粗糙表面的接觸壓力峰值維持在一個定值范圍內，而其他接觸節(jié)點的接觸壓力繼續(xù)增加。換句話說，在本文的研究參數范圍內，接觸壓力峰值是“自限性”的，這種特征保證端面在較大的外載荷下的接觸是由點接觸到局部面接觸、一個節(jié)點接觸到眾多節(jié)點接觸。更多的節(jié)點參與接觸，端面受力條件得到改善，這種“自限性”對于密封環(huán)端面微凸體受力狀況的改善是有利的。

圖5 不同外載荷下滑動60 μs時粗糙表面接觸壓力PB云圖Fig. 5 Contourmap of contact pressure on the rough surface when sliding at 60 μs under different external loads

圖6（a）～（d）所示為P=30 MPa時，V=20、30、40、50 m/s下，滑動60 μs時刻分形粗糙面的接觸壓力分布云圖。

由圖6可以發(fā)現，在外載荷相同時，接觸壓力的分布基本相同，最大接觸壓力分別為430.34、423.25、428.86、420.26 MPa，均位于粗糙表面節(jié)點13538，并隨著速度增加具有微小的下降趨勢，這是因為滑動期間速度增大時引起的微凸體局部摩擦熱流密度變大，微凸體熱變形更大使接觸節(jié)點局部接觸面積略微增大，導致節(jié)點最大接觸壓力略有減小。

圖6 不同速度下滑動60 μs時粗糙表面接觸壓力Pb云圖Fig. 6 Contourmap of contact pressure on rough surface when sliding at 60 μs under different sliding velocity

2.2.3 溫度分布

圖7（a）～（d）所示為V=30 m/s，P=10、20、30、40 MPa下，滑動60 μs時刻分形粗糙面的溫度分布云圖。由于粗糙表面微凸體的形狀大小、高低不一，滑動接觸正是發(fā)生在這些不連續(xù)的接觸微凸體上。摩擦所引起的熱量產生、釋放也發(fā)生這些區(qū)域，使得接觸節(jié)點處溫度較周圍較高，形成“熱點”。

由圖7可以看出，P=10 MPa時高溫區(qū)主要分布在粗糙表面微凸體接觸點位置，且端面上存在明顯的溫度梯度分布，會有熱應力產生。隨著載荷增大，粗糙表面接觸微凸體及接觸點越多，在滑動瞬間形成的“熱點”也越多，接觸中心區(qū)域的溫度也越高，粗糙表面整體溫度也增高。另外，可以發(fā)現，不同載荷下最高溫度值所在的節(jié)點不同，不是單純地位于接觸壓力最大的位置。這表明粗糙表面的最高溫度還受到局部熱流密度的影響；即最高溫度受到局部接觸壓力和局部滑動速度和局部熱傳遞的共同影響。

圖7 不同外載荷下滑動60 μs時粗糙表面溫度云圖Fig. 7 Contourmap of temperature on rough surface when sliding at 60 μs under different external loads

隨著外載荷的增加，雖然粗糙表面最大接觸壓力改變較小，但載荷變化使得粗糙表面接觸點的數量明顯增多，滑動過程中粗糙表面上的摩擦熱源增多，導致了端面整體溫度升高，故材料的耐摩擦磨損性能也受到影響。

圖8所示為外載荷為10 MPa、滑動速度為30 m/s，摩擦滑動300 μs條件下的粗糙面各節(jié)點的溫度變化曲線（節(jié)點分布如圖7（a）所示）。由圖8可以看出：0～40 μs為外載荷加載階段，粗糙面各節(jié)點溫度保持不變；當兩表面開始滑動，40～42 μs內由于速度突增及熱傳導的滯后性出現了溫度急升；隨后42～340 μs，在自身熱傳導的作用下，各節(jié)點溫度開始緩慢上升，溫度上升斜率也是一致的，且滑動60 μs的溫升斜率與滑動300 μs的溫升斜率也相同。這也驗證了本文采用摩擦滑動60 μs的計算時間進行研究是可行的。

圖8 摩擦滑動300 μs內粗糙面節(jié)點溫度隨時間的變化Fig. 8 Temperature change of rough surface node with time within 300 μs of friction sliding

API682中規(guī)定串聯式機械密封第2級接觸式干運轉機械密封的介質壓力應控制在0.07 MPa以內。為了闡釋該密封在有氣體冷卻的配合下實現25 000 h運轉的可能性，本文模擬計算了外載荷0.2 MPa、滑動速度30 m/s條件下滑動60 μs時粗糙面溫度分布云圖，如圖9所示。由圖9可以看到，摩擦滑動60 μs后，粗糙面的最高溫度僅從環(huán)境溫度25 ℃升高到25.42 ℃，溫升僅有0.42 ℃，粗糙面溫度變化不明顯。圖10為滑動速度30 m/s時，不同外載荷下粗糙面最高溫度的變化，隨著外載荷的增加，最高溫度是線性增大的，與式（3）描述的規(guī)律相符?？梢娊档蜋C械密封的外載荷可以有效防止其干摩擦時石墨環(huán)溫升過高，API682標準中建議介質壓力≤0.07 MPa是有必要的。

圖9 0.2 MPa下摩擦滑動60 μs時粗糙表面溫度云圖Fig. 9 Contourmap of temperature on rough surface when sliding at 60 μs under 0.2 MPa

圖10 不同外載荷下滑動60 μs時粗糙面最高溫度的變化Fig. 10 Change of maximum temperature on rough surface under different load when sliding at 60 μs

圖11（a）～（d）所示為P=30 MPa，V=20、30、40、50 m/s下，滑動60 μs時刻分形粗糙面的溫度分布云圖。由圖11可以發(fā)現，在相同的外載荷下，對應的粗糙表面最高溫度值分別為44.84、49.12、52.99和57.52 ℃，隨著滑動速度增加溫度呈上升趨勢，粗糙表面的溫度不均勻趨勢明顯加劇，各區(qū)域溫度值均有所提高。根據式（3），速度增加，摩擦發(fā)熱增加，粗糙表面的溫度梯度也隨之變大，熱應力也會加劇。

圖11 不同速度下滑動60 μs時粗糙表面溫度云圖Fig. 11 Contourmap of temperature on rough surface when sliding at 60 μs under different sliding velocity

3 結 論

在P=0.2～40 MPa、V=20～50 m/s，f=0.25和粗糙表面Ra=0.8 μm的工況下，針對碳石墨–硬質合金配對摩擦副在干摩擦狀態(tài)下的接觸特性進行了分析，研究發(fā)現：

1）粗糙端面真實接觸面積的大小主要取決于外載荷。外載荷增加真實接觸面積近似線性增大，但滑動速度影響較小。

2）隨著外載荷與滑動速度的增加，石墨粗糙表面微凸體接觸壓力峰值在415～432 MPa范圍內，表現出一種“自限性”。

3）粗糙端面溫度分布不均勻，局部溫度較高；端面最高溫度并非位于接觸壓力最大的位置，而是受到局部接觸壓力和局部滑動速度及局部熱傳遞的共同影響。

4）外載荷、速度增大均會使端面整體溫度增加。

5）API682標準中要求抑制機械密封的介質壓力≤0.07 MPa是有必要的。降低介質壓力可以降低外載荷，從而可以有效防止干摩擦時密封環(huán)溫升過高。