非均質(zhì)物理雙擺的混沌特性研究

孔令輝 劉丁楊 蹇開林

摘要：為了解決工程實際中材料質(zhì)量不均勻分布對雙擺系統(tǒng)運(yùn)動的影響，在均質(zhì)物理雙擺模型的基礎(chǔ)上，將擺的質(zhì)心位置和擺的轉(zhuǎn)動慣量提取為變量，建立非均質(zhì)雙擺模型。將非均質(zhì)雙擺系統(tǒng)由Hamilton系統(tǒng)近似為擬Hamilton系統(tǒng)，運(yùn)用雙自由度的Melnikov法，得到擬Hamilton系統(tǒng)存在Smale馬蹄意義下混沌的能量閾值，以此作為Hamilton系統(tǒng)的混沌條件。利用最大Lyapunov指數(shù)圖、分岔圖、Poincaré截面圖等數(shù)值方法驗證混沌條件的正確性，并詳細(xì)分析了各參數(shù)對系統(tǒng)運(yùn)動狀態(tài)的影響和作用機(jī)制。結(jié)果表明，非均質(zhì)雙擺的混沌閾值有較高復(fù)雜性，而且擺長、擺重、第一擺的質(zhì)心位置同時影響著系統(tǒng)的能量與混沌閾值，解釋了質(zhì)心位置和轉(zhuǎn)動慣量等參數(shù)發(fā)生變化時，系統(tǒng)在混沌和擬周期之間交替變換的原因。進(jìn)一步研究了參數(shù)取值與Melnikov法適用性之間的關(guān)系，通過數(shù)值仿真分類討論了Melnikov法不適用時的參數(shù)取值情況。

關(guān)鍵詞：非均質(zhì)物理雙擺；混沌；擬周期；Melnikov法

中圖分類號：O322????????? 文獻(xiàn)標(biāo)志碼：A????? 文章編號：1000-582X（2024）02-106-13

Study on chaotic characteristics of heterogeneous physical double pendulum

KONG Linghui1a， LIU Dingyang2， JIAN Kailin1a，1b

（1a. School of Aeronautics and Astronautics; 1b. Chongqing Key Laobratory of Heterogenous Materisal Mechanics， Chongqing University， Chongqing 400044， P. R. China; 2. CSSC Haizhuang Windpower Co.， Ltd.， Chongqing 401122， P. R. China）

Abstract： To address the impact of materials with uneven mass distribution on the motion of double pendulum systems in engineering practice， a heterogeneous double pendulum model was established based on the homogeneous physical double pendulum model. This model incorporated variables， such as the position of the center of mass and the moment of inertia of the pendulum. To further explore the chaotic characteristics of the system， the heterogeneous double pendulum system was approximated from a Hamilton system to a quasi-Hamilton system. The dynamics equation of the double pendulum system was obtained using the Euler-Lagrange equation of the second kind. The energy threshold for Smale horseshoe chaos in the quasi-Hamiltonian system was determined using Melnikov method with two degrees of freedom， serving as the chaos condition for the Hamiltonian system. Through programming in Matlab， the correctness of the chaos condition was verified using numerical methods such as maximum Lyapunov exponential diagram， bifurcation diagram and Poincare section diagram. The influence of each parameter on the systems motion state and action mechanism was analyzed in detail. The results show that the chaos threshold of heterogeneous double pendulums depends on various factors， including the initial energy of the double pendulums， the position of the center of mass of the first pendulums and the ratio of the moment of inertia of the two pendulums. With different values for each parameter， the system changes from a regular motion state with a chaos threshold to an irregular and complex motion state without a chaos threshold. The reasons for the system alternating between chaotic and quasi-periodic states were explained， and theoretical predictions were validated. The differences between the theoretical threshold and actual numerical simulation results were explained when the center of mass of the first pendulum， the mass ratio of the two pendulums， and the ratio of moment of inertia of the two pendulums were set to their limit values. On this basis， the relationship between parameter values and the applicability of Melnikov method were further explored， and the parameters under which Melnikov method was no longer applicable were discussed by numerical simulation classification.

Keywords： heterogeneous physical double pendulum; chaos; quasi-period; Melnikov method

雙擺廣泛地應(yīng)用于生活、工程與科學(xué)中。生活中，人的身體擺動[1]、高爾夫球揮桿[2]與網(wǎng)球揮拍[3]等常見動作以及工程上的機(jī)械臂[4]、機(jī)械足[5]、各種類型的起重機(jī)[6?8]、隔振器[9]、能量收集裝置[10]都可以簡化為雙擺模型。在科學(xué)上，范洪義等[11]通過擺、小車和彈簧的互相牽制效應(yīng)研究了量子糾纏的經(jīng)典類比問題。陳漢軍等[12]通過對雙物理擺混沌結(jié)果的研究，給一種新型振動機(jī)的參數(shù)設(shè)計提供了理論依據(jù)。Ford等[13]相信雙擺作為一個簡單的二自由度模型，很適用于測試量子力學(xué)能否描述經(jīng)典混沌系統(tǒng)的觀測實驗，對經(jīng)典混沌的量子對應(yīng)研究意義重大。因此，雙擺的混沌問題備受關(guān)注。

Stachowiak等[14]通過Poincaré截面圖研究平面數(shù)學(xué)雙擺，研究表明系統(tǒng)的運(yùn)動狀態(tài)與能量密切相關(guān)，在零能量極限時，雙擺做周期運(yùn)動，隨著能量的增加，雙擺將從周期運(yùn)動轉(zhuǎn)變?yōu)闇?zhǔn)周期運(yùn)動，進(jìn)而轉(zhuǎn)變?yōu)榛煦邕\(yùn)動。趙武等[15]、劉丁楊等[16]和Han等[17]應(yīng)用Melnikov法研究了不同復(fù)雜擺系統(tǒng)的混沌閾值，并通過數(shù)值模擬加以驗證。Maiti等[18]研究了旋轉(zhuǎn)雙擺的混沌特性，發(fā)現(xiàn)系統(tǒng)的內(nèi)共振將導(dǎo)致混沌。Bi等[19]使用Normal Form理論分析了自治雙擺系統(tǒng)在1：1內(nèi)共振時的分岔特性，研究表明當(dāng)存在周期激勵時，系統(tǒng)會按周期倍化的道路通向混沌；沒有周期激勵時，參數(shù)的改變也會影響其解的穩(wěn)定性。Martynyuk等[20]利用動態(tài)對稱原理分析對雙擺產(chǎn)生混沌運(yùn)動的條件，證明了雙擺質(zhì)量比較大時，雙擺存在有條件的周期性軌跡和混沌軌跡。Calv?o等[21]不僅比較了分岔圖、最大Lyapunov指數(shù)譜圖、功率譜圖、傅里葉變換圖、時間歷程圖等不同方法在混沌研究中的優(yōu)缺點(diǎn)，而且研究了兩擺初始角度對雙擺混沌的影響。Oiwa等[22]通過Kosambi-Cartan-Chern理論研究了雙擺的雅可比穩(wěn)定性，給出了雅可比穩(wěn)定性條件，通過Lyapunov指數(shù)和Poincaré截面圖，發(fā)現(xiàn)雅可比不穩(wěn)定區(qū)域與混沌行為的開始有關(guān)。Kovacic等[23]研究了物理雙擺在低能級時的非線性正則模態(tài)，通過諧波平衡法和Lindstedt-Poincaré法分析了系統(tǒng)非線性振動的情況，結(jié)合數(shù)值仿真比較2種方法在不同情況下的優(yōu)缺點(diǎn)。Dudkowski等[24]通過數(shù)值方法和實驗研究了具有參數(shù)激勵的機(jī)械雙擺，研究表明系統(tǒng)不規(guī)則運(yùn)動的距離對參數(shù)和初始條件都非常敏感，而動態(tài)響應(yīng)的特性嚴(yán)格依賴于激勵參數(shù)。

綜上所述，對于雙擺大都集中在采用Lyapunov指數(shù)、分岔圖和Poincaré圖等數(shù)值方法研究系統(tǒng)的混沌性，一部分理論是在研究雙擺穩(wěn)定性之后，再考慮系統(tǒng)穩(wěn)定性和混沌性的關(guān)系，少有對雙擺混沌性的理論研究以及能判斷系統(tǒng)出現(xiàn)混沌的理論閾值的研究。筆者考慮實際工程中材料非均質(zhì)的特點(diǎn)，建立非均質(zhì)物理雙擺模型，應(yīng)用二自由度Melnikov法[25]研究非均質(zhì)雙擺系統(tǒng)的混沌性質(zhì)，從理論上得到系統(tǒng)的混沌閾值，更加深入地探究雙擺轉(zhuǎn)動慣量對混沌閾值的影響，揭示雙擺系統(tǒng)在各參數(shù)變化下誘發(fā)能量改變來對系統(tǒng)運(yùn)動狀態(tài)產(chǎn)生的影響，并通過數(shù)值方法驗證其正確性。

1 非均質(zhì)物理雙擺的動力學(xué)方程

圖1為非均質(zhì)物理雙擺模型。

該模型僅由2桿連接而成，假設(shè)2桿為剛性非均質(zhì)直桿，建立如圖1所示的平面直角坐標(biāo)系，記第二擺端點(diǎn)為P點(diǎn)，其坐標(biāo)表示為（x，y）；2桿與豎直面的夾角為θ1、θ2，2桿角速度為θ ˙_1、θ ˙_2；桿長分別為l1、l2，假定鉸接處到質(zhì)心的距離分別為l_1/k_1? 、 l_2/k_2 ，其中（k1，k2）∈（1，+∞）；質(zhì)量大小分別為m1、m2，繞質(zhì)心軸的轉(zhuǎn)動慣量分別為J1、J2；取重力加速度g=9.81 m/s2，其余參數(shù)的單位與g的單位相一致，系統(tǒng)總動能為

T=1/2 （m_1/（k_1^2 ）+m_2 ） l_1^2 θ ˙_1^2+1/2 J_1 θ ˙_1^2+1/2 J_2 θ ˙_2^2+1/2 m_2 [（l_2^2）/（k_2^2 ） θ ˙_2^2+（2l_1 l_2 θ ˙_1^? θ ˙_2^ ）/k_2? cos（θ_1-θ_2 ） ] 。?? （1）

系統(tǒng)總勢能為

V=（m_1/k_1 +m_2 ）gl_1 （1-cosθ_1 ）+1/k_2? m_2 gl_2 （1-cosθ_2 ） 。??? （2）

Euler-Lagrange第二類方程為

d/dt （（?L_g）/（?θ ˙_j^? ））-（?L_g）/（?θ_j ）=0， （3）

式中，L_g=T-V，因此，由式（3）得非均質(zhì)物理雙擺的微分方程組

{（（m_1/（k_1^2 ）+m_2 ） l_1^2 θ ¨_1^ +（m_1/k_1 +m_2 ）gl_1 sinθ_1+J_1 θ ¨_1+（m_2 l_1 l_2）/k_2? [θ ¨_2 cos（θ_1-θ_2 ）+θ ˙_2^2 sin（θ_1-θ_2 ） ]=0，@（m_2 l_2^2）/（k_2^2 ） θ ¨_2+1/k_2? m_2 gl_2 sinθ_2+J_2 θ ¨_2+（m_2 l_1 l_2）/k_2? [θ ¨_1 cos（θ_1-θ_2 ）-θ ˙_1^2 sin（θ_1-θ_2 ） ]=0 。）┤???? （4）

將式（4）化為一階微分方程組，并令m_1/m_2 =m， l_1/k_1 =L_1， l_2/k_2 =L_2， L_1/L_2 =L， （J_1 k_1^2）/（m_1 l_1^2 ）=a， （J_2 k_2^2）/（m_2 l_2^2 ）=b，則式（4）可以改寫為式（5）。

{（θ ˙_1=w_1 "" ，@θ ˙_2=w_2 "" ，@w ˙_1=（（1+b）[L_2 sin（θ_1-θ_2 ） w_2^2+（1+m/k_1 ）gsinθ_1 ]+cos（θ_1-θ_2 ）[k_1 L_1 sin（θ_1-θ_2 ） w_1^2-gsinθ_2 ]）/（L_1 {k_1 cos^2 （θ_1-θ_2 ）-（1+b）[（a+1）m/k_1 +k_1 ] } ），@w ˙_2=（[（a+1）m/（k_1^2 ）+1]（gsinθ_2-k_1 L_1 sin（θ_1-θ_2 ） w_1^2 ）-cos（θ_1-θ_2 ）[L_2 sin（θ_1-θ_2 ） w_2^2+（1+m/k_1 ）gsinθ_1 ]）/（L_2 {cos^2 （θ_1-θ_2 ）-（1+b）[（a+1）m/（k_1^2 ）+1] } ） 。）┤???? （5）

2 Melnikov法分析非均質(zhì)物理雙擺

Melnikov法是通過構(gòu)建擬Hamilton系統(tǒng)的Melnikov函數(shù)來獲取系統(tǒng)混沌閾值的解析方法。雙擺系統(tǒng)的Hamilton函數(shù)就是系統(tǒng)的機(jī)械能函數(shù)，則雙擺系統(tǒng)的Hamilton函數(shù)H可以表示為式（6）。

（H=1/2 （m_1/（k_1^2 ）+m_2 ） l_1^2 θ ˙_1^2+1/2 m_2 （（l_2^2）/（k_2^2 ） θ ˙_2^2+（2l_1 l_2 θ ˙_1 θ ˙_2）/k_2? cos（θ_1-θ_2 ） ）+@?????? 1/2 J_1 θ ˙_1^2+1/2 J_2 θ ˙_2^2+1/k_2? m_2 gl_2 （1-cosθ_2 ）+（m_1/k_1 +m_2 ）gl_1 （1-cosθ_1 ）。）?? （6）

對式（6）中的余弦項cos（θ_1-θ_2 ）進(jìn)行放縮處理，取cos（θ_1-θ_2 ）/（k_2 g）為ε_1；對cosθ_2進(jìn)行泰勒展開有cosθ_2=1-（θ_2^2）/2！+cosξ/4！ θ_2^4，取Lagrange余項的常數(shù)項cos（ξ）/24 （0<ξ<θ_2 ）為ε_2，因此式（6）化為式（7）。

（H=1/2 （m_1/（k_1^2 ）+m_2 ） l_1^2 θ ˙_1^2+1/2 J_1 θ ˙_1^2+1/2 J_2 θ ˙_2^2+1/k_2? m_2 gl_2 （1/2 θ_2^2-ε_2 θ_2^4 ）+@?????? 1/2 m_2 （（l_2^2）/（k_2^2 ） θ ˙_2^2+2ε_1 l_1 l_2 gθ ˙_1 θ ˙_2 ）+（m_1/k_1 +m_2 ）gl_1 （1-cosθ_1 ），）?? （7）

式中，ε_1 、ε_2都遠(yuǎn)小于1，可以看作微小擾動，且ε_1 、ε_2數(shù)值上相差不大，為方便計算，統(tǒng)一用ε代替。因此，擬Hamilton系統(tǒng)能量函數(shù)H^e可整理為兩相平面函數(shù)F1（θ1，θ ˙_1）、F2（θ2，θ ˙_2）與微擾函數(shù)H^'的和，如式（8）所示。

{（H^e=F_1 （θ_1，θ ˙_1 ）+F_2 （θ_2，θ ˙_2 ）+εH^'，@F_1 （θ_1，θ ˙_1 ）=1/2 （m_1/（k_1^2 ）+m_2 ） l_1^2 θ ˙_1^2+1/2 J_1 θ ˙_1^2+（m_1/k_1 +m_2 ）gl_1 （1-cosθ_1 ），@F_2 （θ_2，θ ˙_2 ）=1/2? （m_2 l_2^2）/（k_2^2 ） θ ˙_2^2+1/2? （m_2 gl_2）/k_2? θ_2^2+1/2 J_2 θ ˙_2^2，@H^'=m_2 l_1 l_2 gθ ˙_1 θ ˙_2-1/k_2? m_2 gl_2 θ_2^4 。）┤? （8）

由F1（θ1，θ ˙_1）可知，相平面（θ1，θ ˙_1）存在中心（0，0），在θ_1∈（-2π，2π）上有雙曲鞍點(diǎn)（-π，0）（π，0），存在2條經(jīng)過2點(diǎn)的異宿軌道，將（-π，0）（π，0）代入F1（θ1，θ ˙_1）得對應(yīng)的Hamilton量為h_1=2（m_1/k_1 +m_2 ）gl_1。

引入變量（I，t），做如式（9）的變換（θ2，θ ˙_2）?（I，t）。將式（9）分別代入式（8）的F2（θ2，θ ˙_2）和H'中，得到式（10），可見相平面（θ2，θ ˙_2）中僅存在以I為參數(shù)的周期軌道。

{（θ_2=√（（2k_2 I）/（m_2 l_2 g） √（（m_2 l_2 k_2 g）/（m_2 l_2^2+k_2^2 J_2 ）））? sin（√（（m_2 l_2 k_2 g）/（m_2 l_2^2+k_2^2 J_2 ）） t），@θ ˙_2=√（（2k_2^2 I）/（m_2 l_2^2+k_2^2 J_2 ） √（（m_2 l_2 k_2 g）/（m_2 l_2^2+k_2^2 J_2 ）））? cos（√（（m_2 l_2 k_2 g）/（m_2 l_2^2+k_2^2 J_2 ）） t），）┤ （9）

{（F_2 （θ_2，θ ˙_2 ）=G（I）=√（（m_2 l_2 k_2 g）/（m_2 l_2^2+k_2^2 J_2 ）） I ，@H^'=m_2 l_1 l_2 k_2 gθ ˙_1 √（2I/（m_2 l_2^2+k_2^2 J_2 ） √（（m_2 l_2 k_2 g）/（m_2 l_2^2+k_2^2 J_2 ））） ?cos（√（（m_2 l_2 k_2 g）/（m_2 l_2^2+k_2^2 J_2 ）） t）-（4k_2^2 I^2）/（m_2 l_2^2+k_2^2 J_2 ） sin^4 （√（（m_2 l_2 k_2 g）/（m_2 l_2^2+k_2^2 J_2 ）） t）。）┤???? （10）

若令1-cosθ_1=1/2 θ_1^2-1/24 θ_1^4+o（θ_1^6 ），略去Peano余項后代入式（8）中的F1（θ1，θ ˙_1），則（θ1，θ ˙_1）上的異宿軌道參數(shù)化為

{（θ_10=±√6 tanh（√（（（m_1 k_1+m_2 k_1^2 ）gl_1）/2（（m_1+m_2 k_1^2 ） l_1^2+J_1 k_1^2 ） ） t），@θ ˙_10=±√（（3（m_1 k_1+m_2 k_1^2 ）gl_1）/（（m_1+m_2 k_1^2 ） l_1^2+J_1 k_1^2 ）） sech^2 （√（（（m_1 k_1+m_2 k_1^2 ）gl_1）/2（（m_1+m_2 k_1^2 ） l_1^2+J_1 k_1^2 ） ） t）。）┤?? （11）

當(dāng)系統(tǒng)的Hamilton量Hh_1時，相平面（θ1，θ ˙_1）上存在不穩(wěn)定流形。穩(wěn)定流形和不穩(wěn)定流形之間的距離由Melnikov函數(shù)[26]表示，如式（12）所示。

（M_± （t_0 ）=∫_（-∞）^（+∞）?{F_1，H^' }? （t-t_0 ）dt=∫_（-∞）^（+∞）?∑_（i=1）^2?（（?F_1）/（?θ_i ）? （?H^'）/（?θ ˙_i ）-（?F_1）/（?θ ˙_i ）? （?H^'）/（?θ_i ））? （t-t_0 ）dt=@??????????????? ?2 （（m_1/k_1 +m_2）m_2^2 l_1 l_2^2 k_2^2 g^2 π）/（（m_2 l_2^2+k_2^2 J_2 ）^2 （m_1 k_1+m_2 k_1^2 ） ） （（m_1+m_2 k_1^2 ） l_1^2+J_1 k_1^2 ）?@???????????????? √（（6m_2 l_2 k_2 （h_1-2（m_1/k_1 +m_2 ）gl_1 ））/（（m_1 k_1+m_2 k_1^2 ） l_1 ）） √（（（m_1+m_2 k_1^2 ） l_1^2+J_1 k_1^2 ） ）?@???????????????? csch（√（（m_2 l_2 k_2 （（m_1+m_2 k_1^2 ） l_1^2+J_1 k_1^2 ））/（2（m_1 k_1+m_2 k_1^2 ） l_1 （m_2 l_2^2+k_2^2 J_2 ） ）） π）sin（√（（m_2 l_2 k_2 g）/（m_2 l_2^2+k_2^2 J_2 ）） t_0 ）。）????? （12）

由式（12）可知，Melnikov函數(shù)為振蕩型，存在簡單零點(diǎn)，則Poincaré映射存在橫截異宿點(diǎn)，由Smale-Birkhoff同宿定理可知，系統(tǒng)存在Smale馬蹄意義下的混沌。綜上可知，當(dāng)ε充分小時，擬Hamilton系統(tǒng)在H≥2（m_1/k_1 +m_2 ）gl_1的能量條件下存在橫截異宿點(diǎn)，系統(tǒng)將呈現(xiàn)混沌的性質(zhì)。擬Hamilton系統(tǒng)在相平面（θ2，θ ˙_2）上沒有不穩(wěn)定流形，也就不存在橫截異宿點(diǎn)，因此只需討論存在不穩(wěn)定流形的相平面（θ1，θ ˙_1）。而對于Hamilton系統(tǒng)，相平面（θ2，θ ˙_2）上的相軌道并非純粹的周期軌道，可能存在不穩(wěn)定流形；不略去Peano余項，F(xiàn)1（θ1，θ ˙_1）也無法參數(shù)化，這些因素或多或少會影響系統(tǒng)的混沌性。但H≥2（m_1/k_1 +m_2 ）gl_1的混沌閾值依然有一定的預(yù)測作用。

前面是在相平面（θ1，θ ˙_1）上求得第1個混沌閾值。若對余弦項cosθ_1進(jìn)行泰勒展開，做變換（θ1，θ ˙_1）?（I，t），重復(fù)上述步驟，也可以得到相似的結(jié)論：當(dāng)ε充分小時，擬Hamilton系統(tǒng)在H≥2/k_2? m_2 gl_2的能量條件下存在橫截同宿點(diǎn)；對于Hamilton系統(tǒng)，在H<2/k_2? m_2 gl_2時也可能處于混沌狀態(tài)。對整個系統(tǒng)而言，無論在哪一相平面上，只要存在橫截同宿點(diǎn)就意味著系統(tǒng)處于混沌狀態(tài)，因此雙擺系統(tǒng)的理論混沌閾值H≥min{2（m_1/k_1 +m_2 ）gl_1，2/k_2? m_2 gl_2 }。

3 數(shù)值仿真驗證混沌條件

本節(jié)將通過數(shù)值方法討論雙擺系統(tǒng)理論混沌閾值的適用條件，并且θ_1、θ_2、θ ˙_1、θ ˙_2特指初始時刻的系統(tǒng)參數(shù)。

3.1 僅考慮動能的混沌條件數(shù)值模擬驗證

觀察Hamilton函數(shù)式與混沌條件H≥min{2（m_1/k_1 +m_2 ）gl_1，2/k_2? m_2 gl_2 }可知，雙擺系統(tǒng)的Hamilton量由（θ1，θ2，θ ˙_1，θ ˙_2，m1，m2，l1，l2）8個物理量構(gòu)成，而混沌條件僅包含（m1，m2，l1，l2），因此可以推測，在（θ1，θ2，m1，m2，l1，l2）相同的參數(shù)條件下，（θ ˙_1，θ ˙_2）越大，越有可能超越混沌閾值。

為驗證上述推測，本節(jié)選取均質(zhì)桿進(jìn)行驗證，均質(zhì)桿可以看作非均質(zhì)桿取k_1=k_2=2，a=b=1/3的特殊情況，可以用來初步驗證所求閾值的正確性。以θ ˙_2為變量，?。é?，θ2，θ ˙_1，m1，m2，l1，l2）=（0，0，0，1，1，1，1）代入混沌條件，可得閾值為θ ˙_2^2≥6g，即|θ ˙_2 |≥7.67，上述條件下，圖2為以θ ˙_2為分岔參數(shù)的分岔圖，可見當(dāng)θ ˙_2=0時，系統(tǒng)始終靜止，因此分岔圖以（0，0）為中心；圖中正、負(fù)半軸近乎對稱，可見速度的方向并不會影響系統(tǒng)的狀態(tài)；在擺的質(zhì)量確定的情況下，速度的大小決定動能的大小，也從側(cè)面反映能量才是決定系統(tǒng)混沌與否的關(guān)鍵。隨著|θ ˙_2 |不斷增大，分岔圖由清晰、連續(xù)的線轉(zhuǎn)變?yōu)槟：?、離散的點(diǎn)，但無法確定具體的混沌閾值。圖3為上述條件下系統(tǒng)的Poincaré截面圖，|θ ˙_2 |=7.6時，Poincaré截面圖依然是封閉的曲線，表現(xiàn)出擬周期的特征；而|θ ˙_2 |=7.7時，Poincaré截面圖已經(jīng)變?yōu)殡s亂無序的點(diǎn)，顯現(xiàn)混沌的特征。圖4為|θ ˙_2 |=7.6和|θ ˙_2 |=7.7時的第二擺上P點(diǎn)運(yùn)動軌跡圖，|θ ˙_2 |=7.6時運(yùn)動軌跡有明顯的區(qū)域范圍，|θ ˙_2 |=7.7時運(yùn)動軌跡則雜亂無序、十分復(fù)雜，說明由Melnikov法獲得的混沌閾值可靠。

3.2 僅考慮勢能的混沌條件數(shù)值模擬驗證

假定雙擺系統(tǒng)無初速度，同時考慮2（m_1/k_1 +m_2 ）gl_1和2/k_2? m_2 gl_2 2個閾值，有

{（H≥2（m_1/k_1 +m_2 ）gl_1@H≥2/k_2? m_2 gl_2 ）┤?{（（k_1+m_1/m_2 ）? （l_1 k_2）/（k_1 l_2 ）≤（1-cosθ_2）/（1+cosθ_1 ），@（k_1+m_1/m_2 ）? （l_1 k_2）/（k_1 l_2 ）≤（1+cosθ_2）/（1-cosθ_1 ） 。）┤? （13）

將式（13）簡化為式（14）

{（（k_1+m）L≤（1-cosθ_2）/（1+cosθ_1 ），@（k_1+m）L≥（1+cosθ_2）/（1-cosθ_1 ） 。）┤??? （14）

由式（14）可知，影響系統(tǒng)混沌的并非擺的總長l_1 、l_2，而是鉸接處到重心的長度L_1 、L_2；代表第一擺質(zhì)心位置的k_1可以單獨(dú)對系統(tǒng)產(chǎn)生影響，而代表第二擺重心位置的k_2必須與l_2組合成L_2來影響系統(tǒng)，可見第一擺重心位置對系統(tǒng)混沌性的影響要大一些。隨角度參數(shù)θ_1 、θ_2的取值不同，會出現(xiàn)兩種假設(shè)情況，一是（1-cosθ_2）/（1+cosθ_1 ）≥（1+cosθ_2）/（1-cosθ_1 ），即cosθ_1+cosθ_2≤0，無論（k_1+m）L取何值，雙擺系統(tǒng)都為混沌運(yùn)動；二是（1-cosθ_2）/（1+cosθ_1 ）<（1+cosθ_2）/（1-cosθ_1 ），即cosθ_1+cosθ_2>0，此時（k_1+m）L在某一范圍內(nèi)，系統(tǒng)存在擬周期運(yùn)動，其余為混沌運(yùn)動。

3.2.1 完全混沌

本節(jié)?。é?，θ2）=（π/3，2π/3）和（θ1，θ2）=（5π/6，2π/3），此時無論（k1+m）L取何值，雙擺系統(tǒng)都應(yīng)為混沌運(yùn)動。同樣選擇均質(zhì)物理雙擺系統(tǒng)進(jìn)行分析，圖5為（θ ˙_1，θ ˙_2，L）=（0，0，2/9）時，隨參數(shù)m變化的最大Lyapunov指數(shù)圖，圖6為在（θ ˙_1，θ ˙_2，m）=（0，0，1）時，隨參數(shù)L變化的最大Lyapunov指數(shù)圖。最大Lyapunov指數(shù)也始終大于0，表明系統(tǒng)與理論推測相同，處于混沌運(yùn)動。兩圖中（b）圖的Lyapunov指數(shù)顯然大于（a）圖的Lyapunov指數(shù)，表明在（b）圖的參數(shù)條件下，系統(tǒng)在相鄰相軌道間發(fā)散得更加嚴(yán)重，系統(tǒng)混沌程度更深。

3.2.2 局部擬周期

當(dāng)（1-cosθ_2）/（1+cosθ_1 ）<（1+cosθ_2）/（1-cosθ_1 ）時，系統(tǒng)應(yīng)存在擬周期運(yùn)動，這里通過數(shù)值仿真，改變第一擺質(zhì)心位置k_1與質(zhì)量比m的取值來驗證理論閾值的正確性。取固定參數(shù)（θ1，θ2，θ ˙_1，θ ˙_2，L，a，b）=（π/3，π/2，0，0，2/9，1/3，1/3），圖7為在不同k1下，以m為分岔參數(shù)的分岔圖。圖7（a）中k1=1.4，圖中m≤1.48處離散點(diǎn)完全發(fā)散，說明系統(tǒng)在（m+k1）≤2.88的區(qū)域做混沌運(yùn)動，略小于理論值；分岔圖在1.48≤m≤7處為十幾條清晰的線相互糾結(jié)在一起，系統(tǒng)做擬周期運(yùn)動；在7≤m≤7.6處清晰的線逐漸轉(zhuǎn)變?yōu)殡x散點(diǎn)，但沒有發(fā)散，可見系統(tǒng)處于擬周期態(tài)到混沌態(tài)的過渡階段；在m=7.6處離散點(diǎn)第一次發(fā)散，可見系統(tǒng)在（m+k1）≥9時開始進(jìn)行混沌運(yùn)動，與理論相符合。圖7（b）中k1=1.7，圖中m≤1.5與m≥7.2處離散點(diǎn)完全發(fā)散，說明系統(tǒng)在（m+k1）≤3.2與（m+k1）≥8.9的區(qū)域做混沌運(yùn)動，混沌區(qū)域略大于理論值。圖7（c）中k1=2，此即均質(zhì)雙擺，圖中m≤1.55與m≥5.67處離散點(diǎn)完全發(fā)散，說明系統(tǒng)在（m+k1）≤3.55與（m+k1）≥7.67處做混沌運(yùn)動，混沌區(qū)域雖略大于理論值，但在理論區(qū)域確實是做混沌運(yùn)動。圖7（d）中k1=2.2，系統(tǒng)在（m+k1）≤3.8時做混沌運(yùn)動，在（m+k1）≥7.2時出現(xiàn)離散點(diǎn)并開始向混沌態(tài)過渡，直到（m+k1）≥7.5時開始做混沌運(yùn)動。

事實表明，當(dāng)（1-cosθ_2）/（1+cosθ_1 ）<（1+cosθ_2）/（1-cosθ_1 ）時，系統(tǒng)存在擬周期態(tài)，且當(dāng)擺的質(zhì)心越接近擺尾，即k_1越小時，系統(tǒng)的擬周期運(yùn)動區(qū)域越大。

4 模型局限性分析

4.1 k_1≈（1-cosθ_2）/（1+cosθ_1 ）L或k_1>（1-cosθ_2）/（1+cosθ_1 ）L

圖8為取參數(shù)（θ1，θ2，θ ˙_1，θ ˙_2，L，a，b）=（π/3，π/2，0，0，2/9，1/3，1/3）時以m為分岔參數(shù)的分岔圖，此時（1-cosθ_2）/（1+cosθ_1 ）L=3。圖8（a）中k1=2.8，按照理論計算，系統(tǒng)應(yīng)在m≤0.2與m≥6.2的區(qū)域內(nèi)做混沌運(yùn)動，圖中0.06≤m≤0.2時系統(tǒng)做擬周期運(yùn)動，與理論預(yù)測不符，但依然滿足m≥6.2的混沌條件。圖8（b）中k1=3.2，圖中5.8≤m≤6.93時系統(tǒng)做擬周期運(yùn)動，與系統(tǒng)應(yīng)在m≥5.8時做混沌運(yùn)動不相符，說明第一擺質(zhì)心位置的k_1在數(shù)值上接近（1-cosθ_2）/（1+cosθ_1 ）L時，理論閾值只能正確預(yù)測（k_1+m）L≥（1+cosθ_2）/（1-cosθ_1 ）的部分，當(dāng)k_1超越（1-cosθ_2）/（1+cosθ_1 ）L時，理論閾值將完全不適用。

4.2 m_1?m_2

圖9為取參數(shù)（θ1，θ2，θ ˙_1，θ ˙_2，k1，a，b）=（π/3，π/2，0，0，2，1/3，1/3）時以m為分岔參數(shù)的分岔圖，該圖表明，若m較大時，在l_1?l_2的情況下，第一擺的能量占系統(tǒng)總能量的絕大部分，第二擺對系統(tǒng)的影響很?。欢趌_2?l_1時，兩擺間難以產(chǎn)生較大的影響，此時可以將系統(tǒng)看作為1個單擺系統(tǒng)加上1個第二擺產(chǎn)生的微小擾動。因此，在這2種條件下系統(tǒng)在參數(shù)m較大的區(qū)域進(jìn)行擬周期運(yùn)動，而非混沌運(yùn)動。

4.3 J_2?J_1或J_1?J_2

圖10為取參數(shù)（θ1，θ2，θ ˙_1，θ ˙_2，L，k1）=（π/3，π/2，0，0，2/9，1.7）時以m為分岔參數(shù)的分岔圖，在此參數(shù)條件下系統(tǒng)將在m≤1.3與m≥7.3的區(qū)域做混沌運(yùn)動，此時J_1/J_2 =4am/81b。

當(dāng)J_2?J_1即b?a時，如圖10（a）（b）中分別?。╝，b）=（0.1，0.4）和（a，b）=（0.1，0.8），系統(tǒng)在0≤m≤10內(nèi)由混沌態(tài)轉(zhuǎn)為擬周期態(tài)后，沒有再出現(xiàn)明確的混沌區(qū)閾值，表明當(dāng)?shù)诙[轉(zhuǎn)動慣量遠(yuǎn)大于第一擺轉(zhuǎn)動慣量時，隨著雙擺質(zhì)量比m的增加，系統(tǒng)將出現(xiàn)趨于穩(wěn)定的擬周期運(yùn)動，其原因是在雙擺運(yùn)動過程中，轉(zhuǎn)動慣量會影響系統(tǒng)的動能，根據(jù)機(jī)械能守恒定理知，第二擺的轉(zhuǎn)動慣量越大，速度就會越小，此時可以把第二擺的運(yùn)動看做在單擺系統(tǒng)上加一個微小的擾動，因此系統(tǒng)會趨于擬周期運(yùn)動。

當(dāng)J_2≈J_1時，如圖10（c）中?。╝，b）=（0.3，0.3）時，系統(tǒng)在0≤m≤1.63時為混沌運(yùn)動，在1.63≤m≤5.76時為擬周期運(yùn)動，在5.76≤m≤10時為混沌運(yùn)動，雖然實際混沌區(qū)域大于理論預(yù)測，但在理論預(yù)測的區(qū)域確實做混沌運(yùn)動。如圖10（d）中取（a，b）=（0.6，0.3）時，系統(tǒng)在0≤m≤1.3時為混沌運(yùn)動，在1.3≤m≤6.73時為擬周期運(yùn)動，在6.73≤m≤10時為混沌運(yùn)動，此時理論閾值具有很好的適用性。

當(dāng)J_1?J_2時，如圖10（e）（f）中分別?。╝，b）=（0.8，0.1）和（a，b）=（0.9，0.1），此時系統(tǒng)在0≤m≤10內(nèi)作由混沌態(tài)轉(zhuǎn)變?yōu)閿M周期態(tài)又轉(zhuǎn)變?yōu)榛煦鐟B(tài)的循環(huán)往復(fù)的運(yùn)動。這說明當(dāng)?shù)谝粩[的轉(zhuǎn)動慣量遠(yuǎn)大于第二擺的轉(zhuǎn)動慣量時，系統(tǒng)的運(yùn)動將受到兩擺質(zhì)量比的影響而出現(xiàn)極其復(fù)雜的運(yùn)動。對比六次仿真結(jié)果可知，系統(tǒng)在運(yùn)動過程中，轉(zhuǎn)動慣量會對系統(tǒng)的動能產(chǎn)生影響從而改變系統(tǒng)混沌條件，當(dāng)兩擺的轉(zhuǎn)動慣量差距過大時，Melnikov法得出的能量閾值不再適用。

5 結(jié)? 論

建立了非均質(zhì)物理雙擺運(yùn)動的擬Hamilton模型，根據(jù)雙自由度的Melnikov法提出了擬Hamilton系統(tǒng)發(fā)生混沌的能量閾值，得到了一般情況下系統(tǒng)發(fā)生混沌的條件。通過改變雙擺模型中的質(zhì)心位置、擺長比以及轉(zhuǎn)動慣量來體現(xiàn)雙擺非均質(zhì)的特點(diǎn)，再根據(jù)不同參數(shù)下系統(tǒng)的分岔圖、Poincaré截面圖對理論閾值進(jìn)行驗證，給出了模型的適用條件，并完成了局限性分析?；诒狙芯浚傻贸鲆韵陆Y(jié)論：

1）分析非均質(zhì)物理雙擺混沌的理論閾值，發(fā)現(xiàn)影響系統(tǒng)混沌的并非雙擺的擺長l1、l2，而是雙擺質(zhì)心到鉸接處的距離L1、L2，并且在雙擺質(zhì)量和長度確定的情況下，L1、L2越大，系統(tǒng)發(fā)生混沌所需要的能量就越大。

2）分析雙擺質(zhì)心位置對系統(tǒng)混沌特性的影響，發(fā)現(xiàn)第一擺的質(zhì)心位置k1可以直接影響系統(tǒng)混沌特性，而第二擺的質(zhì)心位置k2則不能單獨(dú)對系統(tǒng)的混沌特性產(chǎn)生影響。利用數(shù)值仿真改變k1的大小，從而進(jìn)一步得到了第一擺質(zhì)心位置對理論閾值的適用性，發(fā)現(xiàn)當(dāng)k1在數(shù)值上接近（1-cosθ_2）/（1+cosθ_1 ）L時，理論閾值的適用性變差，當(dāng)k1大于（1-cosθ_2）/（1+cosθ_1 ）L時，理論閾值將完全不適用。

3）分析雙擺轉(zhuǎn)動慣量對系統(tǒng)混沌特性的影響，發(fā)現(xiàn)雙擺通過轉(zhuǎn)動慣量的改變使得運(yùn)動時動能發(fā)生改變，進(jìn)而影響系統(tǒng)的混沌閾值，并且知道當(dāng)?shù)诙[的轉(zhuǎn)動慣量遠(yuǎn)大于第一擺的轉(zhuǎn)動慣量時，系統(tǒng)將趨于擬周期運(yùn)動；當(dāng)?shù)谝粩[的轉(zhuǎn)動慣量遠(yuǎn)大于第二擺時，系統(tǒng)作混沌態(tài)、擬周期態(tài)循環(huán)往復(fù)的復(fù)雜運(yùn)動；當(dāng)雙擺的轉(zhuǎn)動慣量相似時，系統(tǒng)能很好地在理論閾值區(qū)域內(nèi)做混沌運(yùn)動。表明在極限條件下，理論閾值存在一定的局限。

4）非均質(zhì)雙擺系統(tǒng)的混沌特性研究對類似系統(tǒng)的混沌抑制有著重要的啟示作用。具體來說，可以通過研究雙擺系統(tǒng)的混沌控制方法來抑制類似系統(tǒng)的混沌特性，實現(xiàn)系統(tǒng)的穩(wěn)定控制。例如，對于機(jī)器人肩部和腿部的運(yùn)動，通過控制雙擺系統(tǒng)的擺臂長度和擺錘重量比等參數(shù)，可以實現(xiàn)關(guān)節(jié)運(yùn)動的控制和規(guī)劃。另外，在航空航天領(lǐng)域中，太陽能帆板的展開涉及到桿件和軸承的不確定性，容易導(dǎo)致非線性振動和失穩(wěn)現(xiàn)象的發(fā)生，因此可以通過研究雙擺混沌控制方法來控制太陽能帆板的展開過程，實現(xiàn)太陽能帆板的準(zhǔn)確展開?？傊ㄟ^研究非均質(zhì)物理雙擺系統(tǒng)的混沌特性，在工程實際中可以針對類似系統(tǒng)的混沌特性進(jìn)行控制和抑制，提高系統(tǒng)的穩(wěn)定性和控制精度。

參考文獻(xiàn)

［1］? Petró B， Kiss B， Kiss R M. Analyzing human balance recovery action using calculated torques of a double pendulum model[J]. Materials Today： Proceedings， 2019， 12： 431-439.

［2］? Betzler N， Monk S， Wallace E， et al. From the double pendulum model to full-body simulation： evolution of golf swing modeling[J]. Sports Technology， 2008， 1（4/5）： 175-188.

［3］? Youn S H. Double pendulum model for a tennis stroke including a collision process[J]. Journal of the Korean Physical Society， 2015， 67（7）： 1110-1117.

［4］? Verduzco F， Alvarez J. Homoclinic chaos in 2-DOF robot manipulators driven by PD controllers[J]. Nonlinear Dynamics， 2000， 21（2）： 157-171.

［5］? 張紅巧， 田瑞蘭， 陳恩利， 等. 參數(shù)激勵下均質(zhì)桿狀雙擺的周期穩(wěn)定振動[J]. 振動與沖擊， 2020， 39（16）： 231-235.

Zhang H Q， Tian R L， Chen E L， et al. Periodically stable vibration of homogeneousrod-shaped double pendulum under parametric excitation[J]. Journal of Vibration and Shock， 2020， 39（16）： 231-235.（in Chinese）

［6］? Shi H T， Yao F X， Yuan Z， et al. Research on nonlinear coupled tracking controller for double pendulum gantry cranes with load hoisting/lowering[J]. Nonlinear Dynamics， 2022， 108（1）： 223-238.

［7］? Sun N， Wu Y M， Chen H， et al. An energy-optimal solution for transportation control of cranes with double pendulum dynamics： design and experiments[J]. Mechanical Systems and Signal Processing， 2018， 102： 87-101.

［8］? Zhang M H， Ma X， Chai H， et al. A novel online motion planning method for double-pendulum overhead cranes[J]. Nonlinear Dynamics， 2016， 85（2）： 1079-1090.

［9］? Kumar R， Gupta S， Ali S F. Energy harvesting from chaos in base excited double pendulum[J]. Mechanical Systems and Signal Processing， 2019， 124： 49-64.

［10］? 王觀， 胡華， 伍康， 等. 基于兩級擺桿結(jié)構(gòu)的超低頻垂直隔振系統(tǒng)[J]. 物理學(xué)報， 2016， 65（20）： 200702.

Wang G， Hu H， Wu K， et al. Ultra-low-frequency vertical vibration isolator based on a two-stage beam structure[J]. Acta Physica Sinica， 2016， 65（20）： 200702.（in Chinese）

［11］? 范洪義， 吳澤. 介觀電路中量子糾纏的經(jīng)典對應(yīng)[J]. 物理學(xué)報， 2022， 71（1）： 010302.

Fan H Y， Wu Z. Classical correspondence of quantum entanglement in mesoscopic circuit[J]. Acta Physica Sinica， 2022， 71（1）： 010302.（in Chinese）

［12］? 陳漢軍， 黃東衛(wèi)， 蘇永福. 雙物理擺的混沌分析[J]. 數(shù)學(xué)的實踐與認(rèn)識， 2007， 37（15）： 49-53.

Chen H J， Huang D W， Su Y F. Chaos analysis of the double-body pendulum[J]. Mathematics in Practice and Theory， 2007， 37（15）： 49-53.（in Chinese）

［13］? Ford J， Ilg M. Eigenfunctions， eigenvalues， and time evolution of finite， bounded， undriven， quantum systems are not chaotic[J]. Physical Review A， 1992， 45（9）： 6165-6173.

［14］? Stachowiak T， Okada T. A numerical analysis of chaos in the double pendulum[J]. Chaos， Solitons & Fractals， 2006， 29（2）： 417-422.

［15］? 趙武， 張鴻斌， 孫超凡， 等. 受垂直激勵和水平約束的單擺系統(tǒng)亞諧共振分岔與混沌[J]. 物理學(xué)報， 2021， 70（24）： 240202.

Zhao W， Zhang H B， Sun C F， et al. Subharmonic resonance bifurcation and chaos of simple pendulum system with vertical excitation and horizontal constraint[J]. Acta Physica Sinica， 2021， 70（24）： 240202.（in Chinese）

［16］? 劉丁楊， 蹇開林， 張亮. 基于Melnikov法的雙擺系統(tǒng)混沌特性研究[J]. 振動與沖擊， 2022， 41（14）： 92-98.

Liu D Y， Jian K L， Zhang L. Chaotic characteristics of a double pendulum system based on the Melnikov method[J]. Journal of Vibration and Shock， 2022， 41（14）： 92-98.（in Chinese）

［17］? Han N， Lu P P. Nonlinear dynamics of a classical rotating pendulum system with multiple excitations[J]. Chinese Physics B， 2020， 29（11）： 110502.

［18］? Maiti S， Roy J， Mallik A K， et al. Nonlinear dynamics of a rotating double pendulum[J]. Physics Letters A， 2016， 380（3）： 408-412.

［19］? Bi Q S， Chen Y S. Bifurcation analysis of a double pendulum with internal resonance[J]. Applied Mathematics and Mechanics， 2000， 21（3）： 255-264.

［20］? Martynyuk A A， Nikitina N V. The theory of motion of a double mathematical pendulum[J]. International Applied Mechanics， 2000， 36（9）： 1252-1258.

［21］? Calv?o A M， Penna T P. The double pendulum： a numerical study[J]. European Journal of Physics， 2015， 36（4）： 045018.

［22］? Oiwa S， Yajima T. Jacobi stability analysis and chaotic behavior of nonlinear double pendulum[J]. International Journal of Geometric Methods in Modern Physics， 2017， 14（12）： 1750176.

［23］? Kovacic I， Zukovic M， Radomirovic D. Normal modes of a double pendulum at low energy levels[J]. Nonlinear Dynamics， 2020， 99（3）： 1893-1908.

［24］? Dudkowski D， Wojewoda J， Czo?czyński K， et al. Is it really chaos？ The complexity of transient dynamics of double pendula[J]. Nonlinear Dynamics， 2020， 102（2）： 759-770.

［25］? 李繼彬， 陳鳳娟. 混沌、Melnikov方法及新發(fā)展[M]. 北京： 科學(xué)出版社， 2012.

Li J B， Chen F J. Chaos， Melnikov method and its new development[M]. Beijing： Science Press， 2012.（in Chinese）

［26］? Holmes P J， Marsden J E. Horseshoes in perturbations of Hamiltonian systems with two degrees of freedom[J]. Communications in Mathematical Physics， 1982， 82（4）： 523-544.

（編輯? 鄭潔）