柔性軸-盤轉子系統的等頻線及其振動特性研究

唐進元， 馮俊易， 關先磊

(中南大學 機電工程學院，長沙 410083)

軸-盤轉子系統建模與分析是轉子動力學的基礎研究內容。例如渦輪轉子系統、計算機硬盤驅動器系統[1-4]，齒輪系統[5-6]等均需要精準的軸-盤轉子系統的建模與分析。

目前多采用有限元法和假設模態法等數值方法建立軸-盤轉子系統的運動微分方程。Shahab等[7]建立了最多3個均勻分布轉盤的軸-盤轉子數學模型，并分析了轉盤參數中“厚度/外半徑”比值對固有頻率的影響，當“厚度/外半徑”比值很小時，耦合的固有頻率接近于未耦合的固有頻率；相反，耦合效應將更加明顯。

Lee等[1]采用子結構綜合法和假設模態法研究了軸-葉盤系統上的葉盤柔性，指出了只有一個節徑的轉盤模式將與軸耦合。后來在另一項研究中，使用假設模態法研究了短軸雙盤系統和長軸雙盤系統。并發現，轉盤的同相模態振型下傾向于強化轉軸的凈轉動慣量，而異相模態振型下則傾向于削弱轉軸的凈轉動慣量。

此外，Lim和Jang等研究了硬盤驅動器系統，前者采用有限元方法和試驗結果對比，指出可以將該系統劃分成少量的有限單元網格，得到系統柔性耦合振動特性；后者采用了有限元法和子結構綜合法研究，考慮了包括滾珠軸承在內的柔性支承結構的硬盤驅動系統，說明了柔性支承結構的剛性連接約束對準確預測硬盤驅動系統的固有頻率具有重要作用。Jia[9]采用假設模態法建立了柔性軸-多轉盤系統，發現轉盤的柔性的增加通常會降低系統的固有頻率。

近年來，軸-盤轉子系統建模與分析的研究是熱點研究內容。但系統中一些參數對動力學特性的影響還有待于拓展和完善。Heydari等[10]分析了柔性轉盤的位置和“徑厚比(轉盤外圓半徑/轉盤厚度)”對軸-盤轉子系統的固有頻率和臨界速度的影響，但是并沒有解釋在結合頻率變化規律時考慮轉盤位置、轉盤外半徑對該轉子系統動力學特性的影響，有些轉盤位置和徑厚比對動力學特性影響的結論值得進一步商榷。因此，本文提出等頻線的概念以描述頻率變化規律(以相對轉盤位置、轉盤相對外半徑為橫坐標和縱坐標繪制固有頻率的變化圖，根據圖中所呈現的線條規律，固有頻率的值在水平范圍內或垂直范圍內保持不變，而形成的線條本文將其定義為“等頻線”。)，在此基礎上分析轉盤位置和轉盤外半徑對系統坎貝爾圖和模態的影響。

本文首先通過無量綱物理量(軸半徑ri、轉盤外半徑rD、轉盤厚度hD和轉盤位置lD)參數化分析，揭示這些參數對柔性軸-盤轉子系統的耦合振動特性中固有頻率的影響，采用“等頻線”描述參數化研究中固有頻率的變化規律，在垂直等頻線和水平等頻線基礎上選取特定的參數變化，判定相應參數對系統動力學特性的影響，本文研究結果可為研究柔性軸-盤轉子系統振動特性時提供結構參數的參考依據。

1 柔性軸-盤轉子系統的動力學模型

1.1 坐標系

具有簡支邊界條件的軸-盤轉子系統的示意圖,如圖1所示。為了描述圓盤的柔性變形和剛性運動，引入了3個坐標系：

圖1 軸-盤轉子系統原理圖及坐標系示意圖

(1)O0X0Y0Z0為固定慣性坐標系，用C0坐標系表示，Z0軸為沿轉軸的軸線。

(2)O1X1Y1Z1為以恒定速度Ω,r/s，繞O0X0Y0Z0的Z0軸建立的旋轉坐標系，用C1坐標系表示，用于描述軸的彎曲振動和轉盤的剛性運動，其原點O1與O0一致，Z1軸與Z0軸一致。

(3)O2X2Y2Z2為以O2為原點建立的擺動坐標系，用C2坐標系表示，可以用來描述圓盤的橫向振動，O2為轉盤部分的中心點同時也是轉軸部分的中心點，Z2軸垂直于轉盤和軸截面。

(a)

e1=Te2

T=

(1)

1.2 轉軸和轉盤上任意一點的速度

圓形截面的Timoshenko軸的柔性變形位移在C1坐標系中定義，分別用沿X1方向的uS和沿Y1方向vS表示， 因此，任意軸截面的中心點的位置矢量為C1rO1O2=[uSvSz]e1,C1rO1O2左上標C1為該矢量在C1坐標系中定義，下標O1O2為該矢量從O1點指向O2點，z為軸截面的坐標。在C2坐標系中，轉軸截面上任意一點相對于轉軸截面中心的位置矢量可以用極坐標的形式表示為C2rO2S=[rcosθrsinθ0]e2，其中位置矢量下標O2S為從O2點產生的速度矢量。然后可以得到轉軸上任意一點相對于O1坐標系的位置矢量為rS=C1rO1O2+C2rO2S，對時間t求導得軸上任意點的絕對速度

(2)

式中,ωC0C1和ωC0C2分別為C1和C2相對于慣性系C0的瞬時角速度，此外，若不考慮表示軸扭轉振動的扭轉角θz，ωC0C2的線性化形式可以表示為

(3)

式中，Ω為軸-盤轉子系統轉速。忽略3階項和3階以上高階項，轉軸中任意點的速度表示為

(4)

注意到，所有的位移都是坐標z的函數。基于經典薄板理論建立的薄環形圓盤，如圖1所示。圓盤中心與軸截面的中心剛性連接，這意味著圓盤的內邊界始終垂直于連接點處的軸。與轉軸類似，轉盤上任意一點位置矢量可寫成

C2rO2D=[rcosθrsinθwD]e2

式中,wD為極坐標r和θ的函數，代表轉盤的橫向位移。與式(2)同理，轉盤上任意一點的速度具體表達式為

(5)

1.3 能量表示

轉軸的動能方程保留二階及以下低階項可以表示為

(6)

轉盤的動能可化簡表示為

(7)

式中:Ri和Ro分別為轉盤的內外半徑,同時Ri也為轉軸的半徑;h為轉盤的厚度。從式(7)可知，第一部分表示剛性平移和轉動的動能，第二部分表示橫向振動的動能，最后一部分描述軸與盤之間的耦合效應。

Timoshenko轉軸的應變能表示為

(8)

式中:E為楊氏模量；其中к=6(1+μ)/(7+6μ)，к為剪切修正系數;μ為泊松比；G為剪切模量。薄環形盤由于橫向振動而產生的應變能為

(9)

此外，由于旋轉引起的應變能為

(10)

式中,σr和σθ為徑向方向和圓周方向上由于旋轉引起的初始應力。

總應變能可表示為VD=VD1+VD2。

1.4 能量表達的離散化

本節將采用有限元法和假設模態法[11]對系統進行離散。具有兩個4自由度節點的軸單元，如圖3所示。兩個平移自由度用沿x軸的uSi和沿y軸的vSi(i=1,2)表示，兩個旋轉自由度用繞x軸的θxi和繞y軸的θyi表示。因此，軸單元的位移矢量e可以表示為

圖3 Timoshenko轉軸單元和環狀轉盤單元

(11)

式中,ω為要確定的圓頻率。轉軸單元的橫向位移函數可以用形函數的形式表示為

(12)

式中:Nui(i=1,2,3,4)為轉軸橫向形狀函數單元;NS,u和NS,v為由Nui組成的轉軸橫向位移形函數。

此外，轉軸旋轉位移表示為

(13)

式中：Nφi(i=1,2,3,4)為轉軸旋轉位移形狀函數單元；NS,θx和NS, θy為由Nφi組成的轉軸旋轉位移形函數。

將形函數代入能量表達式，即可得到軸單元能量表達式的矩陣形式。動能表示為

(14)

其中,

(15)

應變能表示為

(16)

其中,

(17)

式中,le為軸單元的長度。

如圖3所示，用環單元來離散化轉盤。環單元中定義了3個節點，其中第一個是假定位于幾何中心的四自由度節點，同時也是軸單元的節點。另外兩個節點位于環單元同一條徑向線上的內邊界和外邊界，其中包含橫向位移wdi(d=s,c和i=2,3)和繞轉軸的旋轉角度φdi(d=s,c和i=2,3)。環單元第一個節點的位移矢量為

(18)

式中,下標“SD”為將軸與轉盤連接的耦合節點。包含余弦和正弦項在內的其他兩個節點的位移矢量為

(19)

環單元的橫向位移表示為

(20)

式中,n為節徑數目。NDw為環單元的形函數，表示為

NDw1=1-3ξ2+2ξ3,NDw2=(Ro-Ri)ξ(ξ-1)2,

NDw3=3ξ2-2ξ3,NDw4=(Ro-Ri)ξ2(ξ-1),

ξ=(r-Ri)/(Ro-Ri)

(21)

為了方便地組裝單元矩陣，將3個節點的位移矢量組合在一起，表示為

(22)

環單元內節點可表示為

(23)

環單元的橫向位移可表示為

(24)

轉盤動能表達的矩陣形式為

(25)

其中,

(26)

(28)

環單元的應變能表示為

(29)

其中，

(30)

1.5 軸-盤轉子系統的運動方程

將能量表達式寫成形函數乘以位移矢量的形式后，應用拉格朗日方程導出有限元法中軸單元和環單元的單元矩陣

(31)

式中，L=T-U，T為動能，U為應變能(i=S或者D，指該參數為轉軸或者轉盤對應的參數)。

組裝軸單元矩陣和環單元矩陣的示意圖，如圖4所示。通過將所有軸單元和環單元組裝在一起，可以得到多盤轉子系統的微分方程組

圖4 軸單元和環單元矩陣組裝到總體矩陣的示意圖

(32)

式中:M，G，K分別為系統的質量、陀螺、剛度矩陣；δ為系統的位移向量。

2 結果與分析

在本章節中，將使用幾個例子來分析軸-盤轉子系統的耦合自由振動特性。根據第1章的柔性軸-盤轉子動力學模型，柔性軸-盤轉子系統的幾何結構可以由軸長LS，m，軸半徑Ri，m，轉盤外半徑Ro，m，轉盤厚度HD，m和軸上的轉盤位置LD，m等參數來定義。材料參數包括楊氏模量E，N/m2，密度ρ，kg/m3和泊松比μ。本文中暫不考慮材料參數的影響，因此在算法程序上和ANSYS軟件仿真時采用相同的固定的材料參數。

在2.1節中，首先研究了數值模擬的計算收斂性，如表1所示。再用ANSYS軟件進一步驗證模型，ANSYS軟件建模過程如下：材料參數設置為與算法程序相同的固定值(E=2×1011N/m2,ρ=7 700 kg/m3,μ=0.3)。定義轉軸單元類型為solid95，定義轉盤單元類型為shell181，定義shell181單元相對應的轉盤厚度HD,m。繪制軸-盤轉子結構實體后進行自定義的網格劃分,如圖5和圖6所示。將轉軸和轉盤連接處的網格進行節點耦合，實現轉軸與轉盤的耦合。ANSYS整體模型，如圖7所示，施加與算法程序相同的邊界條件后，用ANSYS軟件的分塊蘭索斯求解方法求解零轉速下的軸-盤轉子模型的模態，得到表2和表3中ANSYS固有頻率計算結果；用ANSYS中QR阻尼求解方法，考慮陀螺效應，求解不同轉速下的軸-盤轉子模型的模態，得到表4中ANSYS固有頻率計算結果。

圖5 轉軸網格劃分的軸端視圖

圖6 ANSYS軸-盤網格單元節點耦合以及轉盤網格劃分示例

圖7 軸-盤轉子系統ANSYS模型局部視圖

圖8～圖15的結果都是無量綱的固有頻率，如果沒有另外說明，軸的邊界條件都是兩端簡支固定。此外，本文中只考慮轉盤模態為一節徑型模態的結果。

2.1 收斂與驗證

首先，研究了所提出數值模型的收斂性，如表1所示。其中ri=0.01，rD=0.2，hD=0.01，lD=0.5。可以看出，即使在4-6-6網格的情況下，前4個固有頻率也收斂到一個穩定值。說明了算法程序不會因網格數量的改變而引起固有頻率計算結果的明顯差異，驗證了程序的有效性。

其次，對于單盤轉子系統，軸一般分為兩段，每一段的單元數量見表1，為了進一步驗證所提出模型的準確性，使用另一個不同參數的例子(LS=0.46 m，Ri=0.01 m，LD=0.18 m，Ro=0.125 m，HD=0.006 m)，在表2中給出自由邊界條件下的有限元結果的比較，在表3中給出了表2例子簡支固定邊界條件下與有限元結果的比較，以及與Heydari等研究的例子結果的比較。從表2和表3可知，ANSYS結果驗證了Heydari等所提出計算方法的有效性，同時驗證了本文所用到計算程序的有效性。

表1 在簡支固定的邊界條件下，軸-盤轉子系統的前8階耦合無量綱固有頻率的收斂性

表2 兩端自由邊界條件下單盤轉子系統的固有頻率與ANSYS結果的比較

表3 簡支固定邊界條件下單盤轉子系統的固有頻率與ANSYS結果和參考文獻結果的比較

在表4中，運用與表2相同的幾何參數，將不同轉速下的固有頻率與ANSYS的固有頻率進行了比較。從表4可知，隨著轉速的增加，結果變化趨勢類似，考慮到兩者都是計算機數值仿真結果，兩個結果之間的誤差是可以接受的，從而驗證了本文所使用繪制坎貝爾圖程序的有效性。

表4 不同轉速下系統固有頻率本文計算結果(上)與ANSYS分析結果(下)的比較

2.2 軸-盤轉子系統的參數分析

在本節中，將研究ri，rD，hD和lD4個幾何參數對軸-盤轉子系統的影響。遵從實際情況，這些參數的范圍并非無限。考慮到軸始終是細長軸且轉盤始終是薄盤這個前提，幾個參數的大致變化范圍是：0.005≤ri≤0.050，0.1≤rD≤1.0，0.005≤hD≤0.050，0≤lD≤0.5，此外，本文中采用Timoshenko梁理論和經典薄板理論，不適用于短而厚的軸和厚圓盤。可以發現，ri和hD的范圍相比rD和lD的變化范圍小得多，為了結果更加合理，主要研究rD和lD這兩個參數的影響。

首先，設置ri=0.010，hD=0.010，令rD從0.05～1.00變化，令lD從0～0.5變化。軸-盤轉子系統的前6個耦合無量綱固有頻率λ相對于lD和rD的趨勢，如圖8所示(圖8中每一個小圖右側縱坐標代表λ的數值大小)。一般來說，隨著轉盤外半徑rD的增大，λ逐漸減小，而隨著rD的增大，λ逐漸增大和減小，這在高階固有頻率上表現得更為明顯。對于圖8中的第1階無量綱固有頻率，λ隨rD的增大而迅速減小，但是下降趨勢較圖8中其他曲線圖平緩，且不存在固有頻率的升降突變現象。然而，對于高階固有頻率，如第4階～第6階無量綱固有頻率，很容易區分頻率逐步衰減現象和升降突變現象。為更好研究頻率升降突變現象，以圖8中的第6階無量綱固有頻率圖為例，當0.13≤rD≤0.20時，由圖中高亮線條描繪可看見有3個凸型高峰和3個凹型低槽。當0.3≤rD≤0.5時，有3個高峰和兩個低槽。當0.5≤rD≤0.9時，雖然不明顯但可以確定有兩個高峰和兩個低槽。此外，lD在數值固定情況下，rD在一定的范圍內(0.13≤rD≤0.20，0.3≤rD≤0.5，0.5≤rD≤0.9)變化，無量綱固有頻率基本保持不變，因此可以用一條垂直等頻線來描述此規律；此外，對于固定值的rD，lD在某些間隔內變化時頻率也保持不變，因此用水平等頻線來描述此規律。為了更好地理解等頻線，圖9中顯示了幾條垂直等頻線和水平等頻線。

例如從圖9(a)圖中實線所示，對應圖8中第5階無量綱頻率，當rD從0.40變化到約0.62時，即垂直等頻線部分，無量綱頻率保持不變。

(a) 第1階無量綱頻率

接下來將考慮參數ri和hD對固有頻率的影響。由于兩個參數的變化范圍相對較窄，因此僅選擇部分指定值進行參數分析。圖10顯示出了當ri=0.010時第4階～第6階無量綱固有頻率，其中上方3個圖hD=0.005，下方3個圖hD=0.050。顯然易見，對于hD=0.005的較薄轉盤，垂直等頻線的長度更短；反之，對于hD=0.050的較厚轉盤則更長。圖11顯示出了當hD= 0.010時第4階～第6階無量綱固有頻率，其中上3個圖ri=0.005，下3個圖ri=0.050。可以發現，對于ri=0.005的細長軸，垂直等頻線的長度更長；反之，對于ri=0.050的較粗軸，其長度更短。特別注意到，當ri=0.050時，垂直等頻線幾乎消失并且僅能找到水平等頻線。可以得出結論，hD和ri對垂直等頻線的長度有相反的影響。圖10、圖11中每一個小圖右側縱坐標代表λ的數值大小。又如圖9(b)中緊靠于縱坐標數值3下方的虛線，對應圖11中下方小圖，描繪了rD=0.150的水平等頻線部分，無量綱頻率不會隨著相對轉盤位置的變化而變化。

(a) lD=0.50垂直等頻線

(a) 第4階無量綱頻率

(a) 第4階無量綱頻率

接下來，將討論當參數位于等頻線上時的軸-盤轉子系統。參考圖9(a)密集點虛線等頻線部分(hD=0.050，ri=0.010)，在0.1≤rD≤0.5間隔中幾乎是一條水平線。然后，分析具有相同參數(hD=0.050，ri=0.010，lD=0.500)軸-盤轉子系統的3個例子，變量是轉盤外徑rD，其中rD1=0.2，rD2=0.3，rD3=0.4。這3個例子的坎貝爾圖在圖12中給出，其中僅給出了第5階無量綱固有頻率用于比較。可以發現，不僅固有頻率幾乎相同，而且任意轉速下的正進動和反進動頻率也幾乎相同。對于圖9 (b)中給出的(rD=0.45，hD=0.010，ri=0.005)水平等頻線，可以看到這條底線在0.1≤lD≤0.5的間隔內幾乎是水平的。同理，選擇3種與圖9(b)具有相同參數的軸-盤轉子系統，變量為lD(lD1=0.2，lD2=0.3，lD3=0.4)。這3個例子的坎貝爾圖在圖13中展示，其中僅給出了第5階無量綱固有頻率用于比較。可以發現，隨著轉速的增加，無量綱頻率彼此分離，這表明，盡管這些例子的第5階自然頻率是相同的，但是它們在不同轉速下的頻率是不同的。因此，頻率隨轉速的變化不再相同。

圖12 (hD=0.050，ri=0.010，lD=0.5)軸-盤轉子系統在不同rD下的第5階無量綱固有頻率的坎貝爾圖

圖13 (rD=0.45，hD=0.010，ri=0.005)軸-盤轉子系統在不同lD下的第5階無量綱固有頻率的坎貝爾圖

參考圖10中下方中間的頻率圖，分析在垂直等頻線的基礎上參數rD的變化對于相應軸-盤轉子系統模態的影響，圖14中給出了3個樣例的第5階模態圖，結果表明在垂直等頻線上參數rD的變化并不能顯著改變軸-盤轉子系統的模態。同樣地，參考圖11(b)、圖11(e)頻率圖，分析在水平等頻線的基礎上參數lD的變化對模態的影響，圖15中給出了3個樣例的第5階模態圖，結合圖11分析可知，雖然在水平等頻線上參數lD的變化并未改變固有頻率，但是卻影響了軸-盤轉子系統模態圖中轉軸模態的形狀。

(a) rD=0.2

(a) lD=0.2

3 結 論

(1) 建立了柔性軸-盤轉子系統的轉子動力學模型，參數分析展示了柔性軸-盤轉子系統固有頻率逐步衰減現象和升降突變現象，采用“等頻線”描述柔性軸-盤轉子系統固有頻率升降突變規律，研究了參數hD和ri的變化對等頻線的影響，并進一步研究了在等頻線基礎上參數rD和lD的變化對系統振動特性的影響。

(2) 若在垂直等頻線上參數rD發生變化，柔性軸-盤轉子系統的固有頻率幾乎保持不變。若在水平等頻線上參數lD發生變化，對應的固有頻率也幾乎保持不變。

(3) 參數hD的增加會凸顯上文提及的頻率突變現象，也就是延長了垂直等頻線的長度，而參數ri則有相反的效果。

(4) 在垂直等頻線的研究基礎上，對于參數rD的不同的柔性軸-盤轉子系統樣例，坎貝爾圖展示了不僅零轉速下的固有頻率相同，而且在不同轉速下也始終相同。在水平等頻線的研究基礎上，而對于參數lD不同的樣例，雖然零轉速下相同，但是非零轉速下則彼此不同。

(5) 在垂直等頻線上(lD保持不變)，參數rD的變化并不能顯著改變柔性軸-盤轉子系統的模態。在水平等頻線上(rD保持不變)，參數lD的變化則影響了柔性軸-盤轉子系統模態圖中轉軸模態的形狀。