王 棟，張艷龍，王 麗

（1.蘭州交通大學(xué) 機(jī)電工程學(xué)院，蘭州 730070； 2.蘭州城市學(xué)院 數(shù)學(xué)學(xué)院，蘭州 730070）

由于系統(tǒng)間隙的存在，動(dòng)力機(jī)械系統(tǒng)中碰撞振動(dòng)現(xiàn)象時(shí)常發(fā)生，如高速運(yùn)行列車的輪軌碰撞、齒輪傳動(dòng)系統(tǒng)的輪齒碰撞、車輛氣動(dòng)控制門關(guān)閉時(shí)和車輛的碰撞、運(yùn)行動(dòng)車上純凈水水桶和固定裝置的碰撞等。研究這類系統(tǒng)的動(dòng)力學(xué)行為對(duì)系統(tǒng)減振降噪等優(yōu)化設(shè)計(jì)具有重要意義。近年來，許多學(xué)者研究了碰撞振動(dòng)系統(tǒng)的非常規(guī)分岔和穩(wěn)定性。

Luo等研究了兩自由度含間隙碰撞振動(dòng)系統(tǒng)的碰撞運(yùn)動(dòng)和參數(shù)匹配[1]。伍新等研究了3自由度含間隙碰撞振動(dòng)系統(tǒng)的Neimark-Sacher分岔，并對(duì)其進(jìn)行了反控制[2]。馮進(jìn)鈐等基于圖胞映射理論，論述了Duffing碰撞振動(dòng)系統(tǒng)中擦邊誘導(dǎo)激變的全局動(dòng)力學(xué)行為[3]。系統(tǒng)穩(wěn)定性分析方面，Lyapunov指數(shù)譜是判定動(dòng)力系統(tǒng)運(yùn)動(dòng)穩(wěn)定性和混沌特性的重要工具。目前對(duì)光滑動(dòng)力系統(tǒng)Lyapunov指數(shù)譜的研究已經(jīng)比較成熟。由于非光滑動(dòng)力系統(tǒng)在非光滑點(diǎn)的Jacobi矩陣不存在，使得光滑動(dòng)力系統(tǒng)Lyapunov指數(shù)譜的計(jì)算方法在非光滑動(dòng)力系統(tǒng)中不再適用。Stefanski等采用同步現(xiàn)象的特性估算了具有不連續(xù)性的兩個(gè)相同離散動(dòng)力系統(tǒng)的最大Lyapunov指數(shù)[4-5]。Souza等通過在碰撞瞬時(shí)附加一定的轉(zhuǎn)換條件，將光滑離散動(dòng)力系統(tǒng)的Lyapunov指數(shù)譜計(jì)算方法運(yùn)用到兩個(gè)特殊碰撞振動(dòng)系統(tǒng)的穩(wěn)定性分析中[6]。金俐等利用Poincaré映射分析法得到了一類單自由度碰撞振動(dòng)系統(tǒng)Lyapunov指數(shù)譜的計(jì)算方法[7]。Li等利用Poincaré映射分析法得到了一類2自由度雙側(cè)約束振動(dòng)系統(tǒng)Lyapunov指數(shù)譜的計(jì)算方法[8]。魏艷輝等通過構(gòu)造的Poincaré映射及其Jacobi矩陣，得到了一類2自由度單碰振動(dòng)系統(tǒng)的Lyapunov指數(shù)譜的計(jì)算方法[9]。周期泡現(xiàn)象分析方面，一些學(xué)者也開展了部分研究。石建飛等研究了雙參數(shù)平面下的Duffing系統(tǒng)的Lyapunov指數(shù)和分岔特性[10]。Peng研究了遲滯人口模型中的多重分岔和周期泡現(xiàn)象，并分析了系統(tǒng)的穩(wěn)定性[11]。Yang等分析了需求與供給相互作用的離散經(jīng)濟(jì)模型中的周期泡現(xiàn)象和系統(tǒng)穩(wěn)定性[12]。Gou等研究了雙參平面上直齒輪副運(yùn)動(dòng)時(shí)出現(xiàn)的周期泡現(xiàn)象[13]。對(duì)于3自由度碰撞振動(dòng)系統(tǒng)，由于系統(tǒng)維數(shù)的增加，系統(tǒng)的動(dòng)力學(xué)行為更加復(fù)雜豐富，Lyapunov指數(shù)譜的計(jì)算也更加繁瑣?，F(xiàn)有文獻(xiàn)中對(duì)3自由度單碰振動(dòng)系統(tǒng)中周期泡現(xiàn)象的發(fā)生以及Lyapunov指數(shù)譜的計(jì)算鮮有報(bào)道。

本文基于車輛氣動(dòng)控制門關(guān)閉時(shí)和車輛的碰撞、運(yùn)行動(dòng)車上純凈水水桶和固定裝置的碰撞、運(yùn)行的機(jī)床工作臺(tái)和行程限位機(jī)構(gòu)的碰撞等工程振動(dòng)問題抽象簡化得到一類3自由度單碰振動(dòng)系統(tǒng)的力學(xué)模型，運(yùn)用Poincaré映射分析法給出了該系統(tǒng)的Lyapunov指數(shù)譜的計(jì)算方法，分析了一定參數(shù)變化范圍下系統(tǒng)的穩(wěn)定性，并通過數(shù)值仿真，結(jié)合系統(tǒng)的分岔圖和相圖，發(fā)現(xiàn)了一定參數(shù)下系統(tǒng)存在周期泡和混沌泡現(xiàn)象。本文的研究方法在理論上可以推廣到更高自由度，但實(shí)際計(jì)算會(huì)很復(fù)雜。

1 力學(xué)模型

如圖1所示。阻尼系數(shù)為C1、C2和C3的線性阻尼器和剛度為K1、K2和K3的線性彈簧分別與質(zhì)量為M1、M2和M3的質(zhì)塊連接，簡諧力Pisin(ΩΤ+τ)(i=1,2,3)分別作用于3個(gè)質(zhì)塊，質(zhì)塊只作水平方向的運(yùn)動(dòng)。當(dāng)質(zhì)塊M1的位移X1(T)等于間隙B時(shí)質(zhì)塊與剛性約束A發(fā)生碰撞，碰撞過程由速度恢復(fù)系數(shù)R確定，并且不計(jì)碰撞持續(xù)時(shí)間。

圖1 3自由度單碰振動(dòng)系統(tǒng)動(dòng)力學(xué)模型

質(zhì)塊M1的沖擊方程為

其中：˙1-和˙1+分別為質(zhì)塊M1碰撞前后的瞬時(shí)速度。

令Ψ表示方程式(1)的正則模態(tài)矩陣，ω1、ω2和ω3分別表示在質(zhì)塊位移(x1＜b)下系統(tǒng)的固有頻率。作坐標(biāo)變換Χ=Ψξ。方程式(1)可解耦為

式中：I是一個(gè)3×3階單位矩陣；C和Λ是3×3階對(duì)角矩陣，且

通過模態(tài)疊加法可以得到方程式(1)的解為

式中：Ψij是正則模態(tài)矩陣Ψ的元素；；aj和bj是積分常數(shù)，Aj和Bj為振幅常數(shù)。令則式(1)可寫成如下 1 階狀態(tài)方程

2 Jacobi矩陣及Lyapunov指數(shù)譜

圖2 軌線在相空間穿越Poincaré截面情況

由此可得這幾類映射及其Jacobi矩陣：

此映射的切換面為Σ-，方程為h(zc)=x1-b=0，由式(8)可得到局部映射Pa的Jacobi矩陣表達(dá)式

(2)映射Pb：Σ-→Σ+

根據(jù)質(zhì)塊1的碰撞關(guān)系式(2)，可得到局部映射Pb的Jacobi矩陣表達(dá)式

此映射的切換面為定相位面Φφc，切換面方程為φ=φc，由式(8)可得到Pc的Jacobi矩陣表達(dá)式

由于該映射所描述的軌線是光滑的，通過數(shù)值方法可求得其數(shù)值解。其Jacobi矩陣可由多元函數(shù)求導(dǎo)法則得到

由式(7)定義Poincaré映射P可由以上四種映射復(fù)合表示為

根據(jù)復(fù)合映射求導(dǎo)的鏈?zhǔn)椒▌t，Poincaré映射P的Jacobi矩陣可寫成

3 數(shù)值模擬

事先通過選取多組系統(tǒng)參數(shù)對(duì)系統(tǒng)的動(dòng)力學(xué)行為進(jìn)行了大量的數(shù)值仿真，最終取系統(tǒng)參數(shù)m1=1.0，m2=0.4，m3=1.64，k1=1.0，k2=2.0，k3=3.0，γ=0.1，f1=1.0，f2=0，f3=0，R=0.8，b=0.25，在該系統(tǒng)中發(fā)現(xiàn)了一種非常規(guī)周期結(jié)構(gòu)—周期泡（這種動(dòng)力學(xué)行為在該系統(tǒng)中不常見也不好發(fā)現(xiàn)，目前鮮有報(bào)道）。為了能夠?qū)υ撓到y(tǒng)的周期泡現(xiàn)象作詳細(xì)分析，用符號(hào)q=p/n表示系統(tǒng)對(duì)應(yīng)的周期碰撞運(yùn)動(dòng)，其中p表示碰撞次數(shù)，n表示周期數(shù)。系統(tǒng)隨參數(shù)ω變化時(shí)相位分岔圖如圖3(a)所示，隨ω減小，當(dāng)ω=2.005 2時(shí)，系統(tǒng)由1/2周期運(yùn)動(dòng)倍化為2/4周期運(yùn)動(dòng)；當(dāng)ω=1.935 4時(shí)，進(jìn)入混沌運(yùn)動(dòng)；當(dāng)ω=1.806 8時(shí)，系統(tǒng)運(yùn)動(dòng)變?yōu)?/2周期運(yùn)動(dòng)；當(dāng)ω=1.804 4時(shí)，倍化為4/4周期運(yùn)動(dòng)，系統(tǒng)呈現(xiàn)出周期泡形態(tài)，最后由逆倍化分岔序列退化為1/1周期運(yùn)動(dòng)。分岔圖與最大Lyapunov指數(shù)譜的對(duì)比圖如圖3(b)所示（灰色的實(shí)線為最大Lyapunov指數(shù)譜，虛線為零刻度線），最大Lyapunov指數(shù)譜如圖3(c)所示，相應(yīng)的各個(gè)Lyapunov指數(shù)譜的變化如圖3(d)所示（λ5與λ6接近）。綜合圖3(a)、圖3(b)和圖3(c)可知，系統(tǒng)處于周期運(yùn)動(dòng)時(shí)最大Lyapunov指數(shù)小于零，處于穩(wěn)定狀態(tài)；系統(tǒng)處于混沌運(yùn)動(dòng)時(shí)，最大Lyapunov指數(shù)大于零，處于不穩(wěn)定狀態(tài)；當(dāng)系統(tǒng)周期運(yùn)動(dòng)發(fā)生分岔時(shí)(ω=2.005 2，1.804 4，1.765 2，1.686)，系統(tǒng)的最大Lyapunov指數(shù)為零。

圖3 系統(tǒng)的分岔圖和Lyapunov指數(shù)譜

現(xiàn)在來進(jìn)一歩分析圖(3)系統(tǒng)的動(dòng)力學(xué)行為與穩(wěn)定性。圖4(a)為質(zhì)塊1在ω=1.776時(shí)對(duì)應(yīng)的4/4周期運(yùn)動(dòng)相圖，圖4(b)為該頻率下的Lyapunov指數(shù)譜的收斂情況（虛線為零刻度線，λ5與λ6接近），可以得到4/4周期運(yùn)動(dòng)對(duì)應(yīng)的最大Lyapunov指數(shù)小于零。圖4(c)為質(zhì)塊1在ω=1.9對(duì)應(yīng)混沌運(yùn)動(dòng)相圖。圖4(d)為該頻率下的Lyapunov指數(shù)譜的收斂情況（λ3與λ4接近，λ5與λ6接近），可以得到混沌運(yùn)動(dòng)對(duì)應(yīng)的最大Lyapunov指數(shù)大于零。

通過數(shù)值仿真發(fā)現(xiàn)，該系統(tǒng)的周期泡現(xiàn)象對(duì)參數(shù)m3很敏感，參數(shù)m3的微小變動(dòng)使得系統(tǒng)的周期泡數(shù)量變化甚至消失。為了更好的了解參數(shù)m3對(duì)系統(tǒng)動(dòng)力學(xué)行為的影響，取系統(tǒng)參數(shù)m3=1.7，m3=1.68，m3=1.62，m3=1.61，m3=1.607，m3=1.606，m3=1.601 4，其余參數(shù)保持不變，系統(tǒng)隨參數(shù)ω變化時(shí)速度分岔圖如圖5所示。

隨著ω的增大，圖5(a)所示系統(tǒng)從1/1周期運(yùn)動(dòng)倍化為2/2周期運(yùn)動(dòng)，隨后進(jìn)入混沌運(yùn)動(dòng)；圖5(b)所示系統(tǒng)從1/1周期運(yùn)動(dòng)倍化為2/2周期運(yùn)動(dòng)，再從2/2周期運(yùn)動(dòng)倍化為4/4周期運(yùn)動(dòng)，隨后由4/4周期運(yùn)動(dòng)退化為2/2周期運(yùn)動(dòng)，最后進(jìn)入混沌運(yùn)動(dòng)，系統(tǒng)呈現(xiàn)出兩個(gè)周期泡結(jié)構(gòu)；圖5(c)所示系統(tǒng)從1/1周期運(yùn)動(dòng)倍化為2/2周期運(yùn)動(dòng)，再從2/2周期運(yùn)動(dòng)倍化為4/4周期運(yùn)動(dòng)，然后從4/4周期運(yùn)動(dòng)倍化為8/8周期運(yùn)動(dòng)，接著經(jīng)逆倍化分岔序列退化為2/2周期運(yùn)動(dòng)，系統(tǒng)呈現(xiàn)出四個(gè)周期泡結(jié)構(gòu)；圖5(d)所示系統(tǒng)從1/1周期運(yùn)動(dòng)經(jīng)倍周期分岔序列倍化為16/16周期運(yùn)動(dòng)，接著從16/16周期運(yùn)動(dòng)經(jīng)逆倍化分岔序列退化為2/2周期運(yùn)動(dòng)，系統(tǒng)呈現(xiàn)8個(gè)周期泡結(jié)構(gòu)；圖5(e)所示系統(tǒng)從1/1周期運(yùn)動(dòng)經(jīng)倍周期分岔序列倍化為32/32周期運(yùn)動(dòng)，然后又從32/32周期運(yùn)動(dòng)經(jīng)逆倍化分岔序列退化為2/2周期運(yùn)動(dòng)，系統(tǒng)呈現(xiàn)16個(gè)周期泡結(jié)構(gòu)；圖5(f)所示系統(tǒng)從1/1周期運(yùn)動(dòng)經(jīng)倍周期分岔序列進(jìn)入混沌運(yùn)動(dòng)，隨后又從混沌運(yùn)動(dòng)經(jīng)逆倍化分岔序列退化為2/2周期運(yùn)動(dòng)，系統(tǒng)呈現(xiàn)16個(gè)混沌泡結(jié)構(gòu)；圖5(h)為圖5(g)的局部放大圖，圖5(g)所示系統(tǒng)從1/1周期運(yùn)動(dòng)經(jīng)倍周期分岔序列進(jìn)入混沌運(yùn)動(dòng)，隨后從混沌運(yùn)動(dòng)進(jìn)入周期運(yùn)動(dòng)，接著再次進(jìn)入混沌運(yùn)動(dòng)，然后又從混沌運(yùn)動(dòng)進(jìn)入周期運(yùn)動(dòng)，最后進(jìn)入混沌運(yùn)動(dòng)后經(jīng)逆倍化分岔序列進(jìn)入2/2周期運(yùn)動(dòng)。

圖4 質(zhì)塊1在ω=1.776與ω=1.9的相圖及Lyapunov指數(shù)譜收斂序列

進(jìn)一歩對(duì)圖5(d)所示系統(tǒng)出現(xiàn)的周期泡現(xiàn)象作詳細(xì)分析。圖6為圖5(d)所示系統(tǒng)中質(zhì)塊1在不同頻率下的相圖。圖6(a)為質(zhì)塊1的1/1周期運(yùn)動(dòng)對(duì)應(yīng)的相圖；圖6(b)為質(zhì)塊1的2/2周期運(yùn)動(dòng)對(duì)應(yīng)的相圖；圖6(c)為質(zhì)塊1的4/4周期運(yùn)動(dòng)對(duì)應(yīng)的相圖；圖6(d)為質(zhì)塊1的8/8周期運(yùn)動(dòng)對(duì)應(yīng)的相圖；圖6(e)為質(zhì)塊1的16/16周期運(yùn)動(dòng)對(duì)應(yīng)的相圖。

圖5 系統(tǒng)隨m3變化的速度分岔圖及局部放大圖

進(jìn)一歩對(duì)圖5(f)所示系統(tǒng)的混沌泡現(xiàn)象作具體分析，圖5(f)的局部分岔圖如圖7(a)所示，對(duì)應(yīng)的最大Lyapunov指數(shù)譜如圖7(b)所示。由圖7(b)可知，在系統(tǒng)周期運(yùn)動(dòng)發(fā)生分岔處（ω=1.761 12，1.773 36，1.778 4，1.780 32，1.785 6，1.788 24，1.794 48）對(duì)應(yīng)的最大Lyapunov指數(shù)為零。當(dāng)系統(tǒng)進(jìn)入混沌運(yùn)動(dòng)，即出現(xiàn)混沌泡現(xiàn)象時(shí)，對(duì)應(yīng)的Lyapunov指數(shù)大于零，系統(tǒng)失穩(wěn)。

4 結(jié)語

本文引入局部映射，得到Poincaré映射和Jacobi矩陣，通過Gram-Schmidt正交化和范數(shù)歸一化的方法得到3自由度單碰振動(dòng)系統(tǒng)Lyapunov指數(shù)譜的計(jì)算方法。結(jié)果表明，利用系統(tǒng)的Lyapunov指數(shù)譜，可以有效地對(duì)系統(tǒng)的穩(wěn)定性進(jìn)行分析。通過數(shù)值仿真，結(jié)合系統(tǒng)的分岔圖和相圖，發(fā)現(xiàn)在一定參數(shù)下系統(tǒng)出現(xiàn)周期泡、混沌泡等豐富的動(dòng)力學(xué)行為。在其他參數(shù)保持不變的情況下，隨著質(zhì)量比的減小，系統(tǒng)結(jié)構(gòu)發(fā)生改變，系統(tǒng)運(yùn)動(dòng)的周期泡結(jié)構(gòu)數(shù)量成倍增多，即系統(tǒng)運(yùn)動(dòng)的周期數(shù)在成倍增加，接著周期泡結(jié)構(gòu)消失出現(xiàn)混沌泡結(jié)構(gòu)，最終混沌泡結(jié)構(gòu)消失系統(tǒng)進(jìn)入混沌運(yùn)動(dòng)。

圖6 質(zhì)塊1在不同頻率下的相圖

圖7 分岔圖局部放大圖和Lyapunov指數(shù)譜