基于MATLAB的微軸承動力學(xué)特性的數(shù)值分析*

蘭正義，伍 林*，呂剛磊

(1.上海應(yīng)用技術(shù)大學(xué) 機(jī)械工程學(xué)院，上海 201418;2.鄭州科技學(xué)院 機(jī)械工程學(xué)院，河南 鄭州 450064)

0 引 言

目前，微機(jī)電系統(tǒng)下微軸承—轉(zhuǎn)子系統(tǒng)穩(wěn)定性的判斷分析，主要是結(jié)合微軸承的動力學(xué)特性系數(shù)進(jìn)行分析而獲得，但大部分的動力學(xué)特性參數(shù)計算是在基于不考慮氣體滑移的連續(xù)模型下獲得的[1-2]。

然而隨著微機(jī)電系統(tǒng)尺寸的不斷微型化，其關(guān)鍵支撐部件軸承間的潤滑氣膜將會變得極薄，氣體將會處于稀薄狀態(tài)，其間的氣體將會在固壁表面發(fā)生滑移，此時若采用不考慮滑移影響的傳統(tǒng)連續(xù)模型分析軸承轉(zhuǎn)子系統(tǒng)的動力學(xué)特性，將影響其準(zhǔn)確性。因此，需考慮滑移效應(yīng)下提出的修正模型對其進(jìn)行動力學(xué)特性分析[3-4]。文獻(xiàn)[5]指出了當(dāng)軸承最小氣膜間隙處的克努森數(shù)Kn值處于氣流過渡流動區(qū)域(0.1

針對以上問題，以MEMS下的微軸承為研究對象，本文推導(dǎo)考慮一階滑移修正及Wu新滑移修正的動力學(xué)雷諾方程，采用MATLAB編程軟件編制不同模型下用于數(shù)值計算動力學(xué)雷諾方程的程序，比較不同軸承數(shù)、偏心率、軸頸擾動頻率下連續(xù)模型、一階滑移模型和Wu新滑移模型的動力學(xué)特性系數(shù)，進(jìn)而分析研究氣體滑移對微軸承動力學(xué)特性的影響規(guī)律；在此數(shù)值計算基礎(chǔ)上，對比分析微軸承剛性轉(zhuǎn)子系統(tǒng)在連續(xù)模型、一階滑移模型和Wu新滑移模型下的穩(wěn)定性。

1 滑移修正下的動力學(xué)方程

1.1 微軸承的滑移修正雷諾方程

微動壓氣體軸承，依靠軸承與軸套之間充入的黏性氣體，形成楔形間隙而發(fā)生相對運(yùn)動，從而自動產(chǎn)生氣膜壓力，實(shí)現(xiàn)軸承轉(zhuǎn)子的動壓懸浮。

微軸承動壓作用示意圖如圖1所示。

圖1 微軸承動壓作用示意圖ω—角速度；e—偏心距；h—軸承間隙

在定常、等溫條件下，傳統(tǒng)連續(xù)、一階滑移及Wu新滑移修正模型下的潤滑Reynolds方程的無量綱形式為：

(1)

為方便計算，筆者將3種模型的系數(shù)做歸納處理。

模型系數(shù)如表1所示。

表1 模型系數(shù)

Kn0—參考克努森數(shù),Kn0=C1/ph=C1/(pahmPH)；C1—最小氣膜間隙，C1=6.8×10-3；α—協(xié)調(diào)系數(shù)；f—取值函數(shù)，f=min[1/Kn0,1]

1.2 滑移修正下的靜、動態(tài)雷諾方程

結(jié)合式(1)、表1可得到在小擾動工況下修正的靜態(tài)雷諾方程，及關(guān)于PE和PΘ的修正偏微分動態(tài)雷諾方程。

修正的靜態(tài)雷諾方程為：

(2)

靜態(tài)氣膜厚度為：

H0=1+εcosφ

(3)

無量綱靜態(tài)求解邊界條件為：

(4)

關(guān)于PE修正后的動態(tài)雷諾方程為：

(5)

相應(yīng)的動態(tài)氣膜厚度為：

HE=cosφ

(6)

相應(yīng)的無量綱邊界條件為：

(7)

關(guān)于PΘ修正后的動態(tài)雷諾方程為：

(8)

相應(yīng)的動態(tài)氣膜厚度為：

HΘ=sinφ

(9)

相應(yīng)的無量綱邊界條件為：

PΘ(φ,±B/2)=0

(10)

1.3 動力學(xué)特性系數(shù)求解表達(dá)式

1.3.1 靜態(tài)承載力、偏位角求解表達(dá)式

聯(lián)立式(2～4),基于MATLAB編制靜態(tài)雷諾方程有限差分程序，迭代求解可得到靜態(tài)氣膜壓力P0分布的數(shù)值解，從而得到軸承靜態(tài)承載力W和軸頸偏位角θ的求解表達(dá)式，即：

(11)

(12)

式中：Wn—軸頸承載力在偏心方向的分量；Wt—軸頸承載力在垂直偏心方向的分量。

1.3.2 動態(tài)剛度系數(shù)和動態(tài)阻尼系數(shù)求解表達(dá)式

根據(jù)式(5～10),基于MATLAB編制動態(tài)雷諾方程有限差分程序，通過迭代求解可得到微軸承的動態(tài)氣膜壓力PE和PΘ，從而可以得到滑移效應(yīng)下，微軸承的動特性系數(shù)求解表達(dá)式：

(13)

轉(zhuǎn)換式(13)中的系數(shù)代入到直角坐標(biāo)系中，可求解直角坐標(biāo)系下的動特性系數(shù)，即：

(14)

2 微軸承轉(zhuǎn)子系統(tǒng)穩(wěn)定分析方法

結(jié)合文獻(xiàn)[8]可知，該系統(tǒng)的力平衡運(yùn)動方程為：

(15)

(16)

將上述微分方程組合并，則該方程組的通解為：

(x,y)=(x0,y0)est

(17)

式中：s—為復(fù)頻率，s=u+iwj；u—復(fù)數(shù)的實(shí)部(阻尼)；ωj—虛部(軸頸自然渦動頻率)。

聯(lián)立式(15～17)，整理可得：

(18)

由線性代數(shù)理論可知，式(18)有非零解的充要條件是該方程組的系數(shù)行列式為0，即：

(19)

展開式(19)的2階行列式，可得到微軸承—剛性轉(zhuǎn)子系統(tǒng)的特征方程式：

a4s4+a3s3+a2s2+a1s+a0=0

(20)

(21)

根據(jù)系統(tǒng)穩(wěn)定性Routh-Hurwitz判別準(zhǔn)則[9]可知，只需要判斷特征值s的解是否都有負(fù)實(shí)數(shù)部，就可以判斷軸承轉(zhuǎn)子系統(tǒng)是否穩(wěn)定。因此，通過計算式(20)中的特征值，就可預(yù)估微軸承剛性轉(zhuǎn)子系統(tǒng)的穩(wěn)定性。

3 數(shù)值仿真實(shí)驗(yàn)及結(jié)果分析

3.1 數(shù)值仿真程序驗(yàn)證

針對如圖1所示的微軸承結(jié)構(gòu)，微軸承設(shè)定最終微小間隙下氣膜收斂精度為10-6。

數(shù)值仿真計算參數(shù)如表2所示。

表2 數(shù)值仿真計算參數(shù)

為進(jìn)行微軸承動力學(xué)特性的數(shù)值分析[10]，本文采用MATALB編程平臺編制程序，進(jìn)行微軸承動力學(xué)特性的數(shù)值分析。

數(shù)值仿真流程圖如圖2所示。

圖2 數(shù)值仿真流程圖

為驗(yàn)證本文用于微軸承動力學(xué)特性分析時程序的正確性，基于如圖2所示的數(shù)值仿真流程圖，筆者采用有限差分法進(jìn)行數(shù)值計算，并將得到的結(jié)果與文獻(xiàn)[11]的數(shù)值結(jié)果進(jìn)行對比。

數(shù)值仿真程序驗(yàn)證圖如圖3所示。

圖3 數(shù)值仿真程序驗(yàn)證圖

對比結(jié)果表明：本文程序得到的動特性系數(shù)與文獻(xiàn)中的數(shù)據(jù)幾乎完全一樣；使用該程序計算得到的動力學(xué)特性系數(shù)與文獻(xiàn)中的數(shù)據(jù)相對誤差均在3%以內(nèi)，由此驗(yàn)證了本文程序的準(zhǔn)確性。

3.2 靜特性系數(shù)數(shù)值仿真結(jié)果分析

為分析微軸承的靜態(tài)特性，筆者從微軸承靜態(tài)雷諾方程出發(fā)，對比分析最小膜厚度處微軸承在3種模型下的靜態(tài)氣膜壓力分布及軸承數(shù)(描述氣體可壓縮性對動壓效應(yīng)影響程度的一個壓縮系數(shù)，標(biāo)量)、偏心率(描述軸頸運(yùn)轉(zhuǎn)時的偏心程度，標(biāo)量)分別對軸承的靜態(tài)承載力、偏位角的影響規(guī)律。

基于如表2所示的軸承結(jié)構(gòu)參數(shù)，在φ=180°、偏心率0.6及軸承數(shù)為1.7時，3種模型下的氣膜壓力分布圖如圖4所示。

圖4 3種模型下的氣膜壓力分布圖

對比3種模型下的氣膜壓力分布發(fā)現(xiàn)：在最小氣膜厚度處，3種模型得到的氣膜壓力差異明顯，這是因?yàn)樵谧钚饽ず穸忍帲蜌怏w滑移效應(yīng)密切相關(guān)的流體粘性剪切應(yīng)力最強(qiáng)，致使在該區(qū)域的氣流速度分布發(fā)生變化，進(jìn)而使該區(qū)域的氣壓分布發(fā)生變化；同時，可以得到連續(xù)模型預(yù)測的氣膜壓力最高，一階滑移模型次之，Wu新滑移模型預(yù)測的氣膜壓力最低。

為進(jìn)一步分析微軸承的靜態(tài)特性，筆者采用數(shù)值分析方法，得到3種模型下承載力與偏心率、軸承數(shù)之間的關(guān)系，如圖5所示。

圖5 3種模型下承載力與偏心率、軸承數(shù)之間的關(guān)系

從圖5可以看出：在相同條件下，考慮氣體滑移影響時，隨著偏心率、軸承數(shù)的變化，氣體滑移影響下的承載力明顯下降，相比滑移模型預(yù)測的結(jié)果而言，連續(xù)模型預(yù)測的承載力偏大；且隨著偏心率、軸承數(shù)的增大，各模型之間的承載力偏差在逐漸變大，這是因?yàn)槠穆屎洼S承數(shù)的增大，直接導(dǎo)致了最小氣膜處的間隙變大，剪切應(yīng)力增強(qiáng)，氣體滑移效應(yīng)變得顯著，進(jìn)而使最小氣膜間隙處的氣膜壓力降低，軸承的承載力下降。但連續(xù)模型并沒有考慮氣體滑移的影響，所以預(yù)測得到的承載力值相對較高。

綜合上述的分析可知，連續(xù)模型預(yù)測的承載力偏差最大、一階滑移模型次之、Wu新滑移模型最??；對比一階滑移模型和Wu新滑移模型預(yù)測的軸承氣膜壓力可知，一階滑移模型較高地預(yù)估了滑移效應(yīng)下軸承的承載力。

因此，在對滑移下微軸承的動力學(xué)特性分析時，需考慮精度更高的Wu新滑移模型。

3.3 動態(tài)特性系數(shù)數(shù)值仿真結(jié)果分析

為分析氣體滑移對微軸承動態(tài)特性的影響，筆者按表2所示的軸承結(jié)構(gòu)參數(shù)(氣膜收斂精度為10-6)，從微軸承動態(tài)雷諾方程出發(fā)，研究了氣體滑移影響下軸頸擾動頻率、偏心率、軸承數(shù)與軸承動態(tài)剛度和動態(tài)阻尼系數(shù)的變化關(guān)系。

3種模型下，動態(tài)特性系數(shù)與偏心率之間的關(guān)系如圖6所示。

圖6 3種模型下動態(tài)特性系數(shù)與偏心率之間的關(guān)系

從圖6中可以看出：

(1)在相同條件下，考慮氣體滑移的影響時，軸承的動態(tài)剛度系數(shù)，受偏心率的影響較大，軸承的動態(tài)阻尼系數(shù)，幾乎不受偏心率的影響；其中，當(dāng)軸頸的擾動頻率大于2.5時，動態(tài)主阻尼系數(shù)幾乎不受偏心率的影響；

(2)在相同軸頸擾動頻率下，大偏心率對應(yīng)著大的動態(tài)主剛度系數(shù)和主阻尼系數(shù)，且由于只有垂直方向受外力作用，垂直方向上大偏心率對應(yīng)的動態(tài)主剛度系數(shù)Kyy、主阻尼系數(shù)Dyy相對小偏心率較大。

基于靜特性數(shù)值分析計算結(jié)果，3種模型下動態(tài)特性系數(shù)與軸承數(shù)之間的關(guān)系如圖7所示。

圖7 3種模型下動態(tài)特性系數(shù)與軸承數(shù)之間的關(guān)系

從圖7中可以看出：在相同條件下，考慮氣體滑移影響時，不同軸承數(shù)下的動態(tài)主剛度系數(shù)Kxx、Kyy均隨軸頸擾動頻率的增大而增大，主阻尼系數(shù)Dxx、Dyy隨軸頸擾動頻率的增大而減小，且氣體滑移對阻尼系數(shù)的影響也在逐漸減弱；大的軸承數(shù)對應(yīng)著大的動態(tài)主剛度系數(shù)、小的動態(tài)主阻尼系數(shù)，這主要是因?yàn)?，大轉(zhuǎn)速意味著高氣膜壓力，會使相對應(yīng)的軸承軸頸位置發(fā)生變化，軸承的動態(tài)主剛度隨著轉(zhuǎn)速的增加而變大，相對應(yīng)的動態(tài)主阻尼系數(shù)隨之變小。

由圖(6，7)可知，在氣體滑移影響下，不同軸承數(shù)、偏心率對應(yīng)的動特性系數(shù)均有所下降，且Wu新滑移模型計算的值遠(yuǎn)低于連續(xù)模型和一階滑移模型，這主要是因?yàn)闅怏w滑移會使?jié)櫥瑲怏w的速度有所降低，進(jìn)而氣膜壓力降低，從而導(dǎo)致軸承動特性系數(shù)減小。

4 微軸承轉(zhuǎn)子系統(tǒng)穩(wěn)定性分析

在偏心率為0.6、軸承數(shù)為1.7時，結(jié)合3種模型下的動特性系數(shù)與軸承轉(zhuǎn)子系統(tǒng)穩(wěn)定性判別理論，可求得各擾動頻率下特征值對應(yīng)的所有子行列式的值，判斷在該擾動頻率下是否滿足穩(wěn)定性條件，為軸承穩(wěn)定性設(shè)計提供理論依據(jù)。

在相應(yīng)擾動頻率下，3種模型各子階行列式的值如表3所示。

表3 在相應(yīng)擾動頻率下3種模型各子階行列式的值

從表3可知，3種模型在各擾動頻率下的各子行列式的值都為正數(shù)，且Wu新滑移模型計算的值相比連續(xù)模型、一階滑移模型低；即在軸承數(shù)為1.7的情況下，微軸承剛性轉(zhuǎn)子在相應(yīng)擾動頻率范圍內(nèi)是穩(wěn)定的，且連續(xù)模型較高地估計了轉(zhuǎn)子系統(tǒng)的穩(wěn)定性。

5 結(jié)束語

針對微機(jī)電系統(tǒng)下傳統(tǒng)連續(xù)模型無法準(zhǔn)確描述氣體滑移對微軸承動力學(xué)特性影響規(guī)律的問題，筆者推導(dǎo)了考慮一階滑移修正及Wu新滑移修正的動力學(xué)雷諾方程，比較了不同軸承數(shù)、偏心率、軸頸擾動頻率下連續(xù)模型、一階滑移模型和Wu新滑移模型的動力學(xué)特性系數(shù)，進(jìn)而研究了氣體滑移對微軸承動力學(xué)特性的影響規(guī)律；在此數(shù)值計算基礎(chǔ)上，筆者對比分析了微軸承剛性轉(zhuǎn)子系統(tǒng)在連續(xù)模型、一階滑移模型和Wu新滑移模型下的穩(wěn)定性。

研究結(jié)果可知，氣體滑移下軸承的動力學(xué)特性與軸承的軸承數(shù)、偏心率及軸頸擾動頻率密切相關(guān)，在微小氣膜間隙、高轉(zhuǎn)速工況下，不考慮氣體滑移影響的連續(xù)模型較高地估計了軸承剛性轉(zhuǎn)子系統(tǒng)的穩(wěn)定性，且Wu滑移模型下得到的值最小。

綜上所述，結(jié)合物理理論依據(jù)和微小間隙下氣體滑移的影響，在設(shè)計微機(jī)電系統(tǒng)下的微軸承時，為得到氣體滑移下更準(zhǔn)確的軸承設(shè)計參數(shù)，應(yīng)使用Wu新滑移模型對傳統(tǒng)連續(xù)模型下的動力學(xué)雷諾方程進(jìn)行滑移修正。

為計算方便，筆者使用的MEMS下的微軸承結(jié)構(gòu)較為簡單。在今后的研究中，筆者將會引入更復(fù)雜的微軸承結(jié)構(gòu)，來驗(yàn)證分析該理論的適用性；同時，考慮使用數(shù)值仿真與試驗(yàn)相結(jié)合的方法，筆者在不同的算例條件下，分析氣體滑移對微軸承動力學(xué)特性的影響，為設(shè)計更優(yōu)的MEMS下微軸承提供幫助。