趙文麗 張 紅 厲桂華 高 峰

（山東農(nóng)業(yè)大學(xué)信息科學(xué)與工程學(xué)院 山東 泰安 271018）

1 引言

彈簧擺是由輕質(zhì)彈簧和懸掛的小球構(gòu)成，在鉛直平面內(nèi)運(yùn)動(dòng)的彈簧可以看成是一個(gè)雙自由度振動(dòng)系統(tǒng)，是典型的非線性運(yùn)動(dòng)之一.關(guān)于彈簧擺的混沌行為引起了人們的廣泛關(guān)注［1～6］.這些報(bào)道從理論求解、數(shù)值分析和計(jì)算機(jī)模擬等方面對(duì)彈簧擺的運(yùn)動(dòng)作了一定的研究.彈簧擺的運(yùn)動(dòng)可以分解為振動(dòng)和擺動(dòng)兩種模式的合成，對(duì)應(yīng)于這兩種運(yùn)動(dòng)模式有兩個(gè)固有頻率，當(dāng)這兩個(gè)頻率滿足一定條件時(shí)，兩個(gè)振動(dòng)模態(tài)強(qiáng)烈地耦合，一種振動(dòng)會(huì)激發(fā)另外一種振動(dòng)，這種現(xiàn)象就是彈簧擺的內(nèi)共振現(xiàn)象.本文從系統(tǒng)的動(dòng)力學(xué)方程出發(fā)，利用Mathematica數(shù)值計(jì)算和模擬功能描繪了彈簧擺運(yùn)動(dòng)的軌跡圖、振動(dòng)曲線和能量曲線，從而使得彈簧擺的運(yùn)動(dòng)更加直觀.

2 彈簧擺的動(dòng)力學(xué)方程

如圖1所示，彈簧擺是由質(zhì)量忽略不計(jì)的輕質(zhì)彈簧與懸掛著的質(zhì)量為m的小球構(gòu)成.彈簧勁度系數(shù)為κ，原長(zhǎng)為l0，彈簧的瞬時(shí)長(zhǎng)度為l，小球的直徑可以忽略不計(jì).取O點(diǎn)為擺球的平衡位置，沿水平方向?yàn)閤軸，豎直方向?yàn)閦軸，建立直角坐標(biāo)系.系統(tǒng)的勢(shì)能為

圖1 彈簧擺

根據(jù)式（1）可以得到系統(tǒng)的動(dòng)力學(xué)方程為

理論上來(lái)講，求解式（4）和式（5）構(gòu)成的方程組，就可以分析彈簧擺的運(yùn)動(dòng)行為，而實(shí)際上，這兩個(gè)方程沒(méi)有嚴(yán)格的解析解，研究該體系的動(dòng)力學(xué)行為需要對(duì)這兩個(gè)方程進(jìn)行數(shù)值計(jì)算.

3 彈簧擺內(nèi)共振的Mathematica仿真

數(shù)學(xué)軟件Mathematica具有強(qiáng)大的數(shù)值計(jì)算功能［7］，利用Mathematica的NDSolve命令可以在不作任何近似的情況下直接求得彈簧擺的動(dòng)力學(xué)微分方程（4）、（5）的數(shù)值解.考慮到小球運(yùn)動(dòng)速度不大，可以認(rèn)為過(guò)程中阻力與速度成正比，因此對(duì)方程（4）、（5）修正如下

其中η是阻力系數(shù)，這主要考慮的是空氣摩擦阻力.為了檢驗(yàn)數(shù)值計(jì)算的精確性，計(jì)算中直接選取文獻(xiàn)［6］中的實(shí)驗(yàn)數(shù)據(jù)，m=4.86×10-2kg，κ=2.30 N／m，z0=0.82 m，此時(shí)，ωp=3.46 rad／s，ωs=6.89 rad／s，滿足內(nèi)共振條件.兩種不同的初始條件為

初始條件1：x（0）=14.1 cm，z（0）=10.9 cm（對(duì)應(yīng)于小擺幅振動(dòng)）

初始條件2：x（0）=19.6 cm，z（0）=5.74 cm（對(duì)應(yīng)于大擺幅振動(dòng)）

分別在兩種不同的初始條件下考察彈簧擺運(yùn)動(dòng)軌跡、兩個(gè)方向的位移和能量隨時(shí)間的變化.

圖2為彈簧擺的運(yùn)動(dòng)軌跡，結(jié)果顯示初始條件不同，軌跡所形成的輪廓不同，彈簧擺的運(yùn)動(dòng)范圍不同，表明彈簧擺的實(shí)際運(yùn)動(dòng)軌跡取決于初始條件.但是無(wú)論對(duì)應(yīng)哪一種初始條件，彈簧擺的伸縮擺動(dòng)都只能在一定的空間范圍內(nèi)，其活動(dòng)范圍都是有界的.

圖2 彈簧擺在不同初始條件下的運(yùn)動(dòng)軌跡

圖3和圖4是兩種初始條件下兩個(gè)模態(tài)的振動(dòng)曲線，從圖中可以看出，彈簧擺無(wú)論是做小擺幅擺動(dòng)還是做大擺幅擺動(dòng)，x模態(tài)和z模態(tài)都有倍頻關(guān)系的內(nèi)共振現(xiàn)象，但是唯有小擺幅的振動(dòng)，才有強(qiáng)烈耦合的情形.由圖3（b）可以看出，當(dāng)在z方向的振幅趨于零時(shí)，z模態(tài)的相位發(fā)生了突變，而x模態(tài)的相位是連續(xù)的.這些現(xiàn)象都與文獻(xiàn)［6］中的實(shí)驗(yàn)結(jié)果一致.

圖3 初始條件為x（0）=14.1 cm，z（0）=10.9 cm下的振動(dòng)曲線及相位情況

圖4 初始條件為x（0）=19.6 cm，z（0）=5.74 cm下的振動(dòng)曲線及相位情況

由于我們只關(guān)心各能量之間的相對(duì)關(guān)系，而振動(dòng)系統(tǒng)的能量與振幅的平方成正比，所以可以用振幅的平方代替能量.圖5和圖6畫(huà)出了不同初始條件下兩種模態(tài)能量隨時(shí)間的變化曲線，需要說(shuō)明的是，圖（5）和圖（6）中的縱坐標(biāo)實(shí)際上是振幅的平方.

圖5 初始條件為x（0）=14.1 cm，z（0）=10.9 cm下的相對(duì)能量

圖6 初始條件為x（0）=19.6 cm，z（0）=5.74 cm下的相對(duì)能量

由圖（5）和圖（6）可知，隨著時(shí)間的推移，能量在兩個(gè)模態(tài)之間不斷傳遞和轉(zhuǎn)換，因?yàn)橛凶枇Φ拇嬖?，能量并不守恒，總能量逐漸減小.在初始條件1的振動(dòng)中，兩個(gè)模態(tài)的能量變化幅度較大，而在初始條件2的振動(dòng)中，能量的變化幅度則相對(duì)較小，能量曲線相對(duì)平緩.圖中黑色曲線為理論模擬結(jié)果，點(diǎn)為文獻(xiàn)［6］的實(shí)驗(yàn)值，可以看出，兩種實(shí)驗(yàn)條件下，理論值和實(shí)驗(yàn)值都吻合得非常好.

4 結(jié)論

本文應(yīng)用一種彈簧擺系統(tǒng)模型，通過(guò)計(jì)算機(jī)數(shù)值模擬方法分析了彈簧擺的內(nèi)共振現(xiàn)象，通過(guò)描繪運(yùn)動(dòng)軌跡和振動(dòng)曲線，將彈簧擺的運(yùn)動(dòng)特征直觀地呈現(xiàn)出來(lái)，數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果表現(xiàn)出高度的一致性，再一次驗(yàn)證了應(yīng)用Mathematica對(duì)系統(tǒng)動(dòng)力學(xué)數(shù)值仿真的正確性與可靠性.