于 昊，喻 勇

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

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

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

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

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

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

1 DDA方程以及模型的建立

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

1.1 DDA方程的建立

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

1.1.1 位移模式

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

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

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

U=T·D

(1)

1.1.2 平衡方程組

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

(2)

對(duì)于二維剛性塊體，系數(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)中每個(gè)塊體受到重力影響，由此形成體積荷載：

(3)

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

(4)

(5)

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

(6)

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

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

(7)

(8)

(9)

(10)

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

(11)

(12)

1.2 建立模型

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

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

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

圖1 兩自由度混沌擺模型

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

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

圖2 四自由度混沌擺模型

1.3 理論分析

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

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

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

系統(tǒng)動(dòng)能為

(13)

系統(tǒng)勢(shì)能為

(14)

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

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

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ù)文獻(xiàn)[2]，因文獻(xiàn)[2]公式輸入錯(cuò)誤，對(duì)其進(jìn)行更改，得到正確的4自由度混沌擺運(yùn)動(dòng)微分方程推導(dǎo)如下.

系統(tǒng)動(dòng)能為

(16)

系統(tǒng)勢(shì)能為

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

(17)

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

得到擺的運(yùn)動(dòng)微分方程組：

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，長(zhǎng)度L=10 m，重力加速度g=9.8 m/s2.DDA程序參數(shù)如下：錨桿的剛度s取為5×1013N/m，錨桿的長(zhǎng)度l取為固定值0.01 m ,時(shí)間步長(zhǎng)為0.01 s.通過給T型擺一初始角位移使整個(gè)系統(tǒng)運(yùn)動(dòng)，初始條件如下：

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

2.1 兩自由度混沌擺仿真

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

圖3 在初始條件θ1下運(yùn)動(dòng)10 s的時(shí)程圖

圖4 在初始條件θ1下運(yùn)動(dòng)20 s的時(shí)程圖

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

圖5 兩自由度混沌擺一點(diǎn)軌跡圖

2.2 四自由度混沌擺仿真

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

圖6 在初始條件θ1下運(yùn)動(dòng)10 s的時(shí)程圖

圖7 在初始條件θ1下運(yùn)動(dòng)20 s的時(shí)程圖

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

圖8 四自由度混沌擺一點(diǎn)軌跡圖

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

3 結(jié)論

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

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

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

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