丁紅勝，王佳曦，張師平，董軍軍

(1. 北京科技大學(xué) 數(shù)理學(xué)院應(yīng)用物理系，北京 100083；2. 北京科技大學(xué) 自然科學(xué)基礎(chǔ)實(shí)驗(yàn)中心，北京 100083)

振動(dòng)是一種普遍存在于客觀世界的物質(zhì)運(yùn)動(dòng)形式，系統(tǒng)在周期性外力的持續(xù)作用下的振動(dòng)即為受迫振動(dòng)[1]．波耳共振儀是高校物理實(shí)驗(yàn)教學(xué)中研究受迫振動(dòng)運(yùn)動(dòng)規(guī)律較為普遍的儀器，也是學(xué)者探討較多的實(shí)驗(yàn)項(xiàng)目[2-4]．實(shí)驗(yàn)測(cè)量中，儀器的軸承摩擦和彈簧非線性效應(yīng)的影響是造成實(shí)驗(yàn)數(shù)據(jù)偏離預(yù)期的主要因素．如果在受迫振動(dòng)的分析中考慮這兩種因素的影響，則非常有助于學(xué)生理解實(shí)驗(yàn)原理和結(jié)果偏差．本文基于受迫振動(dòng)常規(guī)實(shí)驗(yàn)，延伸探討其非線性特性，通過(guò)Matlab數(shù)值分析的方法來(lái)研究杜芬方程的混沌特性，以此提高學(xué)生對(duì)受迫振動(dòng)的認(rèn)知，并加深對(duì)其非線性特性的理解．

1 受迫振動(dòng)

1.1 受迫振動(dòng)方程

波耳共振儀[5](如圖1所示)的擺輪可在彈性力矩作用下自由擺動(dòng)，若同時(shí)加上阻尼力矩和驅(qū)動(dòng)力矩，擺輪可做受迫振動(dòng)．

1. 光電門(mén)A；2. 長(zhǎng)凹槽；3. 短凹槽；4. 銅制擺輪；5. 搖桿； 6. 蝸卷彈簧；7. 機(jī)架；8. 阻尼線圈；9. 連桿；10. 搖桿調(diào)節(jié)螺絲； 11. 光電門(mén)B；12. 角度盤(pán)；13. 有機(jī)玻璃轉(zhuǎn)盤(pán)；14. 底座； 15. 彈簧夾持螺絲；16閃光燈圖1 波耳共振儀的結(jié)構(gòu)示意圖[5]

(1)

式中，J為擺輪的轉(zhuǎn)動(dòng)慣量，-kθ為彈性力矩，k為彈簧的勁度系數(shù)，M0為驅(qū)動(dòng)力矩的幅值，ω為驅(qū)動(dòng)力的角頻率．

(2)

當(dāng)mcosωt=0時(shí)，即在無(wú)周期性外力矩作用時(shí)，式(2)為阻尼振動(dòng)方程；當(dāng)阻尼系數(shù)為零，即無(wú)阻尼時(shí)，式(2)為簡(jiǎn)諧振動(dòng)方程，ω0即為振動(dòng)系統(tǒng)的固有頻率．

方程(2)的通解為

θ=θ1e-βtcos(ω1t+α)+θ2cos(ωt+φ)

(3)

通解含有兩個(gè)部分：θ1e-βtcos(ω1t+α)表示阻尼振動(dòng)，隨著時(shí)間的演化，振動(dòng)會(huì)逐漸衰減，直至忽略不計(jì)；θ2cos(ωt+φ)表示簡(jiǎn)諧振動(dòng)，驅(qū)動(dòng)力矩持續(xù)對(duì)擺輪做功，振動(dòng)系統(tǒng)接收驅(qū)動(dòng)力所傳送的能量，使振動(dòng)系統(tǒng)最終達(dá)到一個(gè)穩(wěn)定的狀態(tài)．

振幅為定值，

(4)

系統(tǒng)相位差為

(5)

由式(4)和式(5)可看出，穩(wěn)定振動(dòng)時(shí)，驅(qū)動(dòng)力的幅值M0、阻尼系數(shù)β、驅(qū)動(dòng)力的頻率ω(或驅(qū)動(dòng)力周期T)和系統(tǒng)固有頻率ω0(或固有周期T0)4個(gè)參量共同決定系統(tǒng)振幅和相位差的大小，與振動(dòng)的起始狀態(tài)無(wú)關(guān)．

(6)

可以看出，阻尼系數(shù)越小，驅(qū)動(dòng)力的角頻率越接近系統(tǒng)的固有角頻率，振動(dòng)幅值也越大．

1.2 阻尼系數(shù)、固有周期、固有頻率的測(cè)量

實(shí)驗(yàn)中通常選擇阻尼開(kāi)關(guān)位置為2，將波耳共振儀有機(jī)玻璃盤(pán)上零度標(biāo)志線對(duì)準(zhǔn)0°，用手將擺輪轉(zhuǎn)動(dòng)至140°～150°之間，放手后控制箱會(huì)自動(dòng)連續(xù)記錄擺輪完成阻尼振動(dòng)10次的振幅θ0～θ9．?dāng)?shù)據(jù)如表1所示．

表1 阻尼系數(shù)測(cè)量數(shù)據(jù)

選擇“自由振蕩”進(jìn)行測(cè)量，將有機(jī)玻璃盤(pán)上零度標(biāo)志線保持在0°，用手將擺輪轉(zhuǎn)動(dòng)至140°～150°之間，放手后記錄每次振動(dòng)振幅值及其相應(yīng)周期，并計(jì)算相應(yīng)頻率值．?dāng)?shù)據(jù)如表2所示．

表2 振幅、周期、頻率關(guān)系數(shù)據(jù)

2 受迫振動(dòng)方程非線性修正

2.1 彈簧非線性效應(yīng)

一般情況下，可以認(rèn)為彈簧的勁度系數(shù)k為常數(shù)，但是在上述彈簧自由振動(dòng)的實(shí)驗(yàn)中，通常隨著振幅的減小，振動(dòng)系統(tǒng)固有周期會(huì)逐漸增大，這意味著振動(dòng)系統(tǒng)的彈性回復(fù)力隨著振幅的減小而減弱[2]．因此不妨假設(shè)彈簧的勁度系數(shù)為振幅的函數(shù)k=k(θ)，則回復(fù)力矩為

Fk=k(θ)θ

(7)

將Fk=k(θ)θ進(jìn)行泰勒展開(kāi)至一階項(xiàng)，

(8)

(9)

2.2 軸承摩擦效應(yīng)

軸承摩擦效應(yīng)普遍存在于定軸轉(zhuǎn)動(dòng)的運(yùn)動(dòng)中，當(dāng)軸承轉(zhuǎn)動(dòng)速度較小時(shí)軸承靜摩擦力較大，而當(dāng)轉(zhuǎn)速達(dá)到一定臨界值后，軸承摩擦力將保持穩(wěn)定[6]．

實(shí)驗(yàn)發(fā)現(xiàn)，軸承摩擦力隨軸承角速度的變化率較緩，因此可采用二次函數(shù)的形式來(lái)描述軸承摩擦力矩隨角速度的關(guān)系

(10)

其中，固定系數(shù)B>0，C>0，在此情況下，即表明摩擦力矩隨角速度的增加而緩慢降低．

2.3 阻尼振動(dòng)方程非線性修正

綜合彈簧的非線性效應(yīng)與軸承摩擦效應(yīng)對(duì)振動(dòng)方程的修正，可以得到修正后的阻尼振動(dòng)方程為

(11)

(12)

式(12)即為非線性項(xiàng)修正后的阻尼振動(dòng)方程，是一個(gè)待定系數(shù)的二階微分方程．

在此情況下，基于Matlab軟件，用數(shù)值方法計(jì)算其數(shù)值解，通過(guò)多次給定a、b、c的數(shù)值進(jìn)行模擬計(jì)算，擬合阻尼振動(dòng)實(shí)驗(yàn)所測(cè)得的數(shù)據(jù)．

(13)

利用Matlab軟件內(nèi)置的ode45函數(shù)，根據(jù)四階龍格-庫(kù)塔法，多次變步長(zhǎng)數(shù)值積分，可以得到，當(dāng)待定系數(shù)滿足a=0.014,b=0.002,c=0.005時(shí)，數(shù)值方法解得的阻尼振動(dòng)方程與實(shí)驗(yàn)數(shù)據(jù)擬合的最為接近，特別是在振幅為90°～140°時(shí)，彈簧形變較大，彈簧非線性效應(yīng)較為明顯的區(qū)域．

圖2中指數(shù)曲線為阻尼振動(dòng)的振幅-時(shí)間實(shí)測(cè)數(shù)據(jù)的擬合曲線，周期性曲線為非線性微分方程的數(shù)值解．

圖2 阻尼振動(dòng)數(shù)值積分與擬合曲線

2.4 受迫振動(dòng)方程非線性修正

在修正后的阻尼振動(dòng)方程(12)的右邊加上一個(gè)周期性力矩，即為受迫振動(dòng)方程，

(14)

該方程中m為未知參量，需要通過(guò)實(shí)驗(yàn)儀器進(jìn)行數(shù)據(jù)測(cè)量加以確定．

選取“阻尼2”檔位，打開(kāi)電機(jī)，設(shè)定在某一轉(zhuǎn)速，當(dāng)觀察到擺輪的受迫振動(dòng)周期與電機(jī)的周期趨于一致時(shí)，即受迫振動(dòng)達(dá)到穩(wěn)定時(shí)，開(kāi)始測(cè)量．讀出每次擺輪受迫振動(dòng)的幅值θ和驅(qū)動(dòng)力矩每10次振動(dòng)周期10T，利用頻閃法測(cè)量擺輪受迫振動(dòng)的位移與驅(qū)動(dòng)力之間的相位差φ，并記錄此時(shí)電動(dòng)機(jī)的轉(zhuǎn)速值，結(jié)果如表3所示.

表3 幅頻特性與相頻特性測(cè)量數(shù)據(jù)

至此獲得了在一定情況下，受迫振動(dòng)動(dòng)力學(xué)方程的全部參量：ω0=3.775 rad·s-1，β=0.081，m=1.02，a=0.014，b=-0.002，c=0.005．

(15)

把上述一階微分方程組輸入程序，利用數(shù)值積分求解其動(dòng)力學(xué)特性，結(jié)果如圖3所示．

圖3 系統(tǒng)時(shí)間響應(yīng)與相平面曲線

圖3中左側(cè)為受迫振動(dòng)系統(tǒng)的時(shí)間響應(yīng)曲線，右側(cè)為系統(tǒng)的相平面曲線．受迫振動(dòng)系統(tǒng)的時(shí)間響應(yīng)曲線表示了振動(dòng)振幅隨時(shí)間變化的特性，可以看出，隨著振動(dòng)時(shí)間的增加，振動(dòng)頻率逐漸穩(wěn)定，振幅逐漸增加直到最大值并且保持穩(wěn)定．

受迫振動(dòng)系統(tǒng)的相平面曲線圖表示了速度隨振幅變化的關(guān)系．可以看出隨著時(shí)間的增加，速度與振幅也會(huì)趨近一個(gè)穩(wěn)定狀態(tài)．

如果改變其中的一個(gè)參量，例如將與彈簧非線性效應(yīng)相關(guān)的參數(shù)a從0.014調(diào)整為10，采用同樣的方法對(duì)系統(tǒng)進(jìn)行數(shù)值分析，得到的結(jié)果如圖4所示．原本穩(wěn)定的受迫振動(dòng)系統(tǒng)在振動(dòng)后期開(kāi)始出現(xiàn)一定的隨機(jī)運(yùn)動(dòng)，即混沌現(xiàn)象．

圖4 系統(tǒng)時(shí)間響應(yīng)與相平面曲線

基于波耳共振儀的受迫振動(dòng)，由于其振動(dòng)方程中涉及的非線性項(xiàng)偏多，不利于定性和定量分析受迫振動(dòng)系統(tǒng)中的混沌現(xiàn)象，因此可以選擇典型的非線性受迫振動(dòng)方程——杜芬方程進(jìn)行探討．

3 受迫振動(dòng)方程的混沌

3.1 杜芬方程

振動(dòng)問(wèn)題中的很多數(shù)學(xué)模型都可以通過(guò)轉(zhuǎn)化為杜芬方程來(lái)進(jìn)行研究和分析，杜芬方程是研究和分析某些非線性動(dòng)力學(xué)系統(tǒng)的基礎(chǔ)．

杜芬方程[7]的標(biāo)準(zhǔn)形式為

(16)

其中，δ>0，g(x)是含x3項(xiàng)的函數(shù)，f(x,t)是周期性函數(shù)．

標(biāo)準(zhǔn)化后的杜芬方程稱之為杜芬系統(tǒng)．杜芬系統(tǒng)是一類從簡(jiǎn)單受迫振動(dòng)的物理模型中歸納出的非線性振動(dòng)系統(tǒng)模型，形式簡(jiǎn)潔，具有代表性，其應(yīng)用范圍十分廣泛．許多動(dòng)力學(xué)過(guò)程都與杜芬系統(tǒng)的模型極其相似，比如化學(xué)鍵的斷裂、建筑結(jié)構(gòu)的震顫和輪船的橫搖運(yùn)動(dòng)、電路周期振蕩、電路信號(hào)檢測(cè)、電路系統(tǒng)模擬分析和電路的故障檢測(cè)等等．

3.2 硬彈簧型杜芬方程

杜芬方程一般可以分成四種基本類型，可以選取其中的硬彈簧型杜芬方程[8]進(jìn)行探討，

(17)

將式(17)轉(zhuǎn)化成一階微分方程組：

(18)

其中ω、f、a、β、c為可變參量．

給定不同的初始條件，利用數(shù)值計(jì)算的方法來(lái)分析硬彈簧型杜芬系統(tǒng)的受迫振動(dòng)．

當(dāng)振動(dòng)時(shí)間為300 s時(shí)，結(jié)果如圖5所示，從上至下依次是不同f取值下的系統(tǒng)時(shí)間響應(yīng)圖、系統(tǒng)相平面軌跡圖、系統(tǒng)功率密度譜．

圖5 系統(tǒng)時(shí)間響應(yīng)、相平面軌跡、功率密度譜隨f值的演化

由圖5可以看出，隨著系統(tǒng)參量f的變化，系統(tǒng)逐漸由周期性運(yùn)動(dòng)進(jìn)入混沌狀態(tài)．

3.3 杜芬方程的混沌

從圖5的相平面軌跡圖可以看出，系統(tǒng)從只有一個(gè)相對(duì)穩(wěn)定的平衡位置(原點(diǎn))，逐漸變?yōu)榫哂袃蓚€(gè)相對(duì)平衡的位置(焦點(diǎn))，這體現(xiàn)了無(wú)序混沌狀態(tài)中存在一定程度上的有序特性．相平面軌跡圖中相對(duì)平衡位置發(fā)生改變的同時(shí)，相軌跡也從圓形逐漸轉(zhuǎn)化為類似于莫比烏斯環(huán)的橢圓形狀．進(jìn)一步推測(cè)，系統(tǒng)會(huì)擁有兩個(gè)相對(duì)平衡的焦點(diǎn)，相軌跡圍繞這兩個(gè)焦點(diǎn)運(yùn)動(dòng)，而且可將這兩個(gè)焦點(diǎn)類比為三維非線性系統(tǒng)中奇異吸引子[9]，如圖6所示．

圖6 奇異吸引子

奇異吸引子是混沌系統(tǒng)中獨(dú)有的一種吸引子，也被稱之為混沌吸引子，它代表著混沌狀態(tài)中的一種定態(tài)．具體表現(xiàn)在系統(tǒng)從任一初始狀態(tài)開(kāi)始，最終都會(huì)被奇異吸引子吸引到相空間中的某一特定區(qū)域．奇異吸引子與一般吸引子有所區(qū)別，當(dāng)混沌狀態(tài)的軌跡線進(jìn)入奇異吸引子之后，軌跡線并不會(huì)像一般吸引子一樣，吸引軌跡線繼續(xù)圍繞其運(yùn)動(dòng)，最終形成一個(gè)相對(duì)封閉的圖形，而是會(huì)讓兩條接近的軌跡線發(fā)生指數(shù)分離．從系統(tǒng)外部看，奇異吸引子是在聚集軌跡線；從系統(tǒng)內(nèi)部看，是軌跡線發(fā)散的過(guò)程．奇異吸引子常常成對(duì)出現(xiàn)，兩個(gè)奇異吸引子同時(shí)對(duì)軌跡線的作用，共同構(gòu)筑了混沌系統(tǒng)一定程度上的有序性．

圖7 從左至右分別為f=1, f=3, f=5時(shí)的相空間軌跡線

從圖7可以看出，杜芬方程中的兩個(gè)焦點(diǎn)在含有時(shí)間項(xiàng)的坐標(biāo)下，展開(kāi)成了兩條軸線，體現(xiàn)了其時(shí)間平移不變性，參量f的變化引起相空間軌跡明顯改變，表明混沌系統(tǒng)對(duì)參量具有很強(qiáng)的敏感性，但在速度-振幅的相平面軌跡圖中，兩個(gè)焦點(diǎn)是穩(wěn)定的，不隨時(shí)間發(fā)生變化，這與奇異吸引子的性質(zhì)非常類似．

4 結(jié)論

基于波耳共振儀受迫振動(dòng)實(shí)驗(yàn)，探討其中的非線性影響因素，對(duì)阻尼振動(dòng)方程進(jìn)行修正，考慮軸承摩擦和彈簧非線性效應(yīng)，結(jié)合Matlab數(shù)值計(jì)算結(jié)果與波耳共振儀實(shí)測(cè)數(shù)據(jù)，確定波耳共振儀系統(tǒng)的非線性項(xiàng)參數(shù)，分析其非線性效應(yīng)．通過(guò)引入非線性理論的思想和數(shù)值分析的計(jì)算方法，有助于學(xué)生準(zhǔn)確認(rèn)知實(shí)驗(yàn)現(xiàn)象，加深對(duì)實(shí)驗(yàn)原理的理解．

為了豐富實(shí)驗(yàn)教學(xué)內(nèi)容，引入彈簧型杜芬方程，分析其非線性特性，利用數(shù)值計(jì)算方法探討受迫振動(dòng)模型進(jìn)入混沌狀態(tài)的條件及其運(yùn)動(dòng)特性．受迫振動(dòng)是否會(huì)進(jìn)入混沌狀態(tài)，可以根據(jù)系統(tǒng)響應(yīng)時(shí)間曲線(振幅-時(shí)間)、系統(tǒng)相平面曲線(速度-振幅)和功率密度譜曲線(功率-頻率)綜合判斷．如果受迫振動(dòng)做周期運(yùn)動(dòng)，其響應(yīng)時(shí)間曲線會(huì)呈現(xiàn)規(guī)則的周期性，其相軌跡圖具有極限環(huán)的特性，其功率圖譜峰值呈現(xiàn)離散的狀態(tài)；如果系統(tǒng)進(jìn)入混沌狀態(tài)，其時(shí)間曲線會(huì)具有無(wú)規(guī)則的、類似于噪聲的特性，其相平面曲線在有限區(qū)域內(nèi)是無(wú)周期性的，隨著振動(dòng)時(shí)間的延長(zhǎng)，其軌跡線會(huì)逐漸鋪滿整個(gè)區(qū)域，功率圖譜含有多個(gè)峰值；如果系統(tǒng)的響應(yīng)時(shí)間曲線和相平面曲線在上述兩種情況之間，則功率圖譜呈現(xiàn)離散狀態(tài)，系統(tǒng)做的是擬周期運(yùn)動(dòng)，如果改變它的系統(tǒng)參量或者初始條件，則有可能轉(zhuǎn)變?yōu)橹芷谶\(yùn)動(dòng)或混沌運(yùn)動(dòng)．一般情況下，受迫振動(dòng)受非線性效應(yīng)影響越大，最終進(jìn)入混沌狀態(tài)的可能性越高．