董永兵，張義民，李鐵軍

（1.沈陽化工大學(xué)裝備可靠性研究所，遼寧 沈陽 110142；2.沈陽化工大學(xué)能源與動力工程學(xué)院，遼寧 沈陽 110142）

1 引言

軸承在現(xiàn)代機(jī)械傳動部件中有十分重要的地位，隨著科學(xué)技術(shù)不斷向前發(fā)展，對現(xiàn)代機(jī)械的使用壽命和運(yùn)行精度的要求越來越高，因此也使得人們對軸承的高精度和長壽命的需求越來越大[1]。但軸承在工作中，尤其是在惡劣的工況下，由于其內(nèi)各部件摩擦產(chǎn)生的巨大熱量使得軸承工作溫度異常升高，各零件受熱變形后直接影響傳動部件的工作精度和軸系的使用壽命，嚴(yán)重的會產(chǎn)生軸承溝道和滾動體燒傷甚至抱死現(xiàn)象。因此準(zhǔn)確及時(shí)的獲得軸承內(nèi)部溫度分布，可以對整個(gè)傳動系統(tǒng)進(jìn)行合理設(shè)計(jì)，對軸系部件實(shí)時(shí)調(diào)整使用工況，使得因溫度升高而易失效處得到合理散熱與潤滑，提高整個(gè)機(jī)械的使用精度、可靠性和壽命。所以建立滾承的傳熱模型，準(zhǔn)確快速的分析出軸承溫度場分布就顯得十分重要。

目前分析方法軸承溫度場的方法主要有熱阻網(wǎng)絡(luò)法、有限元法以及實(shí)驗(yàn)法，文獻(xiàn)[2]應(yīng)用熱阻網(wǎng)絡(luò)法建立傳熱模型計(jì)算軸承溫度場，得出部分節(jié)點(diǎn)的溫度變化規(guī)律。但該種方法不便處理復(fù)雜形狀，并對內(nèi)部溫度分布不能很好反應(yīng)；文獻(xiàn)[3]以傳熱學(xué)為理論基礎(chǔ)，應(yīng)用有限元法（FEM）分析出軸承穩(wěn)態(tài)溫度場，由于需要對網(wǎng)格進(jìn)行合理劃分，對問題的求解難度和計(jì)算量較大，求解速度較慢?？紤]采用無網(wǎng)格法對軸承溫度場進(jìn)行計(jì)算。

近幾年來，無網(wǎng)格方法快速向前發(fā)展，從開始應(yīng)用于計(jì)算力學(xué)逐步拓展到包括數(shù)值傳熱學(xué)等諸多領(lǐng)域。對比于FEM來說，無網(wǎng)格法只依靠離散節(jié)點(diǎn)的信息，避免了復(fù)雜的網(wǎng)格生成過程，同時(shí)也避免了低質(zhì)量網(wǎng)格對計(jì)算結(jié)果的影響，可以較為簡單的描述不規(guī)則的幾何形狀，保證計(jì)算精度的基礎(chǔ)上，節(jié)省了大量時(shí)間。文獻(xiàn)[4]用無網(wǎng)格Galerkin 方法（EFG）計(jì)算熱傳導(dǎo)問題，計(jì)算結(jié)果表明其具有較快的收斂速度，但由于采用背景網(wǎng)格積分，所以EFG不是真正的無網(wǎng)格方法。文獻(xiàn)[5]采用移動最小二乘法（MLS）構(gòu)建試探函數(shù)，將無網(wǎng)格Petrov-Galerkin法（MLPG）應(yīng)用在非穩(wěn)態(tài)導(dǎo)熱中，得到了較好的計(jì)算結(jié)果。文獻(xiàn)建立了加權(quán)最小二乘無網(wǎng)格法（MWLS）并成功應(yīng)用到穩(wěn)態(tài)與非穩(wěn)態(tài)數(shù)值傳熱中，計(jì)算結(jié)果表明該方法是一種精度高、速度快的數(shù)值傳熱分析方法。目前尚未發(fā)現(xiàn)將無網(wǎng)格傳熱方法應(yīng)用于軸承非穩(wěn)態(tài)溫度場計(jì)算。

根據(jù)滾動軸承溫度場分析和無網(wǎng)格法發(fā)展現(xiàn)狀，將MWLS應(yīng)用于深溝球軸承溫度場構(gòu)建，推導(dǎo)了相應(yīng)的計(jì)算公式，并編輯了相應(yīng)的程序，解決了有限元方法計(jì)算量大和計(jì)算時(shí)間較長的問題，得到了較為理想的結(jié)果。

2 無網(wǎng)格方法的實(shí)現(xiàn)

2.1 移動最小二乘近似

設(shè)待求函數(shù)u（x）在全域Ω內(nèi)分布有N個(gè)點(diǎn)，I=1，2，…N，且給出各節(jié)點(diǎn)處函數(shù)值，即uI=u（xI）。函數(shù)u（x）的近似表達(dá)式為uh（x），函數(shù)uh（x）在計(jì)算點(diǎn)x處的求解鄰域Ωx內(nèi)可以用式（1）表示［6-7］：

式中：pi（x）—基底函數(shù)，基底函數(shù)的構(gòu)造常采用多項(xiàng)式形式；m—基底函數(shù)的數(shù)量；ai（x）—基底函數(shù)對應(yīng)待定系數(shù)—計(jì)算點(diǎn)x的鄰域Ωx內(nèi)節(jié)點(diǎn)坐標(biāo)，p(x)=

在MLS近似中，ai（x）是在點(diǎn)x處的鄰域Ωx里通過使局部近似表達(dá)式uh(x，)與函數(shù)u（x）之間差的平方最小化得到，即對各個(gè)節(jié)點(diǎn)差值的加權(quán)平方之和取最小值，假設(shè)計(jì)算點(diǎn)鄰域有n個(gè)節(jié)點(diǎn)，由此可以列出如下表達(dá)式：

式中：ωI—非負(fù)權(quán)函數(shù)，與xI點(diǎn)相關(guān)，在點(diǎn)xI處權(quán)函數(shù)值取最大值，在支撐域邊界及域外權(quán)函數(shù)取值最小，即ω（r）=0，其中r=d/dmI，d—點(diǎn)x到節(jié)點(diǎn)xI之間的距離；dmI—權(quán)函數(shù)支撐域半徑，有dmI=s·c；s—支撐域半徑系數(shù)；c—節(jié)點(diǎn)之間距離。選取三次樣條權(quán)函數(shù)，即：

對J（x）其取極小值：

整理可得：

由式（5）可以解出待定系數(shù)向量a（x）：

其中，

將a（x）帶回式（1）可得到：

如果對全域Ω中所有點(diǎn)x都在鄰域Ωx內(nèi)構(gòu)造上述局部最佳近似表達(dá)式，將其組合后就得到了待求函數(shù)u（x）在求解域Ω內(nèi)的全局最佳近似函數(shù)uh（x），即：

其中，形函數(shù)ΦI（x）為：

其中，EI（x）—矩陣E（x）的第I列。

2.2 非穩(wěn)態(tài)傳熱中MWLS模型構(gòu)建

計(jì)算對象的溫度場隨時(shí)間而發(fā)生變化的導(dǎo)熱過程為非穩(wěn)態(tài)導(dǎo)熱，如假定材料導(dǎo)熱系數(shù)不隨溫度變動，即為常數(shù)時(shí)，則其熱微分方程為［8-9］：

邊界約束條件可以表示為：

對于非穩(wěn)態(tài)情況，指定初始條件：

式中：k—材料的導(dǎo)熱系數(shù)；ρQ—熱源項(xiàng)；n—各邊界面外法線方向—第一類邊界Γ1上的規(guī)定的溫度—第二類邊界Γ2上沿n方向的熱流密度；h—表面散熱系數(shù)；Tw—周圍溫度；ρ為材料密度；c—材料比熱。

對于非穩(wěn)態(tài)導(dǎo)熱問題涉及對時(shí)間求導(dǎo)數(shù)，因此采用向后差分改寫控制式（10）可得：

假設(shè)熱源項(xiàng)不隨時(shí)間改變，經(jīng)過整理可將式（13）改寫為：

式中：n—時(shí)間步長，式（14）左面第（n+1）個(gè)時(shí)間步長的結(jié)果可以完全由第n個(gè)時(shí)間步長的結(jié)果唯一確定，可以按照MWLS方法求解。

MWLS是以加權(quán)殘量的方式來求解偏微分方程，定義非穩(wěn)態(tài)傳熱的全域和各邊界的殘量為R和Ri（i=1，2，3），其殘量表達(dá)式為：

將式（15）中各殘量分別平方并加權(quán)積分后相加，可得到如下泛函：

由于積分的計(jì)算量較大，采取離散的方法將其泛函改為：

式中：λs（s=1，2，3）表示三種約束邊界條件所對應(yīng)的罰函數(shù)；n—控制式（10）的殘量的計(jì)算節(jié)點(diǎn)個(gè)數(shù)；nl（l=1，2，3）為各自對應(yīng)約束條件殘量的計(jì)算節(jié)點(diǎn)個(gè)數(shù)，對式（17）的泛函取極小值，即：

可以得到求解非穩(wěn)態(tài)傳熱方程的加權(quán)最小二乘無網(wǎng)格方法（MWLS）的系統(tǒng)方程，并用簡單的迭代可以解出各個(gè)時(shí)刻的溫度分布：

其中，

上述計(jì)算過程中，將三種邊界條件以罰參數(shù)的方式引入求解，為了平衡不同約束條件殘量的數(shù)量級，并且使邊界條件的殘差相對于控制方程來說在求解中占主導(dǎo)地位，經(jīng)過量綱分析定義λ1=λ2（k/L）2，λ3=λ2[min（1，k/hL）]2，其中參數(shù)λ2推薦值為（105～108）。

3 實(shí)例分析及試驗(yàn)驗(yàn)證

3.1 軸承溫度場分析

以CNC機(jī)床滾珠絲杠進(jìn)給驅(qū)動中預(yù)拉伸端軸承25TAC62B（NSK）為研究對象。采用MWLS建立非穩(wěn)態(tài)溫度場模型。設(shè)置外部環(huán)境溫度23℃，軸承工作轉(zhuǎn)速1500r/min。滾動軸承的發(fā)熱大部分來源于滾珠和內(nèi)外滾道的相互滑動和摩擦，根據(jù)文獻(xiàn)[10]計(jì)算出軸承發(fā)熱量為Htot=2.79W。由Steph 和Burton 所提出的觀點(diǎn)，其軸承運(yùn)行時(shí)產(chǎn)生的熱量50%進(jìn)入滾動體，另50%傳入外圈和內(nèi)圈。25TAC62B 軸承結(jié)構(gòu)參數(shù)為：公稱外徑De=62mm；公稱內(nèi)徑Di=25mm；外圈內(nèi)徑de=51.57mm；內(nèi)圈外徑di=37.37mm；滾動體直徑D=6.9mm；內(nèi)溝道溝半徑ri=3.65mm；外滾到溝半徑re=3.65mm；滾動體數(shù)目Z=18；接觸角α=60°。

由于軸承外圈外表面與軸承座內(nèi)表面緊密連接，忽略其接觸熱阻，所以將研究對象簡化為三個(gè)部分，分別為軸承內(nèi)圈、球和軸承外圈加軸承座，如圖1所示。

圖1 軸承節(jié)點(diǎn)分布簡圖Fig.1 Schematic Diagram of Bearing Node Distribution

對軸承內(nèi)圈進(jìn)行分析，其受到軸承摩擦熱和軸端散熱影響，內(nèi)圈溝道為該部分摩擦熱的影響面，經(jīng)過計(jì)算內(nèi)圈溝道面積s1=711.93mm2。沿圖1中x軸方向均勻分布8個(gè)計(jì)算節(jié)點(diǎn)，得到內(nèi)圈溝道溫度隨時(shí)間變化曲線，如圖2所示。得到各節(jié)點(diǎn)在不同時(shí)刻溫度場分布，如圖3所示。

圖2 軸承內(nèi)圈溝道溫度隨時(shí)間變化圖Fig.2 Variation of Temperature of Bearing Inner Ring Groove with Time

圖3 軸承內(nèi)圈各計(jì)算點(diǎn)在不同時(shí)刻溫度變化圖Fig.3 Temperature Variation Diagram of Each Calculation Point of Bearing Inner Ring at Different Time

計(jì)算結(jié)果表明軸承內(nèi)圈溫度在運(yùn)行初期溫升較大，在35min附近接近穩(wěn)態(tài)，最終達(dá)到穩(wěn)態(tài)溫度為45.03℃。由于軸承內(nèi)圈較小，軸端散熱率低，各個(gè)計(jì)算點(diǎn)溫度趨于一致比較容易，所以在同一時(shí)間各節(jié)點(diǎn)間溫差較小。對滾珠進(jìn)行分析，滾珠受到兩端軸承摩擦熱和潤滑油持續(xù)散熱影響，將滾珠影響面簡化為滾珠的一半表面，經(jīng)計(jì)算s2=s3=172.93mm2。將潤滑油造成的球表面散熱轉(zhuǎn)化為隨各點(diǎn)溫度變化的負(fù)的內(nèi)熱源。經(jīng)計(jì)算得滾動體中心點(diǎn)溫度隨時(shí)間增長曲線，如圖4所示。其在35min左右達(dá)到38.04℃后趨于穩(wěn)定。沿圖1中x軸所示方向均勻分布14個(gè)計(jì)算節(jié)點(diǎn)，得到各節(jié)點(diǎn)在不同時(shí)間下的溫度分布規(guī)律，如圖5所示。

圖4 滾珠中心點(diǎn)溫度隨時(shí)間變化示意圖Fig.4 Temperature Changes of Ball Center Point with Time

圖5 滾珠各計(jì)算點(diǎn)在不同時(shí)刻溫度變化圖Fig.5 Temperature Variation Diagram of Each Calculation Point of the Ball at Different Moments

對軸承外圈及軸承座部分進(jìn)行計(jì)算，該整體主要受到軸承摩擦生熱和軸承座對流散熱的影響，外圈滾道為該部分摩擦熱的影響面，經(jīng)過計(jì)算外圈滾道面積s4=1187.06mm2，對整個(gè)研究對象均勻劃分18個(gè)計(jì)算節(jié)點(diǎn)，得到軸承外圈滾道和軸承座表面溫度隨時(shí)間變化規(guī)律，如圖6所示。以及不同時(shí)間軸承溫度場分布狀況，如圖7所示。結(jié)果表明外圈溝道和軸承座在35min后溫度分別在29.51℃和27.5℃達(dá)到穩(wěn)態(tài)。

圖6 軸承外圈溝道與軸承座溫度隨時(shí)間變化圖Fig.6 Variation of the Temperature of Outer Ring Channel and Bearing Seat with Time

圖7 軸承外圈與軸承座各計(jì)算點(diǎn)在不同時(shí)刻溫度變化圖Fig.7 Temperature Changes of the Calculated Points of the Bearing Outer Ring and Bearing Seat at Different Times

3.2 試驗(yàn)驗(yàn)證

該試驗(yàn)由滾珠絲杠驅(qū)動系統(tǒng)和精度為0.1℃的FLUKE熱成像設(shè)備組成，實(shí)驗(yàn)裝置，如圖8所示。將拉伸端軸承作為測試對象，測點(diǎn)位置為軸承座外頂面中點(diǎn)。該軸承材料的物理特性參數(shù)為，密度ρ=7800kg/m3；比熱容c=460J/kg·K；材料導(dǎo)熱系數(shù)k=39W/（m·K）。其中機(jī)床滾珠絲杠驅(qū)動系統(tǒng)的絲杠工作轉(zhuǎn)速為1500r/min。

圖8 實(shí)驗(yàn)裝置Fig.8 Experimental Device

測試時(shí)，首先測量試驗(yàn)的環(huán)境溫度等于23.0℃，電機(jī)開機(jī)并以轉(zhuǎn)速為1500r/min驅(qū)動絲杠暖機(jī)，試驗(yàn)時(shí)間累計(jì)為1h，采用熱成像儀對溫度測點(diǎn)時(shí)間間隔10min采樣一次，采樣熱成像圖以及測點(diǎn)位置示意，如圖9所示。最后使用SmartView讀取并保存的測試結(jié)果。

圖9 熱成像圖及測點(diǎn)位置Fig.9 Heat Image and Position of Measuring Point

將軸承座計(jì)算結(jié)果與試驗(yàn)數(shù)據(jù)匯總，如圖10所示。經(jīng)對比驗(yàn)證可得該溫度辨識結(jié)果與試驗(yàn)結(jié)果均呈指數(shù)上升趨勢，趨勢合理。軸承座外圈計(jì)算溫度與試驗(yàn)溫度最大差值為0.637℃，其最大誤差為2.46%。

圖10 軸承座溫度實(shí)驗(yàn)值和計(jì)算值對比Fig.10 Comparison of Experimental and Calculated Values of the Bearing Seat Temperature

4 結(jié)論

提出了基于無網(wǎng)格傳熱的軸承溫度場建模方法，得到了一般工況下軸承各元件的瞬態(tài)溫度場。該計(jì)算方法獲得了軸承各部分沿徑向溫度分布規(guī)律，由于不需進(jìn)行網(wǎng)格的劃分，并且該方法所得到的溫度剛度矩陣為帶狀矩陣，計(jì)算規(guī)模減少，計(jì)算速度提高明顯，并且與實(shí)驗(yàn)結(jié)果進(jìn)行對比，證明該模型是一種快速準(zhǔn)確的計(jì)算軸承瞬態(tài)溫度場的數(shù)值方法。