一種高效的電主軸系統(tǒng)復(fù)頻率計(jì)算方法

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

（1．湖南大學(xué) 機(jī)械與運(yùn)載工程學(xué)院，長(zhǎng)沙 410082；2．湖南工程學(xué)院 機(jī)械工程學(xué)院，湘潭 411101）

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

1 Riccati傳遞矩陣法的改進(jìn)

1.1 電主軸系統(tǒng)的模塊化處理

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

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

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

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

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

1.2 特征多項(xiàng)式

引入Riccati變換并經(jīng)推導(dǎo)，展開式（2）可得：

參考 Riccati推導(dǎo)［5］和陸頌元的多項(xiàng)式法［7］，引入表征同一截面從位移向量到力向量的轉(zhuǎn)換關(guān)系矩陣B，在每一截面上有：

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

引入

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

整理可得狀態(tài)向量之間的遞推關(guān)系式：

引入邊界條件，如具有N個(gè)節(jié)點(diǎn)的轉(zhuǎn)子，兩端自由，則左節(jié)點(diǎn)的左邊及右節(jié)點(diǎn)的右邊兩截面力向量為零，可得：＝0＝0，由BL1＝0

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

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

2 高效電主軸系統(tǒng)復(fù)頻率計(jì)算方法

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

圖2 高效電主軸系統(tǒng)復(fù)頻率計(jì)算方法流程Fig．2 A high efficient method for calculating complex frequency of the motorized spindle

2.1 復(fù)函數(shù)方程根的分布理論

定理 幅角原理［13］

若f（z）在正向簡(jiǎn)單閉曲線C上每點(diǎn)都非零解析，在C的內(nèi)部是亞純函數(shù)，那么f（z）在C內(nèi)零點(diǎn)的個(gè)數(shù)等于1／（2π）乘以當(dāng) z沿 C的正向繞行一周 f（z）的幅角的改變量。

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

通過幅角原理對(duì)所定義的區(qū)間進(jìn)行根的個(gè)數(shù)的確定。在本文中共振區(qū)域內(nèi)根的個(gè)數(shù)為電主軸的特征多項(xiàng)式的值在原點(diǎn)右邊的圈數(shù)。

2.2 雙種群進(jìn)化策略尋根

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

根的個(gè)體由目標(biāo)變量和標(biāo)準(zhǔn)差兩部分組成，每個(gè)個(gè)體有2個(gè)分量，分別代表根的實(shí)部和虛部。種群1和2的解代入特征多項(xiàng)式中，獲得其函數(shù)值，函數(shù)值越接近零，則根越接近最優(yōu)解。規(guī)定適應(yīng)度為最大值，取適應(yīng)度函數(shù)為：

并終止條件為一個(gè)很接近1的數(shù)ε。當(dāng)適應(yīng)度值大于ε時(shí)終止，從而求解出一個(gè)復(fù)頻率。

2．3 綜合除法降次

用雙種群進(jìn)化策略算法求出特征多項(xiàng)式一個(gè)根后，考慮轉(zhuǎn)子特征多項(xiàng)式根的共軛性，用綜合除法進(jìn)行降次，方程（8）變?yōu)槿缦拢?/p>

式中：Sr為已求根；為已求根的共軛根；F為已求出的根的個(gè)數(shù)。

獲得下次計(jì)算的特征多項(xiàng)式，可知求一個(gè)根后特征多項(xiàng)式可降低二次。

3 計(jì)算實(shí)例

3.1 電主軸系統(tǒng)

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

圖3 電主軸各向異性軸承—轉(zhuǎn)子系統(tǒng)模型Fig．3 Anisotropic bearings-rotor systems of motorized spindle

表1 電主軸系統(tǒng)參數(shù)Tab．1 Parameters of motorized spindle system

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

3.2 高效電主軸系統(tǒng)復(fù)頻率計(jì)算方法

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

圖5 基于輻角原理獲得共振區(qū)內(nèi)復(fù)頻率的個(gè)數(shù)Fig．5 Number of complex frequency among resonance region based on argument principle

表2 復(fù)頻率計(jì)算結(jié)果及時(shí)間對(duì)比Tab．2 Comparison of complex frequency and time consumed

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

4 結(jié) 論

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

（2）針對(duì)復(fù)頻率的平面域掃描求解，傳統(tǒng)的Muller法求解根的無序性。依據(jù)復(fù)函數(shù)方程根的分布理論，由幅角原理判定工程中所關(guān)心的共振區(qū)域內(nèi)根的個(gè)數(shù)。

（3）本文中采用雙種群進(jìn)化策略，進(jìn)化分別在兩個(gè)不同變異算子的種群中并行進(jìn)行。避免了大規(guī)模范圍搜索時(shí)易產(chǎn)生的局部極值問題，是一種全局優(yōu)化搜索方法。計(jì)算精度高、自適應(yīng)強(qiáng)。

（4）以某高速主軸系統(tǒng)為例，采用本文方法對(duì)其進(jìn)行了復(fù)頻率求解，并與Muller法進(jìn)行對(duì)比，表明本方法計(jì)算結(jié)果準(zhǔn)確，收斂速度很快，只須求幾個(gè)復(fù)頻率就能獲得工程中所關(guān)心的臨界轉(zhuǎn)速，工作量大大減小。這對(duì)做出可靠的穩(wěn)定性分析是很重要的。

［1］熊萬里，陽雪兵，呂浪，等．液體動(dòng)靜壓電主軸關(guān)鍵技術(shù)綜述［J］．機(jī)械工程學(xué)報(bào)，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］孟杰，陳小安，合燁．高速電主軸電動(dòng)機(jī)－主軸系統(tǒng)的機(jī)電耦合動(dòng)力學(xué)建模［J］．機(jī)械工程學(xué)報(bào)，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ìn)明，等．高速電主軸驅(qū)動(dòng)控制技術(shù)研究綜述［J］．振動(dòng)與沖擊，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傳遞矩陣法的奇點(diǎn)及其消除方法［J］．振動(dòng)與沖擊，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］陸頌元．轉(zhuǎn)子系統(tǒng)復(fù)頻率的傳遞矩陣－多項(xiàng)式計(jì)算法［J］．應(yīng)用力學(xué)學(xué)報(bào)，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］聞邦椿，顧家柳，夏松波，等．高等轉(zhuǎn)子動(dòng)力學(xué)［M］．北京：機(jī)械工業(yè)出版社，2000．

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

［11］夏慧明，梁華，周永權(quán)．用雙種群進(jìn)化策略算法求解復(fù)函數(shù)方程的根［J］．計(jì)算機(jī)工程與應(yīng)用，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］張才穩(wěn)，瞿紅春，陸頌元．轉(zhuǎn)子系統(tǒng)復(fù)頻率算法中數(shù)值溢出問題的研究［J］．振動(dòng)、測(cè)試與診斷，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］熊萬里，侯志泉，呂浪，等．基于動(dòng)網(wǎng)格模型的液體動(dòng)靜壓軸承剛度阻尼計(jì)算方法［J］．機(jī)械工程學(xué)報(bào)，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．