柔性轉(zhuǎn)子滾動(dòng)軸承系統(tǒng)混沌行為研究

崔 立 鄭建榮

1.上海第二工業(yè)大學(xué)，上海，201209 2.華東理工大學(xué)，上海，200237

0 引言

滾動(dòng)軸承支承的轉(zhuǎn)子系統(tǒng)由于軸承非線性剛度而出現(xiàn)復(fù)雜的非線性動(dòng)力學(xué)特性。近年來，剛性轉(zhuǎn)子滾動(dòng)軸承系統(tǒng)的分叉及混沌特性研究得到了廣泛的重視。

Jang等［1］考慮滾動(dòng)軸承的5個(gè)自由度以及表面波紋度建立了剛性轉(zhuǎn)子滾動(dòng)軸承系統(tǒng)的非線性動(dòng)力學(xué)方程，發(fā)現(xiàn)套圈的波紋會(huì)導(dǎo)致徑向位移和角位移響應(yīng)峰值附近出現(xiàn)邊頻帶。Harsha［2－3］考慮軸承間隙、不平衡力研究了轉(zhuǎn)子系統(tǒng)的動(dòng)力學(xué)響應(yīng)，發(fā)現(xiàn)系統(tǒng)的動(dòng)力學(xué)行為與滾動(dòng)體的通過頻率有關(guān)，且當(dāng)通過頻率及其諧波與固有頻率相等、外圈波紋度階數(shù)與滾動(dòng)體數(shù)目相等時(shí)振動(dòng)幅值將變得極大。Bai等［4］研究了主軸－滾動(dòng)軸承系統(tǒng)非線性動(dòng)力學(xué)行為，發(fā)現(xiàn)軸承間隙的減小有利于提高主軸軸承系統(tǒng)的穩(wěn)定性。高尚晗等［5］研究了主軸－滾動(dòng)軸承在負(fù)游隙情況下的機(jī)床主軸－滾動(dòng)軸承系統(tǒng)的非線性動(dòng)力學(xué)特性，揭示了主軸系統(tǒng)的混沌演化過程。崔立等［6］研究了圓柱滾子軸承剛性轉(zhuǎn)子系統(tǒng)周期運(yùn)動(dòng)分岔特性，發(fā)現(xiàn)隨著徑向間隙、阻尼和力矩的變化，周期運(yùn)動(dòng)將產(chǎn)生倍周期或Hopf分岔，分岔轉(zhuǎn)速隨參數(shù)變化而改變。

以上的研究對象均為剛性轉(zhuǎn)子，模型涉及轉(zhuǎn)軸的彎曲變形和陀螺力矩等參數(shù)。隨著旋轉(zhuǎn)機(jī)械轉(zhuǎn)速的提高，柔性轉(zhuǎn)子系統(tǒng)的設(shè)計(jì)與分析變得越來越重要。近年也有一些研究柔性轉(zhuǎn)子軸承系統(tǒng)的文章發(fā)表，如：Laha等［7］研究了油膜軸承支承的柔性轉(zhuǎn)子系統(tǒng)的分叉行為，分析了轉(zhuǎn)軸的材料、剛度和質(zhì)量等參數(shù)對轉(zhuǎn)子系統(tǒng)分叉行為的影響。Villa等［8］、Sinou［9］研究了球軸承支承的柔性轉(zhuǎn)子系統(tǒng)的非線性動(dòng)力學(xué)行為，發(fā)現(xiàn)系統(tǒng)響應(yīng)中存在跳躍現(xiàn)象和超諧波，但其研究僅考慮了4個(gè)自由度的轉(zhuǎn)軸節(jié)點(diǎn)和2個(gè)軸承自由度，難以滿足實(shí)際工況的需求。

本文采用12自由度Euler－Bernoulli桿單元建立柔性轉(zhuǎn)子滾子軸承系統(tǒng)的非線性動(dòng)力學(xué)模型，研究系統(tǒng)的混沌行為，分析軸承結(jié)構(gòu)參數(shù)與轉(zhuǎn)軸結(jié)構(gòu)參數(shù)對轉(zhuǎn)子系統(tǒng)混沌行為的影響規(guī)律。

1 計(jì)算模型

圖1所示為柔性轉(zhuǎn)子軸承系統(tǒng)模型，采用12自由度的Euler－Bernoulli桿單元，圓盤簡化為質(zhì)點(diǎn)并考慮其質(zhì)量與轉(zhuǎn)動(dòng)慣量，轉(zhuǎn)子由兩個(gè)滾動(dòng)軸承支承。

圖1 柔性轉(zhuǎn)子滾動(dòng)軸承系統(tǒng)模型

1.1 柔性轉(zhuǎn)子軸承系統(tǒng)動(dòng)力學(xué)模型

柔性轉(zhuǎn)子系統(tǒng)的動(dòng)力學(xué)方程為

式中，M為包含轉(zhuǎn)軸、圓盤和軸承的總質(zhì)量矩陣；C為總阻尼矩陣；G為總陀螺矩陣；K為總剛度矩陣；X為轉(zhuǎn)子各節(jié)點(diǎn)的位移向量；f（X，t）為包括轉(zhuǎn)子各節(jié)點(diǎn)軸承力、重力、轉(zhuǎn)子不平衡力和外載荷的矩陣。

先對轉(zhuǎn)子系統(tǒng)進(jìn)行軸段與節(jié)點(diǎn)劃分，然后進(jìn)行轉(zhuǎn)子系統(tǒng)的質(zhì)量矩陣、剛度矩陣、阻尼矩陣、陀螺矩陣、載荷矩陣求解，并按照節(jié)點(diǎn)的順序?qū)Ω骶仃囘M(jìn)行組裝，具體過程如下：首先計(jì)算桿單元的質(zhì)量矩陣、陀螺矩陣、剛度矩陣、載荷矩陣并組裝；然后將剛性圓盤的質(zhì)量矩陣、陀螺矩陣、不平衡力矩陣疊加到所在節(jié)點(diǎn)的相應(yīng)矩陣中；之后將支承軸承的質(zhì)量矩陣、非線性軸承力疊加到所在節(jié)點(diǎn)的相應(yīng)矩陣中；最后計(jì)算系統(tǒng)的結(jié)構(gòu)阻尼、軸承阻尼矩陣并組裝。

1.2 桿單元

圖2所示為針對空間桿單元分析獲得的空間桿單元兩端節(jié)點(diǎn)位移。

圖2 二節(jié)點(diǎn)空間桿單元

針對圖2中12自由度Euler－Bernoulli桿單元，若每個(gè)節(jié)點(diǎn)考慮6個(gè)自由度，則共有6個(gè)廣義位移和6個(gè)廣義力，其表達(dá)式為

式中，xi、yi、zi分別為節(jié)點(diǎn)i沿x、y、z方向的線位移；θxi、θyi、θzi分別為節(jié)點(diǎn)i處截面繞x、y、z軸的轉(zhuǎn)角位移；Fxi、Fyi、Fzi分別為節(jié)點(diǎn)i在x、y、z方向受到的軸向力和剪切力；Mxi、Myi、Mzi分別為節(jié)點(diǎn)i在x、y、z方向所受的扭矩和彎矩。

假設(shè)已知桿單元橫截面面積、截面慣性矩、單元的扭轉(zhuǎn)慣性矩、長度、材料彈性模量和剪切模量，則可根據(jù)有限元理論［10］求出桿單元的質(zhì)量矩陣Ms、剛度矩陣Ks、阻尼矩陣Cs、陀螺矩陣Gs。

1.3 剛性圓盤模型

當(dāng)軸上安裝有圓盤時(shí)，將其視為剛性圓盤，并將其質(zhì)量矩陣、陀螺矩陣、不平衡力矩陣疊加到所在節(jié)點(diǎn)的相應(yīng)矩陣中。

假設(shè)已知圓盤的質(zhì)量、半徑、不平衡質(zhì)徑積，則可建立剛性圓盤的質(zhì)量矩陣Md、陀螺矩陣Gd、不平衡力矩陣Fd，其中Fd為

式中，me為圓盤的不平衡質(zhì)徑積；ω為轉(zhuǎn)速。

1.4 滾子軸承模型

考慮普遍受載的圓柱滾子軸承，其模型如圖3所示，圖中，Dr為滾子直徑，D1、D2分別為外圈外徑和內(nèi)圈內(nèi)徑，δ0為徑向間隙，le為帶凸度滾子的長度，ls為帶凸度滾子直線部分的長度。

假設(shè)滾子數(shù)目為N，使用切片法將滾子分成nr個(gè)圓片。對第j個(gè)滾子進(jìn)行受力分析，假設(shè)第j個(gè)滾子方位角為φj，根據(jù)赫茲接觸理論，并使用擬動(dòng)力學(xué)方法建立滾子軸承的非線性平衡方程組，使用Newton－Raphson法可求出滾子與套圈的接觸力［11］。將各滾子的接觸力分解，即可求出圓柱滾子軸承的非線性軸承力矩陣：

圖3 圓柱滾子軸承結(jié)構(gòu)簡圖

式中，F(xiàn)2jk為第j個(gè)滾子的第k個(gè)切片與滾子軸承內(nèi)圈的作用力。

1.5 結(jié)構(gòu)阻尼求解

Rayleigh提出的結(jié)構(gòu)阻尼模型計(jì)算表達(dá)式為

式中，ε1、ε2為兩個(gè)振型的阻尼比，根據(jù)經(jīng)驗(yàn)取ε1＝0.005，ε2＝0.01；ω1、ω2為轉(zhuǎn)子系統(tǒng)的二階固有頻率。

式（7）中轉(zhuǎn)子系統(tǒng)的固有頻率ω1、ω2可根據(jù)計(jì)算得到的質(zhì)量矩陣、剛度矩陣解｜K－Mω2｜＝0得到。

2 計(jì)算方法

采用Runge－Kutta法、Newton－Raphson法進(jìn)行非線性動(dòng)力學(xué)方程組求解，根據(jù)FPA修正法確定求解周期，然后求解最大Lyapunov指數(shù)，判斷系統(tǒng)的動(dòng)力學(xué)行為。

2.1 求解周期的確定

轉(zhuǎn)子軸承系統(tǒng)中存在軸承變剛度激勵(lì)，還可能存在不平衡力激勵(lì)。不平衡力產(chǎn)生的激勵(lì)周期為軸轉(zhuǎn)動(dòng)周期的整數(shù)倍，但軸承的變剛度激勵(lì)周期往往不是軸轉(zhuǎn)動(dòng)周期的整數(shù)倍，所以在判斷和求解時(shí)，采用修正的FPA法建立統(tǒng)一的求解周期［12］，其表達(dá)式定義為

式中，T為求解周期；Td為不平衡力激勵(lì)周期；TVC為軸承變剛度激勵(lì)周期；ε為常數(shù)，取ε＝0.01；K為比例系數(shù)，K ＝1，2，…，nk。

根據(jù)式（9）進(jìn)行循環(huán)計(jì)算，直至找到滿足其要求的K值，代入式（8）即得求解周期。

2.2 Lyapunov指數(shù)計(jì)算

Lyapunov指數(shù)表示在相平面中2條相鄰軌線間的距離隨時(shí)間的平均指數(shù)發(fā)散率，它明確地區(qū)分了確定性運(yùn)動(dòng)和混沌運(yùn)動(dòng)［13］。

對于連續(xù)系統(tǒng)有

設(shè)在τ0時(shí)刻‖δX（τ0）‖充分小，于是1維的Lyapunov指數(shù)可定義為

在n維連續(xù)系統(tǒng)中，δX（τ）在每個(gè)基底上有分量，每一個(gè)分量均可按上式求出一個(gè)λ，因此共存在n個(gè)Lyapunov指數(shù)λi，稱為Lyapunov指數(shù)譜。當(dāng)任意選取矩陣δX（τ）時(shí)，Lyapunov指數(shù)以概率1可能取得最大值，如果其中最大的Lyapunov指數(shù)λmax＞0，則該系統(tǒng)一定存在混沌運(yùn)動(dòng)。因此，只要計(jì)算出系統(tǒng)的最大Lyapunov指數(shù)，就可以判斷系統(tǒng)是否處于混沌狀態(tài)。

3 計(jì)算與分析

圖4為柔性轉(zhuǎn)子系統(tǒng)簡圖，轉(zhuǎn)軸由兩個(gè)型號(hào)相同的滾子軸承支承，轉(zhuǎn)軸中有剛性圓盤，該圓盤可施加不平衡力。

圖4 柔性轉(zhuǎn)子系統(tǒng)簡圖

表1所示為滾子軸承的結(jié)構(gòu)參數(shù)，其中，軸承的彈性模量為204GPa，泊松比為0.3，阻尼為200N·s／m，軸承載荷為｛0，2000N，2000N，0，0｝。

表1 滾子軸承的結(jié)構(gòu)參數(shù)

柔性轉(zhuǎn)子系統(tǒng)的結(jié)構(gòu)參數(shù)如表2所示，轉(zhuǎn)軸被劃分為4個(gè)軸段，轉(zhuǎn)子系統(tǒng)劃分為5個(gè)節(jié)點(diǎn)。轉(zhuǎn)軸的彈性模量為204GPa，泊松比為0.3。

表2 柔性轉(zhuǎn)子系統(tǒng)結(jié)構(gòu)參數(shù)

對圖4所示的柔性轉(zhuǎn)子系統(tǒng)進(jìn)行計(jì)算，判斷圓盤節(jié)點(diǎn)處的混沌行為，并分析軸承徑向間隙、圓盤不平衡力、轉(zhuǎn)軸剛度比等參數(shù)對系統(tǒng)混沌行為的影響規(guī)律。

3.1 徑向間隙影響分析

假設(shè)系統(tǒng)不平衡力為0，則軸承徑向間隙對系統(tǒng)混沌特性的影響如圖5所示。分別取40μm、60μm和80μm的徑向間隙計(jì)算系統(tǒng)的最大Lyapunov指數(shù)。

圖5 不同節(jié)點(diǎn)處的最大Lyapunov指數(shù)

系統(tǒng)最大Lyapunov指數(shù)小于0，表明系統(tǒng)運(yùn)動(dòng)是穩(wěn)定的；當(dāng)最大Lyapunov指數(shù)等于0時(shí)，系統(tǒng)為倍周期或擬周期分叉運(yùn)動(dòng)；當(dāng)最大Lyapunov指數(shù)大于0時(shí)，系統(tǒng)為混沌運(yùn)動(dòng)。

為了驗(yàn)證計(jì)算的準(zhǔn)確性，使用Runge－Kutta法進(jìn)行非線性動(dòng)力學(xué)方程組求解，得到轉(zhuǎn)子系統(tǒng)響應(yīng)的頻譜圖和Poincaré截面。從圖5可知，當(dāng)徑向間隙為80μm、轉(zhuǎn)速為4200r／min時(shí)，系統(tǒng)的最大Lyapunov指數(shù)等于0（圖中fVC為滾子軸承的變剛度振動(dòng)頻率）。圖6所示為圖5工況下圓盤節(jié)點(diǎn)處的頻譜圖和Poincaré截面，圖形表明系統(tǒng)為擬周期振動(dòng)。

圖6 轉(zhuǎn)速為4200r／min、徑向間隙為80μm時(shí)的響應(yīng)

圖5表明，徑向間隙為80μm、轉(zhuǎn)速為8000 r／min時(shí)，系統(tǒng)的最大Lyapunov指數(shù)大于0，圖7所示為對應(yīng)該工況的圓盤節(jié)點(diǎn)處的頻譜圖和Poincaré截面，其表明系統(tǒng)為混沌運(yùn)動(dòng)。

圖5還表明，當(dāng)徑向間隙為40μm時(shí)，系統(tǒng)的最大Lyapunov指數(shù)小于0，未出現(xiàn)混沌行為；當(dāng)徑向間隙為60μm 時(shí)，在0～2500r／min、6600～7200r／min轉(zhuǎn)速下的最大Lyapunov指數(shù)大于0，系統(tǒng)出現(xiàn)混沌行為；當(dāng)徑向間隙為80μm時(shí)，在0～3100r／min、5000～5400r／min、7500～8500 r／min轉(zhuǎn)速下出現(xiàn)混沌行為。

3.2 不平衡力影響分析

圖8所示為在其他條件不變前提下的不平衡力對系統(tǒng)混沌特性的影響。分別取圓盤的不平衡質(zhì)徑積為5×10－6kg·m、1×10－5kg·m、2×10－5kg·m，計(jì)算系統(tǒng)的最大Lyapunov指數(shù)。

從圖8可知，當(dāng)不平衡質(zhì)徑積為5×10－6kg·m時(shí)，在0～2000r／min，系統(tǒng)出現(xiàn)混沌行為；當(dāng)不平衡質(zhì)徑積為1×10－5kg·m時(shí)，在0～3700r／min，系統(tǒng)出現(xiàn)混沌行為；當(dāng)不平衡質(zhì)徑積為2×10－5kg·m時(shí)，在0～4400r／min，系統(tǒng)出現(xiàn)混沌行為。可以看出，不平衡力的存在會(huì)改變系統(tǒng)的混沌行為，系統(tǒng)出現(xiàn)混沌的轉(zhuǎn)速及范圍隨著不平衡力的增大逐漸增大。

圖7 轉(zhuǎn)速為8000r／min、徑向間隙為80μm時(shí)的響應(yīng)

圖8 不同不平衡力時(shí)的系統(tǒng)最大Lyapunov指數(shù)

3.3 剛度比影響分析

定義柔性轉(zhuǎn)子系統(tǒng)的剛度比為

式中，Kshaft為柔性轉(zhuǎn)軸的剛度；Kbearing為支承軸承剛度。

在其他參數(shù)不變的情況下，分別取長度為500mm、400mm、300mm的轉(zhuǎn)軸進(jìn)行剛度比計(jì)算，其對應(yīng)所得的剛度比分別為7.38×10－3、1.43×10－2、3.31×10－2，據(jù)此，可得出剛度對于柔性轉(zhuǎn)子系統(tǒng)動(dòng)力學(xué)特性的影響。

圖9所示為不同剛度比時(shí)圓盤節(jié)點(diǎn)處的振幅隨轉(zhuǎn)速變化曲線，可以看出，隨著剛度比的增大，振幅峰值對應(yīng)的轉(zhuǎn)速逐漸增大，即臨界轉(zhuǎn)速增大；振幅峰值也隨著剛度比增大而增大，且臨界轉(zhuǎn)速附近的峰值逐漸變多，可以看出隨著轉(zhuǎn)軸剛度增大，軸承非線性振動(dòng)對轉(zhuǎn)子系統(tǒng)的影響增大，導(dǎo)致系統(tǒng)響應(yīng)復(fù)雜。

圖9 不同剛度比時(shí)的圓盤節(jié)點(diǎn)振幅

圖10所示為不同剛度比時(shí)最大Lyapunov指數(shù)隨轉(zhuǎn)速的變化曲線。當(dāng)剛度比為7.38×10－3時(shí)，在0～800r／min出現(xiàn)混沌行為。當(dāng)剛度比為1.43×10－2時(shí)，在0～1200r／min、5200～5400r／min出現(xiàn)混沌行為。當(dāng)剛度比為3.31×10－2時(shí)，在0～3200r／min、4800～5900r／min出現(xiàn)混沌行為。

圖10 不同轉(zhuǎn)軸剛度比時(shí)的系統(tǒng)最大Lyapunov指數(shù)

可見，隨著轉(zhuǎn)軸剛度增大即剛度比增大，混沌運(yùn)動(dòng)的區(qū)間發(fā)生改變。隨著剛度比的增大，系統(tǒng)的混沌區(qū)間增大，軸承引起的轉(zhuǎn)子系統(tǒng)非線性行為明顯。對比柔性轉(zhuǎn)子和剛性轉(zhuǎn)子的混沌特性，發(fā)現(xiàn)軸承的非線性接觸力對剛性轉(zhuǎn)子系統(tǒng)的混沌特性影響大于柔性轉(zhuǎn)子。

4 實(shí)驗(yàn)

為驗(yàn)證本文方法對柔性轉(zhuǎn)子系統(tǒng)動(dòng)力學(xué)行為預(yù)測的準(zhǔn)確性，使用圖11所示的高速滾子軸承柔性轉(zhuǎn)子實(shí)驗(yàn)器進(jìn)行測試。動(dòng)力裝置的高速電主軸，最高轉(zhuǎn)速可達(dá)24 000r／min，實(shí)驗(yàn)滾子軸承參數(shù)如表1所示，實(shí)驗(yàn)過程中滾子軸承承受的徑向載荷為2000N。采用非接觸式的電渦流傳感器測量滾子軸承及圓盤處的徑向振動(dòng)位移。

圖12所示為轉(zhuǎn)子圓盤處振動(dòng)位移幅值計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果對比，可以看出，計(jì)算得到的一階臨界轉(zhuǎn)速為4900r／min，實(shí)驗(yàn)得到的一階臨界轉(zhuǎn)速為4700r／min，與計(jì)算結(jié)果較為接近；實(shí)驗(yàn)測試振幅與計(jì)算結(jié)果也較為接近，證明了本文方法的準(zhǔn)確性。

圖11 滾動(dòng)軸承轉(zhuǎn)子系統(tǒng)實(shí)驗(yàn)器

圖12 幅值計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果對比

5 結(jié)論

（1）軸承徑向間隙是影響轉(zhuǎn)子系統(tǒng)非線性振動(dòng)特性的重要參數(shù)。隨著軸承徑向間隙的增大，系統(tǒng)的混沌區(qū)間逐漸增大、變多。

（2）不平衡力的存在對系統(tǒng)混沌行為也有較大的影響。系統(tǒng)出現(xiàn)混沌的轉(zhuǎn)速及范圍隨著不平衡力的增大逐漸增大。

（3）隨著轉(zhuǎn)軸剛度比增大，轉(zhuǎn)子振幅峰值以及對應(yīng)的轉(zhuǎn)速均逐漸增大。共振轉(zhuǎn)速附近產(chǎn)生混沌運(yùn)動(dòng)對應(yīng)的轉(zhuǎn)速增大時(shí)，系統(tǒng)混沌區(qū)間也增多，軸承非線性對系統(tǒng)響應(yīng)影響明顯。

［1］Jang G H，Jeong S W.Nonlinear Excitation Model of Ball Bearing Waviness in a Rigid Rotor Supported by Two or More Ball Bearings Considering Five Degrees of Freedom［J］.Transactions of the ASME，2002，124：82－90.

［2］Harsha S P.Non－linear Dynamic Analysis of an Unbalanced Rotor Supported by Roller Bearing［J］.Chaos，Solitons and Fractals，2005，26：47－66.

［3］Harsha S P.Non－linear Dynamic Response of a Balanced Rotor Supported by Rolling Element Bearings Due to Radial Internal Clearance Effect［J］.Mechanism and Machine Theory，2006，41：688－706.

［4］Bai Changqing，Xu Qingyu.Dynamic Model of Ball Bearings with Internal Clearance and Waviness［J］.Journal of Sound and Vibration，2006，294：23－48.

［5］高尚晗，龍新華，孟光.主軸－滾動(dòng)軸承系統(tǒng)三種分岔形式［J］.振動(dòng)與沖擊，2009，28（4）：59－64.Gao Shanghan，Long Xinhua，Meng Guang.Three Bifurcation Types in the Spindle－ball Bearing System［J］.Journal of Vibration and Shock，2009，28（4）：59－64.

［6］崔立，劉長利，鄭建榮.圓柱滾子軸承剛性轉(zhuǎn)子系統(tǒng)周期運(yùn)動(dòng)分岔研究［J］.振動(dòng)、測試與診斷，2011，31（5）：637－641.Cui Li，Liu Changli，Zheng Jianrong.Periodic Motion Bifurcation of Rigid Rotor System in Cylinder Roller Bearing［J］.Journal of Vibration，Measurement ＆ Diagnosis，2011，31（5）：637－641.

［7］Laha S K，Kakoty S K.Non－linear Dynamic Analysis of a Flexible Rotor Supported on Porous Oil Journal Bearings［J］.Communications in Nonlinear Science and Numerical Simulation，2011，16：1617－1631.

［8］Villa C，Sinou J J，Thouverez F.Stability and Vibration Analysis of a Complex Flexible Rotor Bearing System［J］.Communications in Non－linear Science and Numerical Simulation，2008，13：804－821.

［9］Sinou J J.Non－linear Dynamics and Contacts of an Unbalanced Flexible Rotor Supported on Ball Bearings［J］.Mechanism and Machine Theory，2009，44：1713－1732.

［10］Wang Maocheng.Finite Element Method［M］.Beijing：Tsinghua University Press，2008.

［11］Harris T A.Rolling Bearing Analysis［M］.New York，USA：John Wiley ＆Sons，2001.

［12］Choi S K，Noah S T.Response and Stability Analysis of Piecewise Linear Oscillations under Multi－forcing Frequencies［J］.Nonlinear Dynamics，1992，3：105－121

［13］Wolf A，Swift J B，Swinney H L.Determining Lyapunov Exponents from a Time Series［J］.Physica D.，1985，16：285－317.