崔 躍， 高 飛， 符 蓉， 李陽杰

(大連交通大學(xué) 連續(xù)擠壓教育部工程研究中心， 遼寧大連 116028)

對于制動摩擦副在制動過程的溫度變化規(guī)律，利用有限元計算是一種重要的研究手段。有限元法就是近似求解一個物理現(xiàn)象的偏微分方程，而有限元法又是求解這些偏微分方程的一種手段，即借助離散網(wǎng)格來近似求解每個單元所代表的線性方程，最后通過節(jié)點(diǎn)處所代表的平衡方程和邊界條件求解整個線性方程組的過程。網(wǎng)格密度代表了求解單元的數(shù)量，其數(shù)量增加意味著求解方程數(shù)量的增加從而提高了計算精度，同時需要的計算時間也增加。因此，網(wǎng)格密度的選擇決定了計算過程的計算精度和效率，在這方面，已得到不同領(lǐng)域?qū)W者的關(guān)注。鄧記松[1]在用ANSYS軟件對筒體接管進(jìn)行應(yīng)力有限元分析時，研究了殼體厚度方向上網(wǎng)格密度(層數(shù)NT=2、3、4、5、6)對計算精度的影響，發(fā)現(xiàn)網(wǎng)格密度越大，總應(yīng)力值的精度越高，同時也關(guān)注到了在厚度方向劃分方法與計算精度的關(guān)系。符雙學(xué)[2]等基于Matlab和APDL軟件發(fā)現(xiàn)了網(wǎng)格密度(在模型齒寬方向分別劃分20、40、80、100、160份)增大，對齒根最大應(yīng)力和齒面平均變形的精度影響較小，計算出的齒根最大應(yīng)力差值僅為0.715 MPa，齒面平均變形最大差值僅為0.088 μm。蔣梅榮[3]等發(fā)現(xiàn)網(wǎng)格單元數(shù)在101 800～403 000范圍內(nèi)，對所計算的溢油傳播形態(tài)的精度沒有明顯提高，但電腦計算時間相差了3～6.5倍。在利用有限元法研究制動盤溫度場方面，也面臨著網(wǎng)格密度選擇的問題，楊鶯，王剛[4]借助ANSYS有限元軟件，在計算制動盤溫度時，劃分了8 458個網(wǎng)格單元，而吳婧斯[5]借助ABAQUS軟件對制動盤劃分了4 380個單元，吳波文[6]借助ADINA軟件對制動盤劃分了5 904個網(wǎng)格單元，可見，不同研究者在計算制動盤溫度時，考慮了不同的網(wǎng)格數(shù)量和計算軟件，但這些不同的計算方法中，采用的網(wǎng)格數(shù)量是否合理，計算精度如何保障，都是值得考慮的。因此，針對制動盤溫度場的計算，探討網(wǎng)格單元數(shù)對計算精度和效率的影響是有意義的。

文中針對盤面溫度場，利用ADINA有限元分析軟件，通過設(shè)置多種網(wǎng)格數(shù)量，研究了網(wǎng)格密度對計算精度和效率的影響。

1 制動摩擦副有限元模型與計算假設(shè)條件

考慮到摩擦副結(jié)構(gòu)的對稱性，取摩擦副結(jié)構(gòu)的一半建立有限元模型，由于每次制動產(chǎn)生的熱量還來不及擴(kuò)散到軸上，制動就已經(jīng)完成，所以數(shù)值模擬過程中將忽略軸的影響，簡化后的三角形摩擦副如圖1所示，其中盤做逆時針轉(zhuǎn)動。

假設(shè)條件如下：

(1)摩擦系數(shù)設(shè)為定值，不隨制動時間變化；

(2)忽略摩擦過程中磨損的影響，假設(shè)制動期間摩擦產(chǎn)生的熱量全部被制動盤和閘片所吸收；

(3)制動壓力均勻地施加在三角形閘片上；

(4)制動期間，制動盤和閘片的材料參數(shù)不隨時間變化；

(5)忽略制動過程中的熱對流及熱輻射的影響。

圖1 三角形摩擦副

2 網(wǎng)格劃分與邊界條件

采用ADINA軟件進(jìn)行有限元模擬，而有限元的核心就是利用離散化的網(wǎng)格進(jìn)行模擬計算。在ADINA的前處理模塊中，一共有3種指定網(wǎng)格密度的方法，即Use End-Point Sizes、Use Length和Use Number of Divisions，采用第3種指定網(wǎng)格細(xì)分份數(shù)的方法，閘片采用自由網(wǎng)格劃分，單元數(shù)為2 007個單元；由于制動盤是對稱結(jié)構(gòu)，故對其進(jìn)行了切分，采用映射網(wǎng)格進(jìn)行計算，通過改變制動盤徑向和周向上的細(xì)分份數(shù)從而探索其網(wǎng)格密度對模擬結(jié)果的影響，具體方案如表1所示，此外，閘片與制動盤都采用六面體8節(jié)點(diǎn)單元。

表1 網(wǎng)格單元數(shù)量方案

三角形摩擦副主要尺寸如表2所示，其中平均摩擦半徑是指三角形閘片的外接圓圓心與制動盤圓心的距離，摩擦面積為三角形閘片與制動盤的接觸面積。制動盤與閘片的物理參數(shù)如表3所示。

表2 制動盤與閘片的幾何尺寸

表3 制動盤與閘片的物理參數(shù)

模擬條件：制動初速度分別為50 km/h、100 km/h、160 km/h，制動壓力為0.5 MPa，制動盤和閘片的初始溫度都是20 ℃，熱流分配系數(shù)方面[7]，制動盤為0.54，閘片為0.46，對流換熱系數(shù)為100，摩擦系數(shù)為0.4。

3 結(jié)果分析與討論

3.1 半徑方向網(wǎng)格密度對制動盤溫度場的影響

表4是在制動速度50 km/h和制動壓力0.5 MPa條件下，針對4種網(wǎng)格密度(A1、A2、A3、A4)模擬的不同時刻盤面最高溫度情況，可見，制動時間為1 s時，4種方案模擬的盤面最高溫度皆為29 ℃，到3 s時，盤面最高溫度皆為34 ℃，到5 s時，盤面最高溫度皆為38 ℃，隨著制動時間的延長，盤面溫度升高，到10 s時，盤面最高溫度都達(dá)到43 ℃，到15 s時，A1、A3、A4盤面最高溫度為44 ℃，A2模擬的盤面最高溫度為43 ℃，19 s時，4種方案模擬的盤面最高溫度皆為41 ℃。

表5是在制動速度50 km/h和制動壓力0.5 MPa條件下，針對4種網(wǎng)格密度(A1、A2、A3、A4)模擬的制動期間盤面最高溫度和出現(xiàn)時間情況，可見，網(wǎng)格密度為A1、A2、A3時，盤面最高溫度為45 ℃，A4模擬的盤面最高溫度為46 ℃；最高溫度出現(xiàn)時間方面，A1和A2為14 s，A3和A4為13 s。

表4和表5的結(jié)果表明，在制動速度50 km/h，制動壓力0.5 MPa條件下，制動盤單元數(shù)在3 000～4 500范圍內(nèi)，從盤半徑方向改變網(wǎng)格密度對計算結(jié)果影響均不明顯，4種方案模擬的不同時刻盤面最高溫度差值是1 ℃，4種方案模擬的制動期間盤面最高溫度差值是1 ℃，4種方案模擬的制動期間盤面最高溫度所出現(xiàn)的時間差值是1 s。

表4 不同時刻模擬的盤面最高溫度 ℃

表5 不同方案模擬的制動期間盤面最高溫度和出現(xiàn)時間

圖2是在制動速度50 km/h和制動壓力0.5 MPa條件下，4種網(wǎng)格密度方案(A1、A2、A3、A4)模擬的不同時刻盤面溫度分布曲線，可見，4種方案對不同時刻的盤面溫度分布曲線影響不大，呈現(xiàn)出相同的變化規(guī)律：盤面高溫區(qū)聚集在91～160 mm范圍內(nèi)，這是因為此摩擦半徑范圍內(nèi)是制動盤與閘片的接觸區(qū)域，制動期間盤的動能轉(zhuǎn)化為熱能被制動盤和閘片所吸收，故盤溫上升較快；盤面最高溫度出現(xiàn)在摩擦半徑125 mm處左右，這是由閘片的幾何形狀決定的，此處摩擦弧長最長，故呈現(xiàn)出此處溫度最高，兩邊溫度逐漸放緩的分布趨勢；在摩擦半徑45～91 mm范圍內(nèi)，盤面熱量的主要來源是熱傳導(dǎo)，制動后期熱傳導(dǎo)作用增強(qiáng)，故制動前期此范圍內(nèi)盤面溫度幾乎等于室溫，制動后期，盤面溫度有所上升。

圖2 不同時刻盤面溫度分布曲線(50 km/h,0.5 MPa)

3.2 圓周方向網(wǎng)格密度對制動盤溫度場的影響

表6是在制動速度50 km/h和制動壓力0.5 MPa條件下，針對4種網(wǎng)格密度(B1、B2、B3、B4)模擬的不同時刻盤面最高溫度結(jié)果，可見，制動時間為1 s時，4種方案模擬的盤面最高溫度皆為29 ℃，到3 s時，盤面最高溫度皆為34 ℃，到5 s時，盤面最高溫度皆為38 ℃，隨著制動時間的延長，盤面溫度不斷升高，到10 s時，盤面最高溫度都達(dá)到43 ℃，到15 s時，盤面最高溫度皆為44 ℃，制動結(jié)束時，盤面最高溫度皆為41 ℃。

表6 不同時刻模擬的盤面最高溫度 ℃

表7是在制動速度50 km/h和制動壓力0.5 MPa條件下，針對4種網(wǎng)格密度(B1、B2、B3、B4)模擬的制動期間盤面最高溫度和出現(xiàn)時間結(jié)果，可見，B1、B2、B3、B4模擬的盤面最高溫度皆為45 ℃；在模擬最高溫度所出現(xiàn)的時間方面，B1、B3、B4為14 s，B2為13 s。

表7 不同方案模擬的盤面最高溫度、出現(xiàn)時間

表6和表7的結(jié)果表明，制動速度50 km/h，制動壓力0.5 MPa，制動盤單元數(shù)在3 000～4 500范圍內(nèi)，從盤圓周方向改變其網(wǎng)格密度對不同時刻盤面最高溫度，制動期間盤面最高溫度及制動期間盤面最高溫度所出現(xiàn)的時間方面影響均不明顯，4種方案模擬的不同時刻盤面最高溫度相同，4種方案模擬的制動期間盤面最高溫度相同，4種方案模擬的制動期間盤面最高溫度所出現(xiàn)的時間差值是1 s。

圖3是在制動速度50 km/h和制動壓力0.5 MPa條件下，4種網(wǎng)格密度方案(B1、B2、B3、B4)模擬的不同時刻盤面溫度分布曲線，可見，4種方案對不同時刻的盤面溫度分布曲線影響不大，4條曲線幾乎重合，不同時刻的曲線變化規(guī)律與圖2類似。

圖3 不同時刻盤面溫度分布曲線(50 km/h,0.5 MPa)

3.3 網(wǎng)格密度對計算時間的影響

圖4是在制動速度50 km/h和制動壓力0.5 MPa條件下，半徑方向4種不同網(wǎng)格密度(A1，A2,A3,A4)和圓周方向4種不同網(wǎng)格密度(B1，B2,B3,B4)的計算時間對比情況，制動盤單元數(shù)相同時，二者差值Δ=(Ai-Bi)/Bi×100%，式中Ai和Bi分別是半徑方向和圓周方向的相關(guān)物理量(i=1、2、3、4)。可見，A1，A2,A3,A4的計算時間關(guān)系為A1

圖4 網(wǎng)格密度對計算時間的影響

3.4 最佳網(wǎng)格密度研究

表8是在制動速度50 km/h和制動壓力0.5 MPa條件下，不同網(wǎng)格密度方案模擬的不同時刻盤面最高溫度，可見，B5至B8方案在不同時刻模擬的盤面最高溫度是相同的，B9和B10方案模擬的盤面最高溫度出現(xiàn)差別。

表8 不同時刻模擬的盤面最高溫度 ℃

表9是在制動速度50 km/h和制動壓力0.5 MPa條件下，不同網(wǎng)格密度方案模擬的制動期間盤面最高溫度、出現(xiàn)時間及計算時間情況，可見，B5至B8的盤面最高溫度相同，數(shù)值是45 ℃；B9和B10的盤面最高溫度分別是48 ℃和55 ℃，出現(xiàn)了一定的偏差；出現(xiàn)時間幾乎相同，計算時間方面，制動盤單元數(shù)越少，計算時間越短。

圖5是在制動速度50 km/h和制動壓力0.5 MPa條件下，6種不同網(wǎng)格密度方案模擬的不同時刻盤面溫度分布曲線，可見，B5、B6、B7模擬的盤面溫度分布曲線幾乎重合，B8、B9、B10模擬的溫度曲線與其比較可以發(fā)現(xiàn)，在摩擦半徑91～160 mm范圍內(nèi)發(fā)生了不同程度的偏移，不符合盤面溫度分布規(guī)律，所以B8、B9、B10方案中制動盤的網(wǎng)格密度是不合理的。制動盤網(wǎng)格密度不同造成盤面溫度分布曲線有差異的原因分析如下：有限元是先把要分析的模型做離散化處理分成有限個單元(網(wǎng)格劃分)，單元與單元之間用節(jié)點(diǎn)相連，一般情況下網(wǎng)格密度越大，劃分網(wǎng)格后的模型與實體模型越接近，相應(yīng)的計算結(jié)果越精確，如圖6是制動盤網(wǎng)格密度不同時的圖案分布，制動盤單元數(shù)為200時，制動盤網(wǎng)格圖案是一個四邊形，制動盤單元數(shù)為300時，制動盤網(wǎng)格圖案是一個六邊形，制動盤單元數(shù)為500時，制動盤網(wǎng)格圖案是一個十邊形，制動盤單元數(shù)為1 000時，制動盤網(wǎng)格圖案較接近于圓形。此外，4種網(wǎng)格圖案在實體制動盤上的覆蓋面積分別是63.7%,82.7%,93.5%和98.4%，由此可知，網(wǎng)格單元數(shù)越少，劃分網(wǎng)格后的圖形在實體制動盤上的覆蓋面積越小，導(dǎo)致實體模型中不參與計算的部分越多，從而導(dǎo)致計算出的結(jié)果與理想結(jié)果相差甚遠(yuǎn)。

表9 不同方案模擬的盤面最高溫度、出現(xiàn)時間及計算時間

圖5 不同時刻盤面溫度分布曲線(50 km/h,0.5 MPa)

圖6 制動盤單元數(shù)不同時圖案分布

圖7 不同時刻盤面溫度分布曲線(100 km/h,0.5 MPa)

圖8 不同時刻盤面溫度分布曲線(160 km/h,0.5 MPa)

綜合表9中計算時間和圖5的結(jié)果，可見，在制動速度50 km/h和制動壓力0.5 MPa條件下，B7方案在滿足精度要求的前提下且用時最少，為了驗證B7是最佳方案，對速度為100 km/h和160 km/h時進(jìn)行了模擬。

由圖7和圖8可知，制動盤單元數(shù)對盤面溫度場的影響規(guī)律與圖5類似，證明了制動盤單元數(shù)為1 000(B7方案)時，能達(dá)到精度與計算時間的完美結(jié)合。

4 結(jié) 論

(1)制動盤單元數(shù)決定了計算效率，當(dāng)單元數(shù)由300增加到4 500，計算時間從6.4 h延長到20.15 h。此外，劃分方法與計算時間也有關(guān)，相比于從盤半徑方向劃分，從盤圓周方向劃分更能提高計算效率，二者計算時間最大差值為7.59%。

(2)網(wǎng)格密度決定了計算精度，網(wǎng)格密度越大，劃分網(wǎng)格后的計算模型與實體模型越接近，計算精度越高。當(dāng)網(wǎng)格在實體模型上的覆蓋率達(dá)到98.4%(制動盤單元數(shù)為1 000)時，計算精度已達(dá)到100%。