任一凡，馮進(jìn)鈐，沈曉娜

(西安工程大學(xué) 理學(xué)院，陜西 西安 710048)

0 引 言

碰撞振動(dòng)現(xiàn)象是生活中常見(jiàn)的一種現(xiàn)象，如列車行駛中車輪與鐵軌之間的碰撞，機(jī)器運(yùn)行中零部件之間的碰撞。這些現(xiàn)象有時(shí)會(huì)給人們的日常生活帶來(lái)負(fù)面影響，所以，人們展開(kāi)了對(duì)碰撞振動(dòng)系統(tǒng)的研究，包括系統(tǒng)的周期運(yùn)動(dòng)、分岔和混沌。近年來(lái)，碰撞等非光滑因素的存在通常使得系統(tǒng)的軌線流出現(xiàn)跳躍和零行為等奇異性，使得在碰撞振動(dòng)系統(tǒng)中出現(xiàn)了新的現(xiàn)象。對(duì)于這些新的現(xiàn)象，許多學(xué)者對(duì)此也做了很多的研究[1]。另一方面，由于不連續(xù)特性的存在，適用于光滑系統(tǒng)的數(shù)值積分方法不再直接有效，不連續(xù)時(shí)刻的數(shù)值定位成為碰撞振動(dòng)系統(tǒng)數(shù)值研究中熱點(diǎn)問(wèn)題之一[2]。

對(duì)于非線性系統(tǒng)的全局分析方法分為解析方法和數(shù)值方法，而胞映射方法作為一種分析動(dòng)力系統(tǒng)全局特性的有效數(shù)值方法，運(yùn)行速度快，提高了工作效率[3]。胞映射方法由HUS在1980年首次提出，為研究動(dòng)力系統(tǒng)的復(fù)雜運(yùn)動(dòng)提供了一種新思路[3-4]。該方法將系統(tǒng)離散化為胞，把感興趣的胞空間分割，利用胞之間的轉(zhuǎn)移關(guān)系研究原動(dòng)力系統(tǒng)的動(dòng)力學(xué)行為。隨后，HSU又提出計(jì)算更加準(zhǔn)確的廣義胞映射方法和圖胞映射方法[5]。為了提高胞映射的計(jì)算精度，JIANG等提出了胞參照點(diǎn)映射法[6]。隨后，HONG等又將廣義胞映射圖論方法擴(kuò)展到模糊動(dòng)力系統(tǒng)問(wèn)題的研究中[7]。賀群等引進(jìn)圖論的四元組，構(gòu)建了瞬態(tài)胞分類方法，該方法在逼近動(dòng)力系統(tǒng)的穩(wěn)定流形和不穩(wěn)定流形方面效果良好[8]。傳統(tǒng)的胞映射方法主要針對(duì)光滑系統(tǒng)，在處理非光滑系統(tǒng)時(shí)必須建立合適的數(shù)值積分方法[9]。關(guān)于非光滑系統(tǒng)的胞映射方法，目前已有的研究較少。李爽等對(duì)邊界胞進(jìn)行了重新細(xì)分，提出了非光滑系統(tǒng)迭代圖胞映射方法[10]。文獻(xiàn)[11-13]基于圖胞映射方法，提出了擦邊流形的逼近算法，建立了適用于隨機(jī)碰撞振動(dòng)系統(tǒng)的自適應(yīng)胞映射方法，并討論了系統(tǒng)的全局結(jié)構(gòu)，研究了噪聲對(duì)全局結(jié)構(gòu)的影響。

針對(duì)碰撞振動(dòng)系統(tǒng)不連續(xù)的特性，建立高效的適用于碰撞振動(dòng)系統(tǒng)的數(shù)值積分方法，難點(diǎn)在于解決非光滑系統(tǒng)的積分運(yùn)算[14-16]。本文結(jié)合牛頓迭代法的思想，構(gòu)建了適用于碰撞振動(dòng)系統(tǒng)的數(shù)值積分方法。將該數(shù)值積分方法運(yùn)用于Duffing型碰撞振動(dòng)系統(tǒng)，驗(yàn)證了該數(shù)值方法的有效性與高效性；同時(shí)，將該數(shù)值積分方法融入胞映射算法中，研究了Duffing型碰撞振動(dòng)系統(tǒng)的全局動(dòng)力學(xué)。

1 碰撞振動(dòng)系統(tǒng)的數(shù)值積分方法

考慮一般的二維碰撞振動(dòng)系統(tǒng)，系統(tǒng)的動(dòng)力學(xué)方程為

(1)

(2)

式中:x1=Δ;上標(biāo)-、+分別表示碰撞前后的時(shí)刻；r表示碰撞恢復(fù)系數(shù)(0

1.1 牛頓迭代法定位碰撞時(shí)刻

在研究碰撞振動(dòng)系統(tǒng)的動(dòng)力學(xué)行為時(shí)，由于碰撞的存在，系統(tǒng)的軌線呈現(xiàn)不連續(xù)或不可微的情形，無(wú)法直接使用傳統(tǒng)的數(shù)值積分方法。因此，建立數(shù)值方法準(zhǔn)確定位碰撞時(shí)刻,是處理不連續(xù)性問(wèn)題的關(guān)鍵之一。為了方便描述數(shù)值定位碰撞時(shí)刻，圖1給出了牛頓迭代法示意圖。

圖 1 牛頓迭代法示意圖Fig.1 Schematic diagram of Newton iteration method

將x1(tk+1)利用泰勒公式展開(kāi)，得

(3)

根據(jù)約束條件，在式(3)中令x1(tk+1)=Δ，則

(4)

代入系統(tǒng)(1)，解得

從而得到時(shí)間序列{tk}，k=0，1，2,…的迭代格式

(5)

由式(5)可產(chǎn)生序列{tk}，當(dāng)|tk+1-tk|<ε時(shí)停止迭代，則t*≈tk，從而得出碰撞時(shí)刻t*的近似數(shù)值解。

1.2 碰撞振動(dòng)系統(tǒng)的數(shù)值積分方法

假設(shè)RK積分步長(zhǎng)為h，Δ表示碰撞面，T表示總積分時(shí)間。考慮一般的碰撞振動(dòng)系統(tǒng)(1)、(2)，采取經(jīng)典的RK數(shù)值積分方法求解。數(shù)值解示意圖如圖2所示。

圖 2 數(shù)值解示意圖Fig.2 Schematic diagram of numerical solution

假設(shè)系統(tǒng)從tk時(shí)刻開(kāi)始積分，在tk時(shí)刻的狀態(tài)位置記為(x1(tk),x2(tk))T。對(duì)系統(tǒng)(1)進(jìn)行RK數(shù)值積分，得到下一時(shí)刻的狀態(tài)位置，判斷這一位置系統(tǒng)(1)是否跨過(guò)碰撞面：若未跨過(guò)碰撞面，則繼續(xù)進(jìn)行下一步RK數(shù)值積分，直至跨越碰撞面；否則，利用迭代格式(5)和RK數(shù)值積分，數(shù)值定位碰撞時(shí)刻t*。然后，利用碰撞映射式(2)得到碰撞后狀態(tài)，繼續(xù)進(jìn)行光滑系統(tǒng)(1)的數(shù)值積分。為了更加清晰準(zhǔn)確地描述這一過(guò)程，給出具體的算法步驟如下：

步驟1：輸入初始條件(tk,x1(tk))，k=0。

步驟2：從tk時(shí)刻開(kāi)始積分，RK一步映射得到tk+1時(shí)刻的狀態(tài)x1(tk+1)。

步驟3：判斷這一時(shí)刻的狀態(tài)是否跨過(guò)碰撞面：若跨過(guò)，則轉(zhuǎn)步驟4；若沒(méi)有跨過(guò)，轉(zhuǎn)步驟5。

步驟5：若tk+1>T，轉(zhuǎn)步驟6；否則，令(tk+1,x1(tk+1))→(tk,x1(tk))，轉(zhuǎn)步驟2。

步驟6：輸出結(jié)果。

2 模型算例

考慮諧和激勵(lì)下的Duffing型碰撞振動(dòng)系統(tǒng)，系統(tǒng)的模型如下：

(6)

(7)

式中:a、c為常數(shù)；b為阻尼系數(shù)；x表示為廣義位移；fcos(ωt)為諧和激勵(lì)。

針對(duì)系統(tǒng)(6)分別作周期運(yùn)動(dòng)和混沌運(yùn)動(dòng)的情形，進(jìn)一步驗(yàn)證上述數(shù)值積分方法的有效性。首先，固定部分系統(tǒng)參數(shù)，a=1.0，b=0.2，c=1.0，r=0.8，Δ=1.0，ω=1.0，取定初始值(x0,y0)T=(0.8,0)T，系統(tǒng)外激勵(lì)周期T=2π/ω。

可知系統(tǒng)(6)在f=0.38時(shí)作周期運(yùn)動(dòng)，在f=0.42時(shí)作混沌運(yùn)動(dòng)[17]。利用上述數(shù)值積分方法，圖3給出了系統(tǒng)(6)在諧和幅值f為0.38時(shí)的相圖和時(shí)間歷程圖?？梢钥闯觯琭=0.38時(shí)，系統(tǒng)作規(guī)則的周期運(yùn)動(dòng)。圖4給出了f=0.42時(shí)的相圖，可以看出，原來(lái)的周期運(yùn)動(dòng)消失，此時(shí)系統(tǒng)作混沌運(yùn)動(dòng)。

(a) 相圖

(b) 時(shí)間歷程圖圖 3 當(dāng)f=0.38時(shí)，系統(tǒng)(6)的數(shù)值解Fig.3 The numerical solution of system (6) when f=0.38

圖 4 f=0.42時(shí)系統(tǒng)(6)的相圖Fig.4 Phase diagram of system (6) when f=0.42

可以得出，圖3和4所得數(shù)值結(jié)果與文獻(xiàn)[17]一致。此外，數(shù)值仿真過(guò)程中，迭代格式(5)中的精度取為1×10-8，最大迭代次數(shù)為3，說(shuō)明該數(shù)值積分方法適用于碰撞振動(dòng)系統(tǒng)的周期運(yùn)動(dòng)以及混沌運(yùn)動(dòng)的數(shù)值研究，并且具有較高的精度與穩(wěn)定性。

3 Duffing雙邊碰撞振動(dòng)系統(tǒng)全局結(jié)構(gòu)

為了更深入地分析非光滑因素在系統(tǒng)全局動(dòng)力學(xué)中的影響，重點(diǎn)考慮諧和激勵(lì)f對(duì)系統(tǒng)(6)全局結(jié)構(gòu)的影響。固定部分參數(shù)a=1.0，b=0.2，c=1.0，r=0.8，Δ=1，ω=1.0；選取龐加萊截面

Π={(x,y,θ)∈R2×S:θ=

ωtmod(2π)}

取狀態(tài)空間的不同初始值(0.75,-0.3)，(-0.75,0.3)，(0,0.8)，(0,-0.7)和(0.8,0)做龐加萊映射。去掉前100個(gè)瞬態(tài)周期點(diǎn)后，保留500個(gè)穩(wěn)態(tài)周期點(diǎn)，得到多初值分岔圖，如圖5所示。

從圖5可以看出:當(dāng)f=0.38時(shí)，系統(tǒng)存在5個(gè)周期為1的共存吸引子;f增加到0.392，其中有3個(gè)周期解由倍周期分岔通往混沌，2個(gè)周期解的拓?fù)湫再|(zhì)保持不變，此時(shí)，系統(tǒng)呈現(xiàn)出周期解與混沌解共存的現(xiàn)象; 當(dāng)f∈[0.392,0.395]時(shí)，系統(tǒng)的混沌解逐漸消失，退化為3個(gè)周期為1的共存解;隨著f逐漸增加到0.409，此時(shí)系統(tǒng)依舊是3個(gè)周期解共存。f繼續(xù)增加至0.41，2個(gè)周期為1的周期解發(fā)生激變，直接進(jìn)入混沌運(yùn)動(dòng)。

圖 5 多初值分岔圖Fig.5 Multi-initial bifurcation diagram

為了更深入地分析該系統(tǒng)的動(dòng)力學(xué)演化，利用簡(jiǎn)單胞映射方法[17]探討系統(tǒng)激發(fā)生前后的全局動(dòng)力學(xué)行為。將感興趣的區(qū)域

D={(x,y)|-1≤x≤1,-2≤y≤2}

劃分成300×400個(gè)大小相同的格子，每個(gè)格子記為正規(guī)胞，區(qū)域D以外的區(qū)域記為陷胞，數(shù)值積分采用上述牛頓迭代積分法。

當(dāng)f=0.38時(shí)，系統(tǒng)(6)的全局圖如6(a)所示，其中Ai(i=1,2,3,4,5)表示與吸引子，Bi(i=1,2,3,4,5)分別表示與吸引子Ai對(duì)應(yīng)的吸引域。此時(shí)系統(tǒng)存在5個(gè)吸引子，對(duì)應(yīng)5個(gè)吸引域。為了更清晰地看出系統(tǒng)的運(yùn)動(dòng)狀態(tài)，吸引子對(duì)應(yīng)的相圖如圖6(b)所示，圖6(b)由左上位置開(kāi)始分別對(duì)應(yīng)吸引子A4、A5、A1、A3、A2。可以看出，此時(shí)系統(tǒng)存在5個(gè)周期為1的運(yùn)動(dòng)，系統(tǒng)呈現(xiàn)出5種不同的周期運(yùn)動(dòng)，并呈現(xiàn)出一定的對(duì)稱性，與圖6(a)的全局結(jié)構(gòu)圖吻合。

(a) 全局結(jié)構(gòu)圖

(b) 相圖圖 6 當(dāng)f=0.38時(shí)，系統(tǒng)(6)的數(shù)值結(jié)果Fig.6 The numerical solution of system (6) When f=0.38

隨著諧和激勵(lì)f逐漸增大到0.392，系統(tǒng)(6)的全局圖如圖7(a)所示。從圖7(a)可以看出：吸引子A1、A2和A3的拓?fù)湫再|(zhì)沒(méi)有變化，而周期吸引子A4和A5消失，突變?yōu)榛煦缥?。此時(shí)，系統(tǒng)呈現(xiàn)出周期解與混沌解共存的現(xiàn)象。當(dāng)f=0.394時(shí)，系統(tǒng)(6)的全局圖如圖7(b)所示。從圖7(b)可知：吸引子A1和A3及其吸引域B1和B3一同消失，原先的混沌吸引子A4和A5退化成2個(gè)共存周期吸引子，系統(tǒng)呈現(xiàn)出3個(gè)周期吸引子共存的現(xiàn)象，系統(tǒng)經(jīng)歷了混沌解到周期解的激變。

繼續(xù)增大f，在f=0.409時(shí)，系統(tǒng)依舊呈現(xiàn)出3個(gè)周期為1的吸引子共存的現(xiàn)象，見(jiàn)圖7(c)。隨著諧和激勵(lì)f繼續(xù)增加到0.41，此時(shí)系統(tǒng)的全局結(jié)構(gòu)如圖7(d)所示。從圖7(d)可以看出：周期吸引子A2的拓?fù)湫再|(zhì)沒(méi)有變化，周期吸引子A4和A5突變?yōu)?個(gè)混沌吸引子，系統(tǒng)呈現(xiàn)出周期吸引子與混沌吸引子共存的現(xiàn)象，系統(tǒng)發(fā)生了周期解到混沌解的激變。

(a) f=0.392

(b) f=0.395

(c) f=0.409

(d) f=0.41圖 7 系統(tǒng)(6)的全局結(jié)構(gòu)圖Fig.7 Global structure diagram of system (6)

4 結(jié) 語(yǔ)

本文構(gòu)建了適用于碰撞振動(dòng)系統(tǒng)的牛頓迭代積分法，并驗(yàn)證了該數(shù)值積分方法的有效性及高效性。同時(shí)，將該方法運(yùn)用于胞映射算法中，研究了雙邊約束下Duffing型碰撞振動(dòng)系統(tǒng)的全局結(jié)構(gòu)。結(jié)果表明，改進(jìn)后的數(shù)值積分方法適用于系統(tǒng)的周期運(yùn)動(dòng)和混沌運(yùn)動(dòng)的數(shù)值研究，并提高了計(jì)算速度。通過(guò)對(duì)系統(tǒng)全局動(dòng)力學(xué)的研究發(fā)現(xiàn)，系統(tǒng)存在著豐富的多吸引子結(jié)構(gòu)，如多周期吸引子、多混沌吸引子的共存現(xiàn)象。此外，諧和激勵(lì)幅值的變化可以誘導(dǎo)系統(tǒng)發(fā)生多解激變現(xiàn)象，包括周期解直接進(jìn)入混沌、混沌激變等。