一種高效的電主軸系統復頻率計算方法

毛文貴，劉桂萍，韓 旭，郭維祺，林祿生

（1．湖南大學 機械與運載工程學院，長沙 410082；2．湖南工程學院 機械工程學院，湘潭 411101）

隨著當代高精密加工的飛速發展，對數控機床加工性能的要求也越來越高，而作為數控機床主要功能部件的電主軸單元是其性能提高的關鍵。電主軸以液體滑動軸承為支撐，將機床主軸和高速電動機功能從結構上融為一體，實現了變頻電機和機床主軸之間的“零傳動”。但高速電主軸技術還待改進，結構設計、工藝、甚至控制問題都未得到有效解決［1－3］。液體滑動軸承支承的高速電主軸可以避免滾動軸承較大的振動與噪聲、磁力軸承過于復雜的控制系統和氣體軸承較低的承載能力。但液體滑動軸承油膜的動特性對整個電主軸系統的動特性有很大影響，因此當前高速電主軸技術研究的一大熱點就是對滑動軸承油膜動特性的研究。滑動軸承油膜動特性通常可用8個系數來描述。Lund［4］就采用4個剛度系數和4個阻尼系數建立了撓性轉子系統模型，揭示了油膜渦動失穩的物理本質。對電主軸各向異性軸承－轉子系統進行動力特性分析時，往往要同時考慮轉子在徑向兩平面內的耦合運動以及軸承的阻尼，此時，系統內各個振動量之間發生了相位差，振動幅值隨時間增加或減小。因此，振動量為一復數，特征值為復頻率，它反映了軸承－轉子系統的穩定性，通常作為評價轉子性能的重要依據。Riccati傳遞矩陣法［5］因可消除奇點的干擾而在轉子動力學的振動計算中被廣泛采用。Murphy等［6］將復變量引入傳遞矩陣法，計算出系統的復頻率。陸頌元［7］對傳遞矩陣－多項式方法算法作了改進，只用一次遞推就求得特征多項式，大大節省了計算時間。故目前數值法求解復頻率一般采用傳遞矩陣一多項式算法［8］。但電主軸不同于一般主軸，電主軸的轉子電流和氣隙磁場相互作用產生電磁轉磁，使轉子沿旋轉磁場方向轉動，所以在對電主軸進行動力學分析時不能忽略轉子部分受到的電磁轉矩。電主軸是電機轉子與主軸過盈熱裝一起的軸，在進行動力學分析時應將其作為一個整體，電機的自重不能忽略。本文對Riccati傳遞矩陣法進行改進，獲得應用于電主軸各向異性軸承－轉子系統的傳遞矩陣，通過傳遞矩陣法直接獲得特征多項式方程隱式表達，即得到轉子系統的頻率方程式，由于這是個具有復數自變量S的復數方程式，因而在計算時需對衰減指數λ和阻尼圓頻率Ω兩個變量作“平面域”掃描，對此問題至今沒有簡便的計算方法。目前可以采用的是拋物線法（Muller法）求復方程根［9］，但Muller法中設定的起始試算頻率不同，求得各個根的先后次序并沒有一定的規律。求解所有的根對于大自由度系統而言工作量太大而不實用，而工程中我們主要關心產生共振的臨界轉速，一般只需知道前幾階。本文由幅角原理［10］判定系統共振區域中根的個數，采用雙種群進化策略［11］求解電主軸系統特征多項式的根，可快速有效地獲得電主軸系統的復頻率。

1 Riccati傳遞矩陣法的改進

1.1 電主軸系統的模塊化處理

電主軸系統被離散為N個集總圓盤，其間用N－1個無質量的彈性軸段相連，由于考慮了油膜力及其他各種阻尼，支承就是各向異性的。所以在X－Z，Y－Z平面（Z為軸向）同時彎曲，則狀態向量的維數是8，為提高傳遞矩陣的數值穩定性，把狀態向量中的元素分為兩部分，一部分為對應于初始截面狀態向量中具有零值的元素組成，另一部分為其余互補元素組成。考慮轉子的材料與主軸不同，在計算用于傳遞參數的轉動慣量等參數時會造成一定的麻煩和誤差，轉子的重力等效為作用于轉子軸段上的外力。同時電主軸轉子部分受到的電磁轉矩不能忽略，等效成一個集中轉矩施加于轉子軸段。因此，電主軸的傳遞矩陣法須考慮外力的加入，為使傳遞矩陣具有通用性并方便建立整體矩陣，本文將狀態向量的維數擴為10。狀態向量變為

式中：Mx，My為彎矩，Qx，Qy為剪力，X，Y為位移，θx，θy為轉角。向量中各元素均為復變量。如y＝Yest，其中Y＝Yc＋i Ys稱為該振動量的復振幅；Yc，Ys分別稱為它的余弦分量和正弦分量，S＝λ＋iΩ稱為復頻率，λ為衰減指數，Ω為阻尼圓頻率。

對于圖1中集總圓盤i和軸段i左右兩端進行受力分析得：

式中：I是單位矩陣，O是零矩陣，上標L、R分別表示集總圓盤和軸段的左端或右端。

式中：m為各節點的質量；kxx，kxy，kyx，kyy，Cxx，Cxy，Cyx，Cyy為支承軸承的剛度和阻尼系數；Jd，JP為直徑和極轉動慣量；ω為自轉角速度；V為考慮軸段的剪切影響，V＝6EJ／aGAL2；a為截面系數（實心圓軸為 0．886，空心圓軸為2／3）；G為剪切彈性模量；E為彈性模量；A為截面面積；L為軸段長；J為軸段的截面慣性矩。

1.2 特征多項式

引入Riccati變換并經推導，展開式（2）可得：

參考 Riccati推導［5］和陸頌元的多項式法［7］，引入表征同一截面從位移向量到力向量的轉換關系矩陣B，在每一截面上有：

式（4）式代入（3）式整理得：

引入

其中：Ti－1＝Ui－1＋Ri－1

整理可得狀態向量之間的遞推關系式：

引入邊界條件，如具有N個節點的轉子，兩端自由，則左節點的左邊及右節點的右邊兩截面力向量為零，可得：＝0＝0，由BL1＝0

可得＝0且C1＝I，D1＝W1，由＝0可得：

從i＝2起反復使用式（7），遞推到 i＝N得到DN（S）特征多項式。采用了4×4階矩陣代替了10×10階矩陣，從而機時進一步減少，具有計算速度快，占用內存少。解此特征多項式便得復頻率S＝λ＋iΩ，其中實部為一系統阻尼，反映轉子軸承系統的穩定性。對數衰減率為δ＝－2πλ／Ω。對數衰減率大于零則系統穩定。特征多項式的最高次數為整個系統的自由度數。對于節點數為N的轉子，有N×4對共扼復根。

2 高效電主軸系統復頻率計算方法

當軸承是各向異性軸承時，S是復數，則特征多項式是復函數，當轉子節點數較多時，易出現數值溢出問題，雖然張才穩提出特征多項式系數放大gi倍的方法使系數處于機器數系所能表達的范圍內［12］。g的取值沒有界定，g的尋找需要大量的時間，而且增加了程序的復雜性。現實中并不實用。工程中我們只對系統的某些低階特征值和特征向量感興趣，并不需要求解所有的根，對于臨界轉速，主要是了解與工作轉速存在共振可能的那幾階臨界轉速。因此計算系統在某一范圍內的部分特征值和特征向量便是一種非常實用而又經濟的方法。Muller法求解復函數，對初值選取較敏感。而且復數的求解不同于實數，由于設定的起始試算頻率不同，求得各個根的先后次序并沒有一定的規律，也不能確保在共振區域中所有的根都已經被找到。因此能了解在求解的范圍里的根的個數才能確保所有的根已經被找到，而不會漏根。以前實際計算中，采用先不計算阻尼和支承的各向異性，獲得剩余量的曲線，初步了解此系統的固有頻率的分布，然后求解復頻率［9］，這樣的過程過于復雜和耗時。本文基于復函數方程根的分布理論－幅角原理［13］判定共振區域中根的個數，在相應的區域提出基于雙種群進化策略求解特征多項式DN（S），獲得電主軸系統的復頻率。其流程如圖2所示。

圖2 高效電主軸系統復頻率計算方法流程Fig．2 A high efficient method for calculating complex frequency of the motorized spindle

2.1 復函數方程根的分布理論

定理 幅角原理［13］

若f（z）在正向簡單閉曲線C上每點都非零解析，在C的內部是亞純函數，那么f（z）在C內零點的個數等于1／（2π）乘以當 z沿 C的正向繞行一周 f（z）的幅角的改變量。

推論 設f（z）＝u＋iv在簡單閉曲線C上和在C內解析，且在C上不等于零，點z0＝x0＋iy0沿C的正向繞行一周，設向量（u，v）作正方向旋轉的次數為N1，作負方向旋轉的次數為N2，那么在封閉曲線C內f（z）＝0的根的個數N＝N1－N2。

通過幅角原理對所定義的區間進行根的個數的確定。在本文中共振區域內根的個數為電主軸的特征多項式的值在原點右邊的圈數。

2.2 雙種群進化策略尋根

進化策略是一種全局優化搜索算法，更新進化很快，特別適合函數優化［14］。在確定了編碼方案，適應度的計算及相應進化算子后，以“優勝劣汰”的方式，不斷逼近最優解。算法中主要用變異算子作為進化算子，變異算子主要有高斯變異和柯西變異，較大的高斯變異算子可以保證種群具有較好的探索能力，較小的高斯變異算子可以保證種群在局部具有良好的搜索能力。柯西分布因有一條很長的尾巴，容易產生一個遠離原點的隨機數，比高斯變異產生的隨機數具有更寬的分布范圍，這樣有可能很快跳出局部極小的區域。雙種群進化策略使得進化分別在兩個不同的種群間并行進行，兩個種群中采用不同的變異算子，種群1的目標變量與決策變量均采用高斯變異，種群2的目標變量采用柯西變異、決策變量采用高斯變異。種群1和種群2均采用（μ，λ）選擇策略，得出優良的個體作為進化的結果。雙種群方法具有計算精度高、自適應強而在求復函數方程中得到較好的應用。

根的個體由目標變量和標準差兩部分組成，每個個體有2個分量，分別代表根的實部和虛部。種群1和2的解代入特征多項式中，獲得其函數值，函數值越接近零，則根越接近最優解。規定適應度為最大值，取適應度函數為：

并終止條件為一個很接近1的數ε。當適應度值大于ε時終止，從而求解出一個復頻率。

2．3 綜合除法降次

用雙種群進化策略算法求出特征多項式一個根后，考慮轉子特征多項式根的共軛性，用綜合除法進行降次，方程（8）變為如下：

式中：Sr為已求根；為已求根的共軛根；F為已求出的根的個數。

獲得下次計算的特征多項式，可知求一個根后特征多項式可降低二次。

3 計算實例

3.1 電主軸系統

根據上述方法，對圖3所示模型轉子進行了計算，電主軸系統最高工作轉速為1 000 rad／s，電機功率為35．3 kW，主軸段密度 ρ為7 650 kg／m3。彈性模量 E為210 GPa。電主軸上各軸長 Lc、軸內徑 dc，、軸外徑Dc和圓盤寬度Lp，圓盤外徑Dp參數如表1，本模型被分為9個節點，液體滑動軸承支承在第3，5節點上，兩支承為相同的圓柱軸承，其參數為：長徑比L／D＝0．5，間隙比ψ＝2．0‰，動力粘度 η＝0．018 4 N．s／m2。本文采用熊萬里提出的基于動網格模型的計算方法［15］獲得液體動靜壓軸承7個偏心率下的8個剛度阻尼系數。軸承的剛度和阻尼無量綱系數隨偏心率的變化關系如圖4。利用改進的Riccati傳遞矩陣法獲得主軸段參數和節點總圓盤參數，得到包含復頻率的特征多項式的隱式表達。

圖3 電主軸各向異性軸承—轉子系統模型Fig．3 Anisotropic bearings-rotor systems of motorized spindle

表1 電主軸系統參數Tab．1 Parameters of motorized spindle system

圖4 軸承剛度阻尼無量綱系數隨偏心率的變化關系Fig．4 Dimensionless stiffness and damping coefficients vs．Journal eccentricity

3.2 高效電主軸系統復頻率計算方法

因本模型最高工作轉速n為1 000 rad／s，對于在第一臨界轉速nc1以下工作的轉子，應使工作轉速小于0．75 nc1，對于在第一臨界轉速nc1以上工作的轉子，應使1．4 nc1＜n＜0．7nc2。由此確定系統關心的臨界轉速不超過 2 000 rad／s，此系統關心的共振區域定為［0，2 000］。對工作頻率所可能激振的共振區域采用幅角原理進行根的個數的確定，如圖5所示，可知在［0，2 000］中電主軸系統特征多項式有5個根。用雙種群進化策略求解復頻率。雙種群進化策略中μ＝20，λ＝140，終止條件ε＝0．999，5個復頻率的計算結果如表2所示。為了驗證結果，采用傳統的Muller法求得本模型的特征值。本模型為有9個節點的軸系，共有36個復頻率，Muller法求復頻率的無序性使得必須求出全部的36個復頻率，然后按復頻率虛部值從小到大排列，取出前6個復頻率，如表2所示。對比表2中兩種方法求解的5個復頻率，可以看出，共振區間上兩種方法的結果誤差很小，說明本文所提方法求解的結果是準確的。第五個復頻率的虛部小于2 000 rad／s，而第六個復頻率的虛部大于2 000 rad／s，說明幅角原理對共振區間里根的個數的判斷是對的。從圖6顯示的適應度值隨迭代次數變化的5條曲線圖可知迭代次數不超過45代，說明雙種群進化策略收斂速度很快。同時采用本文方法只須求5個復頻率就能獲得工程中所關心的臨界轉速，而采用傳統的Muller法要求36個復頻率（本模型求出全部解需要9分鐘）。因此采用本文方法可以大大節省工作量。

圖5 基于輻角原理獲得共振區內復頻率的個數Fig．5 Number of complex frequency among resonance region based on argument principle

表2 復頻率計算結果及時間對比Tab．2 Comparison of complex frequency and time consumed

圖6 適應度值隨迭代次數變化的5條曲線圖Fig．6 five fitness value and iterative number of times change tendency

4 結 論

（1）本文通過對Riccati傳遞矩陣法進行改進，考慮電機外力的加入，將狀態向量的維數擴為10。采用Riccati分塊矩陣遞推，對所求的解不賦以具體數值而將其作為未知數保存，一次遞推求出以此為變量的特征多項式方程。

（2）針對復頻率的平面域掃描求解，傳統的Muller法求解根的無序性。依據復函數方程根的分布理論，由幅角原理判定工程中所關心的共振區域內根的個數。

（3）本文中采用雙種群進化策略，進化分別在兩個不同變異算子的種群中并行進行。避免了大規模范圍搜索時易產生的局部極值問題，是一種全局優化搜索方法。計算精度高、自適應強。

（4）以某高速主軸系統為例，采用本文方法對其進行了復頻率求解，并與Muller法進行對比，表明本方法計算結果準確，收斂速度很快，只須求幾個復頻率就能獲得工程中所關心的臨界轉速，工作量大大減小。這對做出可靠的穩定性分析是很重要的。

［1］熊萬里，陽雪兵，呂浪，等．液體動靜壓電主軸關鍵技術綜述［J］．機械工程學報，2009，45（9）：2－18．XIONG Wan-li，YANG Xue-bing，LULang，et al．Review on key technology of hydrodynamic and hydrostatic highfrequency motor spindles［J］． Journal of Mechanical Engineering，2009，45（9）：2－18．

［2］孟杰，陳小安，合燁．高速電主軸電動機－主軸系統的機電耦合動力學建模［J］．機械工程學報，2007，43（12）：160－165．MENG Jie， CHEN Xiao-an， HE Ye． Electromechanical coupled dynamical modeling of high speed motorized spindle’s motor-spindle subsystems［J］． Journal of Mechanical Engineering，2007，43（12）：160－165．

［3］陳小安，單文桃，周進明，等．高速電主軸驅動控制技術研究綜述［J］．振動與沖擊，2013，32（8）：39－47．CHEN Xiao-an，SHAN Wen-tao，ZHOU Jin-ming，et al．A survey of drive and control technology for high-speed motorized spindle systems［J］．Journal of vibration and shock，2013，32（8）：39－47．

［4］Lund JW．The stability of an elastic rotor in journal bearing with flexible damped supports［J］．ASME Journal of Applied Mechanics，1965，32（4）：911－920．

［5］王正．Riccati傳遞矩陣法的奇點及其消除方法［J］．振動與沖擊，1987，（2）：74－78．WANG Zheng．Singular point and elimination method of Riccati transfer matrix method［J］．Journal of Vibration and Shock，1987，（2）：74－78．

［6］Murphy B T，Vance JM．An improved method for calculating critical speeds and rotor dynamic stability of turbo machinery［C］．Trans of ASME for Power，1983，（5）：141－145．

［7］陸頌元．轉子系統復頻率的傳遞矩陣－多項式計算法［J］．應用力學學報，1986，31（1）：73－82．LU Song-yuan．A method for calculating damped critical speeds and stability of rotor-bearing systems［J］．Journal of Applied Mechanics，1986，31（1）：73－82．

［8］Zu JW，Ji Z．An improved transfer matrix method for steadystate analysis of nonlinear rotor-bearing systems ［J］．Transactions of the ASME，APRIL，2002（124）：302－310．

［9］聞邦椿，顧家柳，夏松波，等．高等轉子動力學［M］．北京：機械工業出版社，2000．

［10］Beardon A F．Complex analysis［J］．The Argument Principle in Analysis and Topology．Chichester，U．K．：Wiley，1979．

［11］夏慧明，梁華，周永權．用雙種群進化策略算法求解復函數方程的根［J］．計算機工程與應用，2008，44（7）：78－81．XIA Hum-ming，LIANG Hua，ZHOU Yong-quan．Novel bigroup evolution strategy algorithm for solving complex functional equation ［J］． Computer Engineering and Applications，2008，44（7）：78－81．

［12］張才穩，瞿紅春，陸頌元．轉子系統復頻率算法中數值溢出問題的研究［J］．振動、測試與診斷，1998，18（3）：206－210．ZHANG Cai-wen，QU Hong-chun，LU Song-yuan．Research on the problem of overflow in calculating the damped critical speeds of rotor-bearing systems［J］．Journal of Vibration，Measurement＆Diagnosis，1998，18（3）：206－210．

［13］Rusu C，Astola J．Note on argument principle［J］．Ieee Signal Processing Letters，2004，11（10）：817－820．

［14］Lee Wei-ping，Chien Wan-jou．A novel differential evolution with co-evolution［J］．Journal of Computers，2011，6（3）：594－602．

［15］熊萬里，侯志泉，呂浪，等．基于動網格模型的液體動靜壓軸承剛度阻尼計算方法［J］．機械工程學報，2012，48（23）：118－126．XIONG Wan-li，HOU Zhi-quan，LULang，et al．Method for calculating stiffness and damping coefficients of hybrid bearings based on dynamic mesh model［J］．Journal of Mechanical Engineering，2012，48（23）：118－126．