任一凡，馮進鈐，沈曉娜

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

0 引 言

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

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

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

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

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

(1)

(2)

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

1.1 牛頓迭代法定位碰撞時刻

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

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

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

(3)

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

(4)

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

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

(5)

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

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

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

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

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

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

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

步驟3：判斷這一時刻的狀態(tài)是否跨過碰撞面：若跨過，則轉步驟4；若沒有跨過，轉步驟5。

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

步驟6：輸出結果。

2 模型算例

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

(6)

(7)

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

針對系統(tǒng)(6)分別作周期運動和混沌運動的情形，進一步驗證上述數(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)外激勵周期T=2π/ω。

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

(a) 相圖

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

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

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

3 Duffing雙邊碰撞振動系統(tǒng)全局結構

為了更深入地分析非光滑因素在系統(tǒng)全局動力學中的影響，重點考慮諧和激勵f對系統(tǒng)(6)全局結構的影響。固定部分參數(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個瞬態(tài)周期點后，保留500個穩(wěn)態(tài)周期點，得到多初值分岔圖，如圖5所示。

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

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

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

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

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

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

(a) 全局結構圖

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

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

繼續(xù)增大f，在f=0.409時，系統(tǒng)依舊呈現(xiàn)出3個周期為1的吸引子共存的現(xiàn)象，見圖7(c)。隨著諧和激勵f繼續(xù)增加到0.41，此時系統(tǒng)的全局結構如圖7(d)所示。從圖7(d)可以看出：周期吸引子A2的拓撲性質沒有變化，周期吸引子A4和A5突變?yōu)?個混沌吸引子，系統(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)的全局結構圖Fig.7 Global structure diagram of system (6)

4 結 語

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