雙參變量下單級齒輪傳動系統(tǒng)分岔/沖擊特性分析

田亞平， 褚衍東， 饒曉波

(1.蘭州交通大學 機電工程學院，蘭州 730070； 2. 蘭州交通大學 甘肅省軌道交通裝備系統(tǒng)動力學與可靠性重點實驗室，蘭州 730070)

因齒側(cè)間隙、時變嚙合剛度、加工制造誤差、質(zhì)量偏心等因素導(dǎo)致齒輪系統(tǒng)表現(xiàn)出強非線性的動力學振動特征。Kahraman等[1]考慮齒側(cè)間隙和誤差激勵因素建立了單自由度單級齒輪系統(tǒng)動力學模型，發(fā)現(xiàn)了次諧和混沌響應(yīng)。Vinayak等[2]研究了齒輪嚙合綜合傳遞誤差對系統(tǒng)非線性動力學行為的影響。劉曉寧等[3]采用增量諧波平衡法對三自由度齒輪系統(tǒng)周期運動進行了分析并判穩(wěn)。茍向鋒等[4]采用相圖、Lyapnuov指數(shù)譜、Poincaré截面法找出了三自由度單級齒輪副隨參數(shù)變化通向混沌運動的途徑。王曉筍等[5]研究了含磨損故障的單級齒輪傳動系統(tǒng)的分岔行為及通向混沌的途徑。

結(jié)合Floquet穩(wěn)定性理論的PNF法(Poincaré-Newton-Floquet)是周期運動求解及判穩(wěn)同步進行的一種高效數(shù)值方法。Luo等[6]運用PNF法研究了含裂紋以及碰磨耦合故障的轉(zhuǎn)子系統(tǒng)的周期運動的穩(wěn)定性問題，Han[7-8]運用PNF法研究了碰磨轉(zhuǎn)子系統(tǒng)的周期運動穩(wěn)定性問題，李同杰等[9]運用PNF法研究了行星齒輪傳動系統(tǒng)周期運動的穩(wěn)定性問題。上述研究中均采用單參變量分析系統(tǒng)的分岔特性，無法全面了解參數(shù)間耦合下系統(tǒng)的非線性特性轉(zhuǎn)遷規(guī)律。雙參平面內(nèi)研究系統(tǒng)分岔/沖擊特性能更清晰的了解系統(tǒng)參數(shù)之間的耦合關(guān)系并確定系統(tǒng)穩(wěn)定運行的參數(shù)范圍，為工程實際中的齒輪結(jié)構(gòu)的參數(shù)優(yōu)化提供理論指導(dǎo)。Gou等[10-11]采用胞映射和逃逸算法研究了雙參平面內(nèi)的純扭轉(zhuǎn)齒輪副的分岔特性。但采用PNF法和延續(xù)算法研究齒輪系統(tǒng)雙參平面分岔特性的文獻還鮮有報道。

本文以三自由度單級齒輪傳動系統(tǒng)的動力學模型為研究對象，采用PNF法和延續(xù)追蹤法求解系統(tǒng)周期運動并判別分岔類型，然后采用雙參平面?zhèn)尾蕡D探尋系統(tǒng)的運動分岔/沖擊特性隨參數(shù)轉(zhuǎn)遷和參數(shù)間的耦合關(guān)系。

1 單級齒輪非線性動力學模型

單級齒輪傳動系統(tǒng)動力學模型如圖1所示。圖中，mgi，Igi，rgi，θgi(i=1,2)分別代表主從動齒輪的質(zhì)量、轉(zhuǎn)動慣量、基圓半徑、扭轉(zhuǎn)角位移；bi，cbi，kbi，F(xiàn)bi，ygi分別表示主從動齒輪支承軸承的間隙、阻尼、支承剛度和作用于軸承的徑向預(yù)載荷和徑向位移。Tg1，Tg2為輸入輸出扭矩均值。齒輪嚙合的時變剛度、齒側(cè)間隙、嚙合線性阻尼、齒輪嚙合誤差分別用kh(t)，bh，ch，e(t)表示。

圖1 單級齒輪傳動系統(tǒng)非線性模型

由文獻[3]知其系統(tǒng)的量綱一化動力學方程為

(1)

(2)

式中:bhi為軸承支承徑向間隙(i=1,2)和半齒側(cè)間隙(i=3)；bi為其對應(yīng)的量綱一化間隙。

y3為量綱一化齒輪嚙合的動態(tài)傳動誤差與靜態(tài)傳動誤差的差值。

y3(τ)=(rg1θg1(τ)-rg2θg2(τ)+yg1-

yg2-e(τ))/bc

(3)

(4)

2 系統(tǒng)周期運動分岔轉(zhuǎn)遷規(guī)律的求解過程

雙參平面內(nèi)，系統(tǒng)(4)的分岔/沖擊特性及其轉(zhuǎn)遷規(guī)律采用PNF法周期不動點的求解和周期數(shù)的統(tǒng)計、周期解延續(xù)算法[12]、分岔類型的判斷(Floquet乘子法[13])、齒輪嚙合沖擊類型的判斷來實現(xiàn)。其具體求解過程如下。

2.1 PNF法周期數(shù)統(tǒng)計及追蹤

PNF法是一種Floquet理論和打靶法相結(jié)合的同時進行周期運動求解并判穩(wěn)的數(shù)值法，在周期運動求解中統(tǒng)計出其周期數(shù)。將系統(tǒng)方程式(4)改寫為周期激勵的n+1維狀態(tài)方程

Rn×R,n=2N,f(X,τ)=f(X,τ+T0),τ>0

(5)

式中：T0=2π/Ω為系統(tǒng)激振力周期；N為系統(tǒng)的自由度，本文N=3。系統(tǒng)的m倍周期運動定義為

X(τ)=X(τ+T), ?τ>0

(6)

式中：T=mT0，m為正整數(shù)。

在n+1維狀態(tài)空間中定義n維Poincaré截面

∑={(X(τ),τ)|mod(τ,T)=0}

(7)

在該截面上定義Poincaré映射

(8)

連續(xù)系統(tǒng)(5)的周期解等價于求解如下代數(shù)方程的不動點XF

Q(X)=P(X)-X=0

(9)

式(9)的求解由打靶法思想進行：事先給出XF的一個近似解Xk∈Σ，然后按照Newton-Raphson方法構(gòu)建迭代公式，以確定更加接近XF的向量Xk+1。

Xk+1=[I-DP(Xk)]-1[P(Xk)-

DP(Xk)·Xk]

(10)

式中：I為單位方陣；P(Xk)為Xk處的Poincaré映射；DP(Xk)為Poincaré映射在Xk處的Jcaobi矩陣(離散轉(zhuǎn)移矩陣)。實際計算P(Xk)的初值由下式迭代

(11)

式(11)在T時刻的解便是Xk的Poincaré映射，即P(Xk)=X(T)。DP(Xk)的n×n個初值求解

(12)

式(12)在T時刻的解就是Poincaré映射在Xk處的Jacobi矩陣，即Φ(T)=DP(Xk)，?fX(Xk,τ)/?Xk為狀態(tài)方程(5)在Xk處對應(yīng)的Jacobi矩陣。聯(lián)立式(11)和式(12)得

(13)

式(13)以(Xk,I)為初值積分一個周期T獲得系統(tǒng)的Xk+1和DP(Xk)。令Xk=Xk+1作為下次迭代的初始值，通過式(13)反復(fù)迭代，直到滿足精度要求的不動點XF，就是滿足等式(6)的m倍周期運動X(τ)的不動點。

找到穩(wěn)定周期運動后，根據(jù)延續(xù)分岔算法引入分岔參數(shù)λ，采用延續(xù)追蹤的方式進行周期追蹤。構(gòu)建代數(shù)方程

G(λ,X)=P(λ,X)-X=0

(14)

分岔參數(shù)λ的延續(xù)追蹤過程就是式(14)關(guān)于λ的解曲線求解過程。式(14)的解曲線問題轉(zhuǎn)化為求解常微分方程

(15)

的問題。設(shè)已求得λ=λk(k=0,1,2,…)時，周期T的Poincaré截面的不動點Xk，采用Euler型數(shù)值積分方法，可得λ=λk+1對應(yīng)不動點Xk+1的初始不動點Xk+1,0的預(yù)測公式

(16)

式中：Δλ為步長，GX(λk,Xk)=I-DP(λk，Xk),Gλ(λk,Xk)=Pλ(λk,Xk)。DP(λk，Xk)由式(12)求得，Pλ(λk,Xk)由打靶法思想建立的如下微分方程組求得

(17)

以(Xk，λk)為初值積分一個Poincaré映射周期T獲(Xk，Pλ(λk,Xk))。

從分岔參數(shù)λk開始，通過式(13)迭代獲得周期不動點后，再以式(16)進行延續(xù)追蹤得λk+Δλ處的初始迭代點，重復(fù)式(13)和(16)的迭代過程，可求得一系列外參數(shù)λ值的周期不動點Xk及其離散轉(zhuǎn)移矩陣DP(Xk)，完成周期m的追蹤。周期解的穩(wěn)定性由離散轉(zhuǎn)移矩陣DP(Xk)的特征值(Floquet乘子模)的最大值|λ|max確定。當判斷周期解失穩(wěn)時判斷分岔類型并重新尋求新的周期數(shù)。

2.2 周期和沖擊類型的表示方法

定義符號“I/P”為系統(tǒng)的齒輪副嚙合沖擊類型和周期數(shù)，I表示齒輪嚙合的沖擊類型：I=0,1,2分別表示無沖擊、單邊沖擊、雙邊沖擊類型，P表示周期數(shù)。設(shè)量綱一化齒側(cè)間隙為2b，依據(jù)文獻[14]，當xmin>b時齒面接觸無沖擊；當|xmin|

3 系統(tǒng)雙參分岔仿真分析

為了分析系統(tǒng)(4)的動力學分岔、沖擊、幅值跳躍等特性，取仿真參數(shù)：ζ11=ζ22=0.1，ζ13=ζ23=0.012 5，ζ33=0.05，k11=k22=1.25，fm=0.2，fah1=0.05，fb1=fb2=0.1，b1=b2=0，b3=1，仿真初值取為0。在系統(tǒng)特定參數(shù)下，雙參平面內(nèi)用不同的顏色代表不同的“I/P”類型，其區(qū)域的邊界線即為分岔曲線。Ω-α雙參平面分岔圖如圖2所示。齒輪系統(tǒng)常穩(wěn)定運行在短周期狀態(tài)下，為簡化分岔圖對P≥65的長周期、擬周期(Qusi-P)和混沌(chaos)運動均用“I/N”表示(白色區(qū)域)。當系統(tǒng)從狀態(tài)“0/P”向“1/P”狀態(tài)轉(zhuǎn)遷中，系統(tǒng)從無沖擊區(qū)域在x=b截面處穿越到脫嚙區(qū)，在該點發(fā)生了擦切(G Bif)，其交線為擦切分岔曲線。當系統(tǒng)從“I/P”向“I/2P”轉(zhuǎn)遷時發(fā)生了倍化分岔(PD Bif)，其交線為倍化分岔曲線。當系統(tǒng)從“I/1”向“I/n”轉(zhuǎn)遷時系統(tǒng)從單周期經(jīng)hopf分岔(PF Bif)為擬周期運動或經(jīng)激變分岔(CIC)為混沌運動，具體由Floquet乘子判斷。

齒輪系統(tǒng)的重合度和轉(zhuǎn)速是影響系統(tǒng)運轉(zhuǎn)穩(wěn)定性的主要因素，而時變嚙合剛度幅值系數(shù)α是表征重合度大小的主要指標，量綱一頻率Ω是表征轉(zhuǎn)速大小的主要指標。因此，選雙參(Ω，α)進行系統(tǒng)分岔特性轉(zhuǎn)遷規(guī)律分析，如圖2所示。橫坐標量綱一化頻率Ω∈[1,2]，縱坐標嚙合剛度幅值系數(shù)α∈[0,0.49]。在雙參平面內(nèi)，系統(tǒng)出現(xiàn)了倍化(PD Bif)、擦切(G Bif)、幅值跳躍(JP Bif)、激變、hopf等分岔現(xiàn)象，周期1、周期2、周期4、周期8、周期16、擬周期、混沌等運動狀態(tài)和三種沖擊狀態(tài)共存。

圖2 Ω-α雙參平面分岔圖

沿Ω方向系統(tǒng)分岔轉(zhuǎn)遷規(guī)律較復(fù)雜。當α∈[0,0.05]范圍內(nèi)(小幅值剛度波動大重合度狀態(tài)下)系統(tǒng)表現(xiàn)出擦切和倍化分岔，系統(tǒng)在“0/1”、“1/1”和“1/2”運動狀態(tài)轉(zhuǎn)遷出現(xiàn)周期泡現(xiàn)象，并在Ω∈[1.4,1.6]區(qū)域內(nèi)系統(tǒng)處于單邊沖擊脫嚙狀態(tài)。當α∈(0.05,0.1]范圍內(nèi)系統(tǒng)倍化分岔出現(xiàn)了“1/4”、“1/8”運動狀態(tài)。當α∈(0.1,0.4]范圍內(nèi)系統(tǒng)通過hopf和鞍結(jié)激變分岔，單周期運動進入混沌運動后經(jīng)逆倍化分岔經(jīng)周期“1/16”、“1/8”、…、到穩(wěn)定“0/1”運動狀態(tài)，單邊沖擊脫嚙的范圍逐漸擴大。當α∈(0.4,0.49]范圍內(nèi)系統(tǒng)經(jīng)歷了單周期、激變分岔進入雙邊沖擊齒背接觸混沌“2/N”狀態(tài)、混沌吸引子激變演變?yōu)閱芜厸_擊脫嚙混沌“1/N”運動、經(jīng)短暫的周期運動后Hopf分岔進入擬周期運動、逆Hopf分岔進入單周期運動狀態(tài)。

(a) Ω-α-x5

(b) α-Ω-x5

沿α方向系統(tǒng)分岔運動的轉(zhuǎn)遷規(guī)律較為簡單。當Ω∈[1,1.35]∪[1.58,2]范圍內(nèi)系統(tǒng)處于單周期無沖擊穩(wěn)定運動狀態(tài)，僅在Ω∈[1.6,1.9]區(qū)域內(nèi)出現(xiàn)了擦切分岔形成了單邊沖擊周期運動狀態(tài)。在Ω∈(1.35,1.58)的共振頻率附近，系統(tǒng)通過倍周期分岔形成周期2、周期4、周期8、…、混沌運動的演化過程。特別在α∈(0.22,0.3)范圍內(nèi)發(fā)生小窗口的多次幅值跳躍雙邊沖擊周期2運動狀態(tài)，即在該區(qū)域內(nèi)系統(tǒng)的系統(tǒng)動力特性非常復(fù)雜，設(shè)計中應(yīng)盡量避開該區(qū)域選取參數(shù)。

為了驗證PNF法和延續(xù)追蹤法求解的雙參平面分岔圖(圖2)的正確性，沿圖2中的α和Ω方向等間距選取截面采用4～5階Runge-kutta數(shù)值法仿真獲得系統(tǒng)(4)的Ω-α-x5三維分岔圖如圖3所示。圖中，齒輪嚙合相對綜合誤差x5隨Ω和α的變化，系統(tǒng)的動力學分岔特性的轉(zhuǎn)遷規(guī)律及分岔點位置與圖2相一致，說明PNF法和延續(xù)追蹤法求解雙參平面分岔圖的方法是有效的。

取圖2中Ω=1.45截面的α-x5分岔圖如圖4(a)所示。系統(tǒng)通過倍化分岔進入混沌運動，其間出現(xiàn)了幅值跳躍分岔現(xiàn)象。取圖2中α=0.45截面的Ω-x5分岔圖如圖4(b)所示。系統(tǒng)由混沌運動激變退化為周期8運動、再經(jīng)過鞍結(jié)分岔幅值跳躍為周期2運動，經(jīng)hopf分岔進入擬周期運動，逆hopf分岔退化為周期2運動，激變分岔為混沌運動、逆倍化分岔退化為周期1運動的轉(zhuǎn)遷規(guī)律。

取圖4(b)中Ω=1.5截面關(guān)于(x5,x6)Poincaré映射圖如圖5(a)所示。圖5(a)中出現(xiàn)近似周期2的黑點，將其放大可看見兩個不穩(wěn)定的焦點向外擴散，當Ω增大到1.505時系統(tǒng)焦點分岔為兩個封閉的hopf圈(圖5(b)所示)，即在Ω=1.5處發(fā)生了hopf分岔，考查其Floquet乘子為[0.791 1, -0.889 2±0.459 5i, -0.624 2±0.553 1i, -0.310 4], |λ|max=|-0.889 2±0.459 5i|=1.049，即從一對共軛復(fù)數(shù)穿出單位圓。在圖2中取α=0.079，Ω=1.632、1.638、1.642值的(x5,x6)相圖如圖5(c)、圖5(d)、圖5(e)所示，相圖展示了系統(tǒng)從擦切前的單邊沖擊運動向擦切后的無沖擊運動的轉(zhuǎn)遷過程，很明顯在x5=1截面處發(fā)生了擦切運動。

采用Runge-kutta數(shù)值法仿真的三維/二維分岔圖、Poincaré映射圖和相圖驗證了綜合PNF法和延續(xù)追蹤法的雙參平面分岔求解的正確性。因延續(xù)追蹤法省去了大量的瞬態(tài)過程求解時間，故是一種有效的雙參平面分岔計算方法。

(a) Ω=1.45

(b) α=0.45

(a) Ω=1.50, α=0.45

(b) Ω=1.505, α=0.45

(c) Ω=1.634, α=0.079

(d) Ω=1.638, α=0.079

(e) Ω=1.642, α=0.079

雙參平面分岔圖揭示了參數(shù)Ω、α范圍內(nèi)系統(tǒng)分岔特性的轉(zhuǎn)遷規(guī)律，獲得了周期運動、擬周期運動、混沌運動的分岔曲線。工程實際應(yīng)用中，為了延長齒輪設(shè)備的運行壽命降低事故的發(fā)生的幾率和設(shè)備維修費用，常常希望齒輪處于“0/1”運動狀態(tài)，因此在齒輪裝備設(shè)計中合理選擇齒輪結(jié)構(gòu)參數(shù)、潤滑油黏度及制造精度使系統(tǒng)量綱一化參數(shù)在雙參平面 “0/1”區(qū)域內(nèi)為佳。

4 結(jié) 論

(1) 在雙參分岔平面內(nèi)，提出了一種PNF法和延續(xù)算法相結(jié)合的非線性動力系統(tǒng)周期運動求解和嚙合沖擊類型判別(I/P分岔)的新方法，并用Runge-Kutta數(shù)值法的相圖、三維/二維分岔圖等驗證了其有效性。

(2) 以三自由度單級齒輪傳動系統(tǒng)為分析模型，在時變嚙合剛度幅值系數(shù)和量綱一化頻率雙參平面內(nèi)分析了系統(tǒng)分岔/沖擊特性及轉(zhuǎn)遷規(guī)律，劃分了雙參耦合下系統(tǒng)穩(wěn)定周期運動和失穩(wěn)混沌運動參數(shù)區(qū)域。

(3) 在時變剛度幅值系數(shù)和量綱一化頻率特定參數(shù)耦合狀態(tài)下，避開共振頻率范圍和小時變嚙合剛度幅值系數(shù)(大重合度)范圍內(nèi)選取參數(shù)能提高齒輪傳動的穩(wěn)定性降低疲勞磨損。