于 昊，喻 勇

(1. 西南交通大學(xué) 力學(xué)與航空航天學(xué)院，四川 成都 611756;2. 西安交通大學(xué) 應(yīng)用力學(xué)與結(jié)構(gòu)安全重點實驗室，四川 成都 611756)

混沌現(xiàn)象是20世紀一項重大的科學(xué)發(fā)現(xiàn)，至今仍然具有很高的研究意義.混沌即在確定性系統(tǒng)中表現(xiàn)出來的隨機的不規(guī)則運動，一個系統(tǒng)雖然可以用確定性理論來描述，但實際行為卻表現(xiàn)出不確定性、不可重復(fù)、不可預(yù)測等性質(zhì)，這就是混沌現(xiàn)象[1].非線性混沌科學(xué)不僅在理論上具有重大意義，其在生態(tài)、醫(yī)療、金融以及決策等問題上也具有很大的價值，混沌對于非常復(fù)雜、系統(tǒng)性疾病的研究很有幫助，也可以進行氣象特征的分析等[2].此外，混沌還可用于信息加密，通過混沌系統(tǒng)得到的密碼，具有了混沌的特征，使人們一時不能破解，也無法預(yù)測出密碼的信息[3].混沌擺對于非線性系統(tǒng)的力學(xué)行為演示更為直接，能把混沌中所蘊含的確定性和不確定性展示出來，對于混沌理論的理解很有幫助[4].混沌擺還是力學(xué)教學(xué)中的熱點.目前來說，雖然在力學(xué)教學(xué)過程當(dāng)中并沒有規(guī)定進行混沌擺理論的介紹，但是向?qū)W生適當(dāng)進行相關(guān)知識的介紹不僅會打破學(xué)生確定性的概念，也會增加其對力學(xué)學(xué)習(xí)的興趣[5].

唐有綺等[5]建立四自由度無阻尼混沌擺實驗?zāi)Ｐ?，進而建立其動力學(xué)方程，通過求解微分方程獲得系統(tǒng)的擺角時程圖，揭示了混沌現(xiàn)象.朱桂萍等[6]分別研究了混沌擺在保守系統(tǒng)和耗散系統(tǒng)下的動力學(xué)性質(zhì)，得到了兩種系統(tǒng)情形下的相圖和角速度時序圖，并分析了系統(tǒng)參數(shù)對系統(tǒng)的影響.孟勇[7]采用Maple軟件對大角度單擺、雙擺、傅科擺的擺動進行了模擬，分析了每個擺的運動情況和物理規(guī)律.李明達等[8]建立混沌擺模型，得到系統(tǒng)的相圖和分叉圖，分析系統(tǒng)對外界參數(shù)的敏感性特征.

石根華[9-10]提出的非連續(xù)變形分析方法(Discontinuous Deformation Analysis，DDA)平行于有限單元法，其在分析非連續(xù)變形問題中，可以說是最有效的方法之一.該方法以塊體系統(tǒng)的位移作為未知量，從而建立系統(tǒng)的平衡方程，進行求解.由于其計算比較高效，目前已經(jīng)被廣泛應(yīng)用到模擬落石、滑坡、隧道等方面.

目前大部分對于混沌擺運動的揭示都是通過建立系統(tǒng)的拉格朗日方程，然后對方程求解，得到角度的解析解，然后得到系統(tǒng)的相圖、頻譜圖或者時程圖來進行分析.雖然也可采用數(shù)值軟件實現(xiàn)方程結(jié)果的可視化，但是前期對于方程的求解略微枯燥，對于不同的模型，可能求解出的結(jié)果非常復(fù)雜.本文采用非連續(xù)變形分析方法，利用Matlab編寫DDA的計算程序進行兩種擺運動的模擬，不需要人工求解復(fù)雜的微分方程，并且可以得到運動的動畫，使顯示更加直觀，更容易理解.

但是在DDA程序運行過程中，需要對程序的具體參數(shù)進行合理取值，如彈簧剛度和時間步長等，如果取值不當(dāng)，將會導(dǎo)致程序不能運行或運行出錯.本文對于程序中具體參數(shù)按照文獻[11-13]所討論進行取值.

1 DDA方程以及模型的建立

本節(jié)介紹DDA方程的構(gòu)造方法以及兩種模型的建立.

1.1 DDA方程的建立

DDA方程的建立關(guān)鍵是確定方程的子矩陣以及力矩陣的組成及其表達式，本節(jié)介紹DDA方程的建立依據(jù)以及子矩陣的確定.

1.1.1 位移模式

在本文中，只考慮二維塊體，并且塊體都為剛體，因此可以將原始DDA中的位移模式進行較大簡化.

在原始DDA中塊體的位移矩陣以及子矩陣都為6階矩陣.在考慮二維剛體情況下，矩陣都將變?yōu)?階矩陣，這將使得形式更加簡單，計算更加方便.

塊體內(nèi)任一點(x,y)的位移U=(u,v)T可通過位移函數(shù)確定，即

U=T·D

(1)

1.1.2 平衡方程組

多個塊體之間是依靠約束相互連接的，因此塊體間形成一聯(lián)立方程組，以本文中的四塊體系統(tǒng)為例

(2)

對于二維剛性塊體，系數(shù)矩陣kij(i,j=1,2,3,4)為3×3矩陣，Di,Fi(i=1,2,3,4)為3×1矩陣.系數(shù)矩陣kij的具體數(shù)值由塊體和塊體間的連接方式確定.

1.1.3 方程組子矩陣的確定

首先由于系統(tǒng)中每個塊體受到重力影響，由此形成體積荷載：

(3)

其中S為塊體面積，fy為體積荷載.該子矩陣加入到總體方程中的[Fi]中去.其次，塊體在運動過程中會受到慣性力的影響，由此產(chǎn)生兩個子矩陣：

(4)

(5)

式(4)加入到總體方程[kii]中去，式(5)加入到總體方程[Fi]中，其中，M為塊體質(zhì)量，Δ為時間步長，Vi(0)為該時間步的初始速度，為上一時間步末的速度，具體計算公式為：

(6)

最后，需要確定主擺和副擺之間的連接方式.在原始DDA中，對塊體可以添加點位移、點荷載和錨桿連接等方式.主擺與復(fù)擺間的連接方式會對模擬結(jié)果產(chǎn)生重要的影響，不同的連接方式會產(chǎn)生不同子矩陣，通過多次實驗，本文繼續(xù)采用原始DDA中所用到的錨桿連接，錨桿連接方式會產(chǎn)生4個子矩陣.

采用錨桿將主擺和復(fù)擺短邊中點相連接，由此會產(chǎn)生4個子矩陣：

(7)

(8)

(9)

(10)

式(7)加入到總體方程的子矩陣[Kii]中去，式(8)加入到總體方程的子矩陣[Kij]中去，式(9)加入到總體方程的子矩陣[Kji]中去，式(10)加入到總體方程的子矩陣[Kjj]中去.其中s是錨桿的剛度，l是錨桿的長度，i、j表示兩個塊體.[Ei]和[Gj]具體計算表達式：

(11)

(12)

1.2 建立模型

分別建立2自由度和4自由度混沌擺模型進行研究，本節(jié)簡要介紹模型構(gòu)成以及具體參數(shù).

1.2.1 兩自由度混沌擺模型的建立

建立如圖1所示的2自由度無阻尼混沌擺模型，由一個T型主擺以及與其相連的一個副擺組成.轉(zhuǎn)動時，通過給T型主擺一角位移，使其繞中間一固定點O旋轉(zhuǎn).T型主擺可看作由3根等大的擺組成，其中每個擺的質(zhì)量為M，副擺的質(zhì)量為M，每根桿的長度均為L.

圖1 兩自由度混沌擺模型

1.2.2 四自由度混沌擺模型的建立

建立4自由度無阻尼混沌擺，如圖2所示.該裝置由一個T型主擺以及所連接的3個副擺所組成.該T型主擺繞中間一固定點O旋轉(zhuǎn)，并通過錨桿與3個副擺相連，這樣四根桿組成四自由度無阻尼自由振動系統(tǒng).當(dāng)進行實驗時，通過給T型主擺一角位移，使整個系統(tǒng)進行運動.T型主擺可看作由3根等大的擺組成，其每個擺的質(zhì)量為M，副擺的質(zhì)量為M，每根桿的長度均為L.

圖2 四自由度混沌擺模型

1.3 理論分析

本節(jié)介紹2自由度和4自由度混沌擺理論結(jié)果的推導(dǎo)以及計算.

1.3.1 兩自由度混沌擺理論結(jié)果

根據(jù)文獻[2]，對兩自由度混沌擺運動微分方程進行如下推導(dǎo).

系統(tǒng)動能為

(13)

系統(tǒng)勢能為

(14)

其中，m1、m2分別為主擺和副擺的質(zhì)量，l1、l2分別為主擺、副擺的長度.

在本文中，不考慮阻力的影響，因此將動能、勢能代入到保守系拉格朗日方程中，對拉格朗日方程采用4階Runge-Kutta法進行求解，利用Matlab編寫了相關(guān)計算程序，其中使用了ode45函數(shù)，得到擺的運動微分方程組：

4gm2sinθ1-3gm2cos(θ1-θ2)sinθ2]/

{l1[4m1+4m2-3m2cos2(θ1-θ2)]}，

2gm1sinθ2+2gm2cos(θ1-θ2)sinθ1+

{l2[4m1+4m2-3m2cos2(θ1-θ2)]}

(15)

1.3.2 四自由度混沌擺理論結(jié)果

根據(jù)文獻[2]，因文獻[2]公式輸入錯誤，對其進行更改，得到正確的4自由度混沌擺運動微分方程推導(dǎo)如下.

系統(tǒng)動能為

(16)

系統(tǒng)勢能為

m2l2(cosθ2+cosθ3+cosθ4)]

(17)

其中，m1、m2分別為主擺和副擺的質(zhì)量，l1、l2分別為主擺、副擺的長度.

得到擺的運動微分方程組：

3m2sin(θ1-2θ2)-(4m1+5m2)sinθ1]-

3l1m2[sin 2(θ1-θ2)-sin 2(θ1-θ3)-

{8l1m1+3l1m2[5-cos 2(θ1-θ2)+

cos 2(θ1-θ3)+cos 2(θ1-θ4)]}，

(18)

2 數(shù)值仿真

在本文中，系統(tǒng)無阻尼且塊體之間的連接都為光滑連接.塊體具體參數(shù)如下，質(zhì)量M=2×103kg，長度L=10 m，重力加速度g=9.8 m/s2.DDA程序參數(shù)如下：錨桿的剛度s取為5×1013N/m，錨桿的長度l取為固定值0.01 m ,時間步長為0.01 s.通過給T型擺一初始角位移使整個系統(tǒng)運動，初始條件如下：

θ′1=0，θ′2=0，θ′3=0，θ′4=0.

2.1 兩自由度混沌擺仿真

各桿無初速度釋放，可以得到兩自由度混沌擺各桿角度隨時間變化的具體數(shù)值，并將該數(shù)值與理論解進行對比，得到如圖3、圖4所示的角度對比圖，可發(fā)現(xiàn)DDA計算的結(jié)果與理論結(jié)果相差不大，吻合性較好.

圖3 在初始條件θ1下運動10 s的時程圖

圖4 在初始條件θ1下運動20 s的時程圖

在運行20 s情況下，選取最下方塊體左下角點，可以畫出其軌跡曲線如圖5所示，可以發(fā)現(xiàn)其運動為無規(guī)律的混沌現(xiàn)象.

圖5 兩自由度混沌擺一點軌跡圖

2.2 四自由度混沌擺仿真

各桿無初速度釋放，可以得到4自由度混沌擺各桿角度隨時間變化的具體數(shù)值，并將該數(shù)值與理論解進行對比，得到如圖6、圖7所示的角度對比圖，可發(fā)現(xiàn)DDA計算的結(jié)果與理論微分方程計算的結(jié)果相差不大，吻合性較好.

圖6 在初始條件θ1下運動10 s的時程圖

圖7 在初始條件θ1下運動20 s的時程圖

在運行20 s情況下，選取最下方塊體左下角點，可以畫出其軌跡曲線如圖8所示，可以發(fā)現(xiàn)其運動為無規(guī)律的混沌現(xiàn)象.

圖8 四自由度混沌擺一點軌跡圖

從圖8中可以發(fā)現(xiàn)，當(dāng)初始條件較大時，系統(tǒng)即呈現(xiàn)出非周期運動的特點，桿的運動將變得不規(guī)律，數(shù)值模擬出的結(jié)果與理論解差別很小，說明了該DDA程序的適用性，也說明了對于塊體間連接方式選擇錨桿連接的正確性.

3 結(jié)論

本文不同于以往用微分方程直接研究混沌擺的運動，采用非連續(xù)變形分析方法，對2自由度和4自由度混沌擺分別進行模擬，得到以下結(jié)論：

1) 運用Matlab軟件編寫DDA程序，分別將模型運行10 s和20 s，并與理論結(jié)果進行對比，發(fā)現(xiàn)角度變化曲線吻合較好，說明了該DDA程序的適用性；

2) 運用DDA方法所給出的錨桿連接方式，可以近似代替混沌擺中主擺與復(fù)擺之間的鉚釘連接，并且結(jié)果相差不大；

3) 隨著時間的延長，DDA計算的曲線與理論解得到的曲線仍會存在一定誤差，是因為在每一時間步內(nèi)加速度為常數(shù)，因此隨著計算時間的增加，引起誤差的積累，從而導(dǎo)致曲線出現(xiàn)偏差. 此外，對于連接塊體的錨桿長度可能也會對結(jié)果產(chǎn)生影響.