盤式制動器熱-應力-磨損耦合行為的數值模擬*

張方宇，桂良進，范子杰

(清華大學，汽車安全與節能國家重點實驗室，北京 100084)

前言

盤式制動器的制動性能主要取決于摩擦片和制動盤之間的摩擦狀況。制動過程中，摩擦力做功生熱，摩擦片被磨損，熱、接觸應力和磨損三者相互耦合，共同影響制動器的制動性能。

涉及磨損的理論有很多，如文獻[1]中的線性磨損理論、文獻[2]中提出的磨粒磨損微切削理論、文獻[3]中提出的粘著理論、文獻[4]中提出的疲勞磨損理論、文獻[5]中提出的能量磨損理論、文獻[6]中提出的剝層磨損理論、文獻[7]中基于“摩擦疲勞”概念發展的相應的理論和文獻[8]中提出的磨損率與磨損循環次數的冪指數關系法則。其中文獻[1]中的線性磨損理論應用最為廣泛。

目前對于制動熱應力和磨損行為的研究主要采用的是有限單元法。文獻[9]中研究了飛機正常著陸時，制動盤在制動過程中的熱應力分布，以及熱應力對C/C復合材料磨損表面形貌的影響，并利用有限元軟件進行了熱應力的仿真，但未對磨損行為進行有限元的模擬。文獻[10]中利用有限元仿真軟件Abaqus對織物襯墊自潤滑關節軸承摩擦副間的應力分布進行了仿真，并根據Archard的磨損理論，計算了單步磨損量，以此循環，利用顯式積分的方法計算總磨損量并描繪了應力分布情況的演變。文獻[11]中使用了類似的方法對盤式制動器磨損行為進行了研究，其中還使用了模型結構矩陣(MSM)[12]的模塊化建模方法，但它們在計算磨損時均未考慮摩擦熱的影響。總之，當前對盤式制動器進行熱-應力-磨損耦合行為的數值模擬，尚無文獻報到。

本文中探索了盤式制動器熱-應力-磨損耦合行為的數值模擬方法。首先建立了盤式制動器的有限元模型；然后利用Archard的磨損理論實現了摩擦片接觸應力-磨損耦合分析的數值計算方法，并模擬了摩擦片和制動盤的磨合過程；最后探索了模型在熱-應力耦合條件下摩擦片的磨損行為，并對汽車100km時速緊急制動工況下的摩擦片磨損過程進行了仿真。

1 有限元模型

盤式制動器組件的結構示意圖如圖1所示。本文中只關心摩擦片和制動盤接觸界面的接觸應力和磨損的分布以及摩擦熱和熱傳導的過程，因此只對活塞、背板、摩擦片和制動盤進行有限元建模。制動盤兩側的制動襯片在幾何形狀、支撐形式和受力上雖有差異，但在與制動盤的作用形式方面還是比較接近，可近似地認為它們關于制動盤中面對稱，因此可以只對一個制動襯片和半個制動盤建模，簡化后的結構如圖2所示。另外，由于活塞和制動盤材料較硬，制動過程中的變形對仿真結果影響很小，其磨損率相比摩擦片小很多，可以忽略，因此兩者均可以當作剛體處理，只將制動襯片(包括背板和摩擦片)視為變形體來建模，且忽略制動盤的磨損，只考慮摩擦片的磨損。

整個有限元模型的爆炸圖如圖3所示。其中，活塞是解析剛體；背板是三維變形體，共包含1 467個節點，876個單元，單元類型為C3D8T；摩擦片是三維變形體，包含2 420個節點，1 662個單元，單元類型為C3D8T；制動盤是離散剛體，包含1 500個節點，960個單元，單元類型為C3D8T。背板的材料是普通鋼鐵，使用線彈性材料特性，彈性模量E=2.09×109Pa，泊松比μ=0.3，其熱力學材料性能參數來自通用鋼材的材料特性數據，其材料性能見表1。

熱導率/(W/(m·K))密度/(kg/m3)彈性模量/GPa泊松比膨脹系數/K-1比熱容/(J/(kg·K))4878002090311×10-5452

摩擦片的材料特性極其復雜，此處假定其力學性能為線彈性，彈性模量和泊松比參考文獻[11]，熱力學材料性能參數參考文獻[13]。本文中邊界條件、載荷和接觸統稱為接口(interface)。制動盤只有繞其軸線的轉動自由度(i2)，其他各個方向的自由度全部為零(i1)。背板受到制動鉗的約束，只能沿z軸方向移動，沿x和y軸被約束(i3)。摩擦片嵌合在背板上(i7)。活塞只能沿z軸移動(i4)，其參考點受到沿z軸負方向的集中載荷作用(i5)，提供制動力。摩擦片和制動盤之間為有摩擦接觸約束(i6)并且包含摩擦生熱特性，生熱效率為100%，分配系數為0.5，接觸時相互之間的熱傳導系數為1×109W/(m2·℃)。背板和活塞之間為無摩擦接觸約束(i8)。

在Abaqus中，模型的散熱邊界條件通過背板、摩擦片和制動盤三者表面的surface film condition實現，其散熱方程為

q=-h(θ-θ0)

(1)

式中：q為邊界熱流密度；h為膜層系數；θ為表面節點溫度；θ0為冷源溫度。

計算時，設定θ0=20℃，制動盤和背板的膜層系數設為100，摩擦片的膜層系數設為10。

2 應力-磨損耦合分析方法

在制動過程中，摩擦片和制動盤的接觸伴隨著摩擦和磨損，其程度受眾多因素的影響，比如接觸應力、滑動速率、材料特性等等。有很多磨損理論揭示了這種關系，其中Archard的摩擦定律相對簡單，使用較為廣泛，并且工程應用比較成熟，因此本文中也采用此理論模型。

Archard的磨損模型通常具有以下形式：

(2)

式中：V為磨損體積；s為磨損位移；k為無量綱的磨損系數；FN為法向載荷；H為磨損表面的硬度。

在無窮小的接觸面積ΔA，無窮小時間增量dt下，將此形式的磨損模型推廣到微分形式：

(3)

式中： dV為磨損體積增量；ds為滑動位移；FN為法向接觸力。

假設ΔA在dt下的磨損厚度增量為dh，則dV=ΔAdh，代入到式(3)得

(4)

式中：FN/ΔA為接觸微面積ΔA處的接觸應力，假設其為pc；k/H為有量綱磨損系數，可以用kd表示，此時式(4)變為

dh=kdpcds

(5)

如果取有限小的時間增量為Δt，相應的磨損厚度增量為Δh，相對滑動位移增量為Δs，kd和pc假設在Δt內變化不大，仍取瞬時值，則式(5)改寫為

Δh=kdpcΔs

(6)

在盤式制動器制動過程中，可以認為摩擦片是靜止不動的，Δs為制動盤相對摩擦片滑動位移，Δh為Δt時間內的磨損深度。

不考慮熱效應時，kd可認為是常數，Δh是pc的函數，但Δh在摩擦片的各個部位并不相同，從而對接觸應力pc又產生影響，所以磨損過程是Δh和pc耦合的過程。為了求解這個耦合過程，本文中采用顯式歐拉積分的方法。首先利用有限元軟件計算得到接觸應力pc。當Δs較小時，可近似認為pc為不變量，從而求得對應的Δh，再根據h變化之后的模型，應用有限元軟件求得變化之后模型的接觸應力pc，如此循環計算下去，求解最終的磨損過程。具體計算過程如下：

假設摩擦片和制動盤之間的接觸點為i=1,…,N，顯式歐拉積分的步數為j=1,…,M,第i個接觸點第j個積分步的接觸應力為pc(i,j)，磨損深度為Δh(i,j)，j次總磨損深度為h(i,j)，制動盤轉動位移為Δs(i,j)。則有如下關系：

Δh(i,j)=kdpc(i,j)Δs(i,j)

(7)

h(i,j+1)=h(i,j)+Δh(i,j)

(8)

利用Abaqus Standard的ALE功能來實現磨損的仿真，其中磨損的數值算法通過用戶子例程UMESHMOTION實現。ALE算法中設置磨損邊界的約束類型為速度約束，將式(7)左右兩側同時除以第j個時間步長Δt(j)，得到如下關系：

(9)

(10)

3 摩擦片-制動盤界面磨合仿真

在全新的制動器進入穩定磨損之前，摩擦片-制動盤接觸界面先要經歷磨合過程。由于整個磨合期滑動位移和接觸時間很長，現有的普通計算機的計算能力幾乎無法模擬整個過程中溫度場、接觸應力場和磨損三者的耦合行為，因此在摩擦片-制動盤接觸界面磨合過程的仿真中忽略了溫度場的影響，只對接觸應力-磨損兩個場的耦合過程進行模擬。

磨合期摩擦片和制動盤之間的摩擦磨損可以看作一個準靜態過程，其工況如圖4(a)所示，對制動器的活塞施加沿制動盤軸線方向1.8kN載荷(相當于1MPa制動液壓力)，制動盤做勻速運動。

此工況的仿真在Abaqus中采用Static General 算法來實現，計算用時約3h，計算得到的摩擦片接觸應力分布和累積磨損量分布隨制動盤轉數的演變過程如圖5和圖6所示。

由圖5和圖6可見：

(1)制動盤剛開始轉動(0轉)時，摩擦片-制動盤接觸界面在受活塞底部環面擠壓的部位接觸應力比較大，其他部位接觸應力相對較小。但隨制動盤轉數的增加，應力大的區域磨損較大，使接觸應力的分布又趨于均勻；

(2)制動盤剛開始轉動(0轉)時，摩擦片-制動盤接觸界面接觸應力的分布左右并不對稱，這是由摩擦力的存在造成的。制動盤從摩擦片的左側向右側滑過，它們之間的摩擦力使摩擦片的左側壓緊制動盤，右側脫離接觸，因此右側的接觸應力為零；但是隨著制動盤轉數的增加，接觸面中間應力較大，磨損較快，而右側沒有接觸，也就沒有磨損，這就導致摩擦片右側逐漸進入接觸，接觸應力由小變大；

(3)當制動盤轉數達到一定值后(10 000～20 000轉)，摩擦片接觸面上的接觸應力分布達到穩定狀態，幾乎不再變化，其應力值沿制動盤半徑方向逐漸減小，即從摩擦片內側向外側逐漸減小。這是由于在Archard磨損模型里磨損量與滑動距離和接觸應力的乘積成正比，當磨損達到穩定時，摩擦片的整個接觸面均勻地向背板方向磨損，也就是說接觸面上各點的磨損速率相等，因此接觸應力與滑動距離成反比。又因為滑動距離和制動盤上相應點的半徑成正比，所以接觸應力和半徑成反比，這也就是接觸面應力沿半徑方向均勻減小的原因。

以上3點分析在一定程度上可以驗證利用Abaqus中ALE算法模擬磨損過程的正確性。

4 摩擦片-制動盤界面緊急制動仿真

在制動盤轉動20 000轉之后，摩擦片和制動盤結束了磨合期，進入穩定磨損階段。本節在上一節摩擦片和制動盤完成磨合期(20 000轉)后的摩擦片輪廓基礎上，研究一次緊急制動過程中熱-應力-磨損三者的耦合行為。一次緊急制動工況如圖4(b)所示，對制動盤的活塞施加沿制動盤軸線方向9kN載荷(相當于5MPa制動液壓力)，制動盤做初速度為120rad/s(相當于汽車100km/h)的勻減速轉動，3s后停止。

本文中利用Abaqus中Coupled temp-displacement (Transient) 算法來實現此工況，同時模擬了考慮和不考慮磨損兩種情況，并對結果做出了比較。

考慮磨損時，計算用時約16h，計算得到的摩擦片接觸應力分布、溫度分布、累積磨損量分布和制動盤溫度分布隨時間的演變過程如圖7～圖10所示。

對圖7～圖10展示的結果做如下分析。

(1)一次緊急制動的制動力較大，所以開始制動 (t=0.001s) 時，摩擦片和制動盤接觸界面在受活塞環擠壓的部位接觸應力較大，但隨著摩擦力做功，接觸界面溫度升高，摩擦片向制動盤方向膨脹，這時摩擦片表面的中心部位與制動盤接觸，接觸應力的分布變為中心大，四周小。

(2)接觸應力大的部位摩擦功率大，溫度升高快，因此摩擦片溫度分布和接觸應力的分布大體一致。但是考慮到摩擦片與制動盤和空氣的熱交換，以及摩擦片內部的熱傳導都會使摩擦片的溫度分布趨于均勻，因此摩擦片的溫度分布又與應力分布不完全一致。

(3)剛開始制動時，制動盤轉速快，摩擦功率大，摩擦片的溫度迅速上升。但在制動過程中隨著制動盤轉速逐漸減小，摩擦功率也會減小，而隨著溫度的上升，散熱功率增加，因此摩擦片的溫升越來越慢，到1.5s左右摩擦片和制動盤的溫度達到最高值，隨后摩擦功率小于散熱功率，摩擦片和制動盤的溫度開始下降并趨于均勻。

以上3點分析在一定程度上說明本次模擬比較符合實際經驗，具有準確性。

不考慮摩擦片磨損時，計算用時約12h，計算得到的摩擦片接觸應力分布、溫度分布和制動盤溫度分布隨時間的演變過程與考慮磨損的情況類似。

選取摩擦片接觸應力最大值附近一點，節點編號為N2588，如圖11(a)所示。在考慮磨損和不考慮磨損兩種情況下將該點接觸應力隨時間變化的曲線進行比較，如圖12(a)所示。由圖可見，考慮磨損時接觸應力值比不考慮磨損時的接觸應力值小，其原因是接觸應力大的地方磨損量大，使接觸應力的分布相對均勻。

選取摩擦片溫度最高值附近一點，節點編號為N2602，如圖11(b)所示。在考慮磨損和不考慮磨損兩種情況下將該點溫度隨時間變化的曲線進行比較，如圖12(b)所示。從圖中可以看到，考慮磨損時溫度值比不考慮磨損時的溫度值略小，其原因是考慮磨損時，接觸應力相對均勻，摩擦生熱的熱流密度也相對均勻。但這種差別不大，原因是一次緊急制動磨損量很小，對溫度分布的影響不顯著。

5 結論

本文中探索了盤式制動器熱-接觸應力-磨損耦合行為的數值模擬方法。首先建立了盤式制動器的有限元模型；然后模擬了盤式制動器磨合期摩擦片的接觸應力-磨損耦合過程；最后模擬了模型在一次緊急制動過程中熱-接觸應力-磨損三者的耦合情況，并將結果與不考慮磨損情況下的結果進行了比較。分析結果表明，所提出的模擬方法是有效的，可為盤式制動器的結構參數設計提供理論指導。

[1] Archard J F. Contact and Rubbing of Flat Surfaces [J]. Journal of Applied Physics，1953, 24: 981-988.

[2] 赫羅紹夫M M, 巴比契夫M A. 金屬的磨損 [M]. 胡紹衣, 佘滬生，譯.北京: 機械工業出版社, 1966.

[3] Bowden F P,Tabor D.The Friction and Lubrication of Solids [M].Oxford at the Clarenden Press, 1964.

[4] 克拉蓋爾斯基иB，陀貝欽M H, 康巴洛夫B C. 摩擦磨損計算原理 [M]. 汪一麟,朱安仁，范明德，譯.北京: 機械工業出版社, 1982.

[5] Fleisher G. Energishe Methode der Bestimung der Verschleiβes [J]. Schmierungstechnik, 1973(4): 9.

[6] Suh N P. An Overview of the Delamination Theory of Wear [J]. Wear, 1977, 44(1):1-16.

[7] СосновскийЛ А.Основы Трибофатики [M]. Гомель, 2003.

[8] Siniawski M T, Harris S J, Wang Q. A Universal Wear Law for Abrasion [J]. Wear, 2007, 262: 883-888.

[9] 徐惠娟, 易茂中, 熊翔，等. 制動過程中熱應力對C/C復合材料磨損表面形貌的影響 [J]. 中國有色金屬學報, 2011, 21(1): 131-137.

[10] Shen X,Cao L, Li R. Numerical Simulation of Sliding Wear Based on Archard Model[C]. Mechanic Automation and Control Engineering (MACE)，2010：325-329.

[11] Soderberg A, Andersson S. Simulation of Wear and Contact Pressure Distribution at the Pad-to-rotor Interface in a Disc Brake Using General Purpose Finite Element Analysis Software [J]. Wear, 2009, 267: 2243-2251.

[12] Sellgren U. Architecting Models of Technical Systems for Non-routine Simulations[C]. Proceedings of ICED 03-14th International Conference on Engineering Design, Stockholm, 19-21 August 2003.

[13] Day A J, Newcomb T P. Dissipation of Frictional Energy from the Interface of an Annular Disc Brake [C]// Proceedings of the Institution of Mechanical Engineers, Part D: Transport Engineering, 1984, 198(11): 201-209.