碰撞振動系統的牛頓迭代積分法與全局動力學

任一凡，馮進鈐，沈曉娜

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

0 引 言

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

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

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

1 碰撞振動系統的數值積分方法

考慮一般的二維碰撞振動系統，系統的動力學方程為

(1)

(2)

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

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

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

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

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

(3)

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

(4)

代入系統(1)，解得

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

(5)

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

1.2 碰撞振動系統的數值積分方法

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

圖 2 數值解示意圖Fig.2 Schematic diagram of numerical solution

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

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

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

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

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

步驟6：輸出結果。

2 模型算例

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

(6)

(7)

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

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

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

(a) 相圖

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

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

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

3 Duffing雙邊碰撞振動系統全局結構

為了更深入地分析非光滑因素在系統全局動力學中的影響，重點考慮諧和激勵f對系統(6)全局結構的影響。固定部分參數a=1.0，b=0.2，c=1.0，r=0.8，Δ=1，ω=1.0；選取龐加萊截面

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

ωtmod(2π)}

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

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

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

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

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

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

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

(a) 全局結構圖

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

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

繼續增大f，在f=0.409時，系統依舊呈現出3個周期為1的吸引子共存的現象，見圖7(c)。隨著諧和激勵f繼續增加到0.41，此時系統的全局結構如圖7(d)所示。從圖7(d)可以看出：周期吸引子A2的拓撲性質沒有變化，周期吸引子A4和A5突變為2個混沌吸引子，系統呈現出周期吸引子與混沌吸引子共存的現象，系統發生了周期解到混沌解的激變。

(a) f=0.392

(b) f=0.395

(c) f=0.409

(d) f=0.41圖 7 系統(6)的全局結構圖Fig.7 Global structure diagram of system (6)

4 結 語

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