吳 磊，于慎波，于言明，竇汝桐，翟鳳晨

(1.沈陽工業(yè)大學機械工程學院, 沈陽 110870；2.中國郵政集團公司遼寧省機要通信局, 沈陽 110000)

0 引言

永磁同步電主軸作為高速機床的核心部件，將電機和主軸組成一個整體，實現(xiàn)了機床主軸系統(tǒng)的高效傳動[1]。隨著永磁材料的不斷創(chuàng)新與發(fā)展，使得永磁同步電主軸在航空航天、機械加工和精密制造等領(lǐng)域得到了越來普遍的應用[2]。永磁同步電主軸轉(zhuǎn)子細長，且轉(zhuǎn)速較高，和其它轉(zhuǎn)子相比，其轉(zhuǎn)子系統(tǒng)接近臨界轉(zhuǎn)速的機率更大。另外，由于外部環(huán)境、負載的變化、加工誤差等原因會導致轉(zhuǎn)子發(fā)生動態(tài)偏心[3]，進而產(chǎn)生不平衡磁拉力。在正常運行過程中，電主軸轉(zhuǎn)子固有頻率和轉(zhuǎn)子不平衡力導致的強迫振動頻率相耦合，當轉(zhuǎn)速接近該階模態(tài)下的臨界轉(zhuǎn)速時，主軸將會發(fā)生多重共振，導致噪聲過大，嚴重時，甚至導致轉(zhuǎn)子不穩(wěn)定或損壞[4]。為了避免轉(zhuǎn)子在運轉(zhuǎn)中產(chǎn)生共振，影響電主軸正常運行，需要對轉(zhuǎn)子偏心引起的不平衡磁拉力以及轉(zhuǎn)子臨界轉(zhuǎn)速和穩(wěn)定性進行相關(guān)的動力學分析[5-6]。

文獻[7-9]在轉(zhuǎn)子系統(tǒng)的不平衡響應分析、臨界轉(zhuǎn)速和穩(wěn)定性等領(lǐng)域做了許多的研究，為旋轉(zhuǎn)機械的發(fā)展奠定了理論基礎(chǔ)，但主要集中在水、汽輪機，對永磁同步電主軸的研究涉及很少。目前求解臨界轉(zhuǎn)速常采用傳統(tǒng)Prohl傳遞矩陣法、Riccati傳遞矩陣法和有限元法等[10]，其中Prohl傳遞矩陣法應用在各項同性轉(zhuǎn)子系統(tǒng)，Riccati應用在各項異性轉(zhuǎn)子系統(tǒng)，但由于復雜的外部環(huán)境以及算法本身的缺陷，計算時會出現(xiàn)根的溢出和丟根的現(xiàn)象，導致臨界轉(zhuǎn)速計算不精確。有限元法計算轉(zhuǎn)子系統(tǒng)臨界轉(zhuǎn)速時，對于復雜轉(zhuǎn)子由于建模復雜，占用計算機內(nèi)存較大，計算速度較慢，在一定程度上限制了其應用范圍。因此，如何提高解析法計算轉(zhuǎn)子臨界轉(zhuǎn)速的精確度以及進行機電耦合的分析，成為了轉(zhuǎn)子動力學研究熱點。

為了解決上述問題，本文中提出了一種用加時間因t的方法構(gòu)建轉(zhuǎn)子傳遞矩陣，運用了雙重QR法[11]求解傳遞矩陣特征值，進而可求得轉(zhuǎn)子系統(tǒng)的臨界轉(zhuǎn)速，既提高了計算速度，又解決了根的上下溢出問題。同時運用編程計算了轉(zhuǎn)子運行中由于動態(tài)偏心而導致的不平衡磁拉力，使其與機械系統(tǒng)進行耦合，分析其對動力學特性的影響，對研究轉(zhuǎn)子系統(tǒng)動力學特性有著重要的參考價值。

1 永磁同步電主軸臨界轉(zhuǎn)速計算

1.1 轉(zhuǎn)子離散建模與傳遞矩陣的建立

根據(jù)文獻[12]中節(jié)點劃分原則，應用集總質(zhì)量法對所求轉(zhuǎn)子進行離散處理[8]。本文中，永磁同步電主軸轉(zhuǎn)子如圖1所示，將其簡化為25個軸段，26個節(jié)點(節(jié)點7和19為軸承支撐點)并在對應節(jié)點附加上相應的質(zhì)量和轉(zhuǎn)動慣量。

圖1 電主軸轉(zhuǎn)子總成模型

對于離散后的第i個單元，將i和i+1分別定義為其左右兩端的截面編號，根據(jù)轉(zhuǎn)子受力特點，考慮軸承支承為各項異性，將截面的x方向和y方向彎矩、剪力、位移和轉(zhuǎn)角組合成該截面的狀態(tài)向量，如式(1)所示：

(1)

式中，Mx、My、θx、θy、x、y、Qx、Qy分別為水平和垂直方向的彎矩、轉(zhuǎn)角、撓度和剪力。

引入以復變量表示的振動量w=τ+σi，其中，w為復渦動頻率；實部τ為衰減指數(shù)；虛部σ為阻尼圓周率。根據(jù)圓盤的受力分析，可得出加時間因子t的剛性薄圓盤的傳遞矩陣，如式(2)所示：

(2)

根據(jù)軸段左右截面狀態(tài)向量關(guān)系，可以得出等截面彈性軸的傳遞矩陣，如式(3)所示：

(3)

將狀態(tài)向量{Z}i分別寫成力向量{Z}g和位移向量{Z}h，根據(jù)邊界條件以及轉(zhuǎn)子系統(tǒng)各截面狀態(tài)向量的關(guān)系，可得出方程(4)：

(4)

其中，[Q] =[B]K[D]K[B]K-1[D]K-1…[B]1[D]1。

由方程(4)可得：

[Q12]{Z}h,1=0

(5)

由于{Z}h,1≠0時，上式有非零解的條件為：

|Q12|=0

(6)

由于式(6)中的元素為多項式，運用代數(shù)余子式法將其展開，頻率方程實際上是關(guān)于w的8K階一元高次方程，即：

e0+e1w+e2w2+···+e8K-1w8K-1+w8K=0

(7)

求方程(7)的根可以轉(zhuǎn)化為求其對應友矩陣全部特征值的問題。

電主軸轉(zhuǎn)子為細長階梯軸，為了準確模擬實際轉(zhuǎn)子工況，需要離散足夠的節(jié)點數(shù)，本文中離散節(jié)點數(shù)K=26，方程最高階會達到208，求解過程中會出現(xiàn)根的溢出，求解難度特別大。文中，為了縮小求解過程中大數(shù)與小數(shù)之間的差異，防止根的溢出，提出加時間因子t的方法建立轉(zhuǎn)子系統(tǒng)傳遞矩陣。計算過程中其他參數(shù)與時間因子相乘，轉(zhuǎn)速與時間因子相除，運用雙重步QR法編程求解出復頻率的實部與虛部，計算完成后將虛部乘以時間因子得到臨界轉(zhuǎn)速。

1.2 解析法計算永磁同步電主軸臨界轉(zhuǎn)速

本文中以額定功率9.5 kW永磁同步電主軸轉(zhuǎn)子為例，取時間因子t等于15 000構(gòu)建傳遞矩陣，取3000 r/min為轉(zhuǎn)速步長，運用編程的方法求解出復頻率并畫出坎貝爾圖。同頻線與各階渦動速度曲線的交點即是轉(zhuǎn)子系統(tǒng)臨界轉(zhuǎn)速，如圖2所示。從圖中可以看出，隨著轉(zhuǎn)子自轉(zhuǎn)速度逐步升高，一階正向和反向渦動速度相近且變化很??；二階反向渦動速度逐步減小，正向渦動速度逐步增大。

圖2 解析法計算轉(zhuǎn)子坎貝爾圖

1.3 有限元法計算永磁同步電主軸臨界轉(zhuǎn)速

利用ANSYS計算臨界轉(zhuǎn)速時，運用實體模型建模復雜，且計算結(jié)果精度達不到要求。因此本文運用有限元法時，以MASS21、BEAM188和COMBI214分別模擬質(zhì)量點、彈性軸和軸承單元。同時質(zhì)量點單元、軸承單元分別施加旋轉(zhuǎn)約束和全約束，并充分考慮了陀螺效應的影響，通過編寫ANSYS APDL命令流的方式建立二維有限元模型如圖3所示。

圖3 ANSYS二維質(zhì)量點單元模型

本文用QR阻尼法進行模態(tài)提取，利用PLCAMP命令畫出軸承-轉(zhuǎn)子系統(tǒng)的CAMP-BELL圖，如圖4所示。

圖4 ANSYS有限元計算轉(zhuǎn)子坎貝爾圖

1.4 解析法和有限元法計算臨界轉(zhuǎn)速結(jié)果對照

各階臨界轉(zhuǎn)速為同頻線與正向渦動速度的交點。轉(zhuǎn)子臨界轉(zhuǎn)速計算結(jié)果對比如表1所示。由表1可知，解析法與有限元計算結(jié)果對照，誤差均在3.2%以內(nèi)，驗證了本文中解析法計算臨界轉(zhuǎn)速的有效性。

表1 臨界轉(zhuǎn)速結(jié)果對比 (r/min)

2 穩(wěn)定性分析

穩(wěn)定性裕度體現(xiàn)在系統(tǒng)自身的修復能力，也就是能夠回到平穩(wěn)狀態(tài)的能力，度量方法主要有對數(shù)衰減率、系統(tǒng)阻尼和系統(tǒng)抗定常干擾界限值法，其中對數(shù)衰減率是以振幅衰減的快慢程度為指標，對系統(tǒng)受干擾后恢復到平衡狀態(tài)過程的度量[8]。

本文以對數(shù)衰減率為依據(jù)判定軸系的穩(wěn)定性裕度。由本文第一部分解析法得到復頻率w=τ+σi，則第i階特征值對應的對數(shù)衰減率：

(8)

γi越大表明系統(tǒng)穩(wěn)定程度越高，當γi<0時，表明系統(tǒng)處于失穩(wěn)狀態(tài)。

計算得到的軸系對數(shù)衰減率曲線如圖5所示。圖中表明，一階反向、二階正向和三階反向?qū)?shù)衰減率呈逐漸下降趨勢，其中三階反向?qū)?shù)衰減率最小且有加速下降的趨勢，最易引起系統(tǒng)失穩(wěn)，設計時應盡量避開三階反向渦動臨界轉(zhuǎn)速。

圖5 永磁同步電主軸轉(zhuǎn)子對數(shù)衰減率

3 不平衡磁拉力

3.1 不平衡磁拉力的計算

通常氣隙磁場分布均勻，轉(zhuǎn)子不承受不平衡磁拉力，但由于轉(zhuǎn)子偏心導致徑向氣隙磁場畸變，會產(chǎn)生不平衡磁拉力，使得轉(zhuǎn)子受力不均勻。傳統(tǒng)公式計算不平衡磁拉力時，將模型簡化為理想狀態(tài)，計算結(jié)果不夠精確。因此本文基于Ansoft有限元軟件對某型號8極36槽永磁同步電主軸轉(zhuǎn)子動態(tài)偏心引起的不平衡磁拉力進行分析。電主軸的基本參數(shù)如表2所示。

表2 永磁同步電主軸部分參數(shù)

電主軸有限元模型建立完成并求解后可提取氣隙磁密，然后由式(9)求出不平衡磁拉力[13-15]。

(9)

其中，F(xiàn)x、Fy分別是不平衡磁拉力在X和Y方向上的分量；l是轉(zhuǎn)子硅鋼片的軸向長度；rp是不平衡磁拉力積分半徑。

轉(zhuǎn)子偏心后，徑向畸變氣隙磁場可看作由多個同心機磁場疊加而成，且切向氣隙磁密對徑向不平衡磁拉力影響極小，可忽略不計。因此可以將氣隙積分圓周離散成N個相等的子區(qū)域，求得各區(qū)域徑向畸變氣隙磁場產(chǎn)生的不平衡磁拉力然后進行疊加，所以式(9)可以表示成式(10)的形式。

(10)

式中，δiω、δiρ分別為第i個區(qū)域的起始和結(jié)束角度。

則不平衡磁拉力的合力為：

(11)

為了防止轉(zhuǎn)子和定子出現(xiàn)掃膛現(xiàn)象，轉(zhuǎn)子偏心距不應大于平均氣隙的10%[16]。本文中電主軸氣隙總長為1.6 mm，因此偏心距不應超過0.16 mm。根據(jù)此原則，以無偏心和動態(tài)偏心0.16 mm為例，將模型積分圓周劃分為1000個區(qū)域，分別提取其徑向氣隙磁密，如圖6、圖7所示。

圖6 無偏心徑向氣隙磁密

圖7 偏心距0.16 mm徑向氣隙磁密

應用式(10)，利用編程求出不平衡磁拉力在X和Y軸的分力，如圖8所示。圖中可看出，當無偏心時，X和Y方向不平衡電磁力均在0 N附近上下波動，且浮動很??；當偏心距為0.16 mm時，X和Y方向上均產(chǎn)生了隨時間周期性變化的幅值為230 N的不平衡磁拉力。

圖8 不平衡磁拉力在X和Y軸的分力

應用上述方法，可求出偏心距分別為0.032 mm、0.064 mm、0.096 mm和0.128 mm時在X和Y方向的不平衡磁拉力，由式(10)、式(11)可求得不同偏心距下不平衡磁拉力在Y方向分量Fy和在X方向分量Fx隨圓周的變化以及不平衡磁拉力的合力，如圖9、圖10所示。

圖9 Fy隨Fx在圓周內(nèi)的變化

圖10 不平衡磁拉力合力隨時間變化

由圖9可看出，相同偏心距下，不同時刻的X和Y方向的不平衡磁拉力在圓周上組成了一個同心圓，且圓的半徑即為該偏心狀態(tài)下的不平衡磁拉力，這是由于動態(tài)偏心時最小氣隙隨著時間在變化。由圖9和圖10可得到不同偏心距下不平衡磁拉力的數(shù)值如表3所示，且不平衡磁拉力隨著偏心距的增大呈線性增加。

表3 不同偏心距下的不平衡磁拉力

3.2 機電耦合

當電主軸轉(zhuǎn)子在高速旋轉(zhuǎn)時，不平衡磁拉力被引入機械系統(tǒng)引發(fā)共振[17]，嚴重情況下會造成定轉(zhuǎn)子接觸，刮蹭永磁體。因此在動力學分析時，考慮不平衡磁拉力的影響十分必要。將不同動態(tài)偏心下產(chǎn)生的不平衡磁拉力與轉(zhuǎn)子機械系統(tǒng)進行耦合后，臨界轉(zhuǎn)速的變化如圖11所示。

圖11 不平衡磁拉力對臨界轉(zhuǎn)速的影響

通過分析圖11可以看出，考慮不平衡磁拉力時的臨界轉(zhuǎn)速低于無不平衡磁拉力時的臨界轉(zhuǎn)速。隨著不平衡磁拉力的增加，一階臨界轉(zhuǎn)速由23 657 r/min降低到23 339 r/min，降幅為1.4%，二階臨界轉(zhuǎn)速幾乎無明顯變化。

4 結(jié)論

(1)運用加時間因子t的方法構(gòu)建轉(zhuǎn)子傳遞矩陣，同時運用雙重步QR法計算轉(zhuǎn)子系統(tǒng)臨界轉(zhuǎn)速，既提高了計算精度又解決了根的上下溢出問題。同時將解析法與有限元法計算的各階臨界轉(zhuǎn)速進行對比，誤差均小于3.2%，驗證了解析法的有效性。

(2)通過分析對數(shù)衰減率曲線，發(fā)現(xiàn)三階反向渦動衰減率曲線最低，且下降速度有逐漸加快的趨勢。說明三階反向渦動穩(wěn)定性裕度最小，最易引起轉(zhuǎn)子系統(tǒng)失穩(wěn)。

(3)分析了動態(tài)偏心工況下的不平衡磁拉力，結(jié)果表明，不平衡磁拉力隨著偏心距的增大呈線性增加。隨著不平衡磁拉力增加，一階臨界轉(zhuǎn)速將下降，二階臨界轉(zhuǎn)速幾乎不受影響。