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

田亞平， 褚衍東， 饒曉波

(1.蘭州交通大學(xué) 機(jī)電工程學(xué)院，蘭州 730070； 2. 蘭州交通大學(xué) 甘肅省軌道交通裝備系統(tǒng)動(dòng)力學(xué)與可靠性重點(diǎn)實(shí)驗(yàn)室，蘭州 730070)

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

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

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

1 單級(jí)齒輪非線性動(dòng)力學(xué)模型

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

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

由文獻(xiàn)[3]知其系統(tǒng)的量綱一化動(dòng)力學(xué)方程為

(1)

(2)

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

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

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

yg2-e(τ))/bc

(3)

(4)

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

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

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

PNF法是一種Floquet理論和打靶法相結(jié)合的同時(shí)進(jìn)行周期運(yùn)動(dòng)求解并判穩(wěn)的數(shù)值法，在周期運(yùn)動(dòng)求解中統(tǒng)計(jì)出其周期數(shù)。將系統(tǒng)方程式(4)改寫為周期激勵(lì)的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倍周期運(yùn)動(dòng)定義為

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)的周期解等價(jià)于求解如下代數(shù)方程的不動(dòng)點(diǎn)XF

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

(9)

式(9)的求解由打靶法思想進(jìn)行：事先給出XF的一個(gè)近似解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)移矩陣)。實(shí)際計(jì)算P(Xk)的初值由下式迭代

(11)

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

(12)

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

(13)

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

找到穩(wěn)定周期運(yùn)動(dòng)后，根據(jù)延續(xù)分岔算法引入分岔參數(shù)λ，采用延續(xù)追蹤的方式進(jìn)行周期追蹤。構(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,…)時(shí)，周期T的Poincaré截面的不動(dòng)點(diǎn)Xk，采用Euler型數(shù)值積分方法，可得λ=λk+1對(duì)應(yīng)不動(dòng)點(diǎn)Xk+1的初始不動(dòng)點(diǎn)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)為初值積分一個(gè)Poincaré映射周期T獲(Xk，Pλ(λk,Xk))。

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

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

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

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

為了分析系統(tǒng)(4)的動(dòng)力學(xué)分岔、沖擊、幅值跳躍等特性，取仿真參數(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)定運(yùn)行在短周期狀態(tài)下，為簡化分岔圖對(duì)P≥65的長周期、擬周期(Qusi-P)和混沌(chaos)運(yùn)動(dòng)均用“I/N”表示(白色區(qū)域)。當(dāng)系統(tǒng)從狀態(tài)“0/P”向“1/P”狀態(tài)轉(zhuǎn)遷中，系統(tǒng)從無沖擊區(qū)域在x=b截面處穿越到脫嚙區(qū)，在該點(diǎn)發(fā)生了擦切(G Bif)，其交線為擦切分岔曲線。當(dāng)系統(tǒng)從“I/P”向“I/2P”轉(zhuǎn)遷時(shí)發(fā)生了倍化分岔(PD Bif)，其交線為倍化分岔曲線。當(dāng)系統(tǒng)從“I/1”向“I/n”轉(zhuǎn)遷時(shí)系統(tǒng)從單周期經(jīng)hopf分岔(PF Bif)為擬周期運(yùn)動(dòng)或經(jīng)激變分岔(CIC)為混沌運(yùn)動(dòng)，具體由Floquet乘子判斷。

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

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

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

(a) Ω-α-x5

(b) α-Ω-x5

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

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

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

取圖4(b)中Ω=1.5截面關(guān)于(x5,x6)Poincaré映射圖如圖5(a)所示。圖5(a)中出現(xiàn)近似周期2的黑點(diǎn)，將其放大可看見兩個(gè)不穩(wěn)定的焦點(diǎn)向外擴(kuò)散，當(dāng)Ω增大到1.505時(shí)系統(tǒng)焦點(diǎn)分岔為兩個(gè)封閉的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，即從一對(duì)共軛復(fù)數(shù)穿出單位圓。在圖2中取α=0.079，Ω=1.632、1.638、1.642值的(x5,x6)相圖如圖5(c)、圖5(d)、圖5(e)所示，相圖展示了系統(tǒng)從擦切前的單邊沖擊運(yùn)動(dòng)向擦切后的無沖擊運(yùn)動(dòng)的轉(zhuǎn)遷過程，很明顯在x5=1截面處發(fā)生了擦切運(yùn)動(dòng)。

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

(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ùn)動(dòng)、擬周期運(yùn)動(dòng)、混沌運(yùn)動(dòng)的分岔曲線。工程實(shí)際應(yīng)用中，為了延長齒輪設(shè)備的運(yùn)行壽命降低事故的發(fā)生的幾率和設(shè)備維修費(fèi)用，常常希望齒輪處于“0/1”運(yùn)動(dòng)狀態(tài)，因此在齒輪裝備設(shè)計(jì)中合理選擇齒輪結(jié)構(gòu)參數(shù)、潤滑油黏度及制造精度使系統(tǒng)量綱一化參數(shù)在雙參平面 “0/1”區(qū)域內(nèi)為佳。

4 結(jié) 論

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

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

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