橢圓曲面擺參數振動的數值計算及仿真*

何烈云(浙江警察學院交通管理工程系 浙江 杭州 310053)李國軍(浙江警察學院公共基礎部 浙江 杭州 310053)

根據傅利葉變換所有周期振動均可分解成幾個(甚至無窮多個)簡諧振動，因此簡諧振動是一種簡單的機械振動形式，任何一本大學物理教材力學部分都將簡諧振動作為機械振動基礎內容介紹.典型的簡諧振動模型裝置有彈簧振子和單擺，除此之外，在習題中還會遇到球面擺問題，小球在一個半徑很大的光滑凹球面上做較小幅度擺動，討論小球擺動規律.教師在球面擺拓展教學中，往往從增大擺角或增加小球擺動的阻尼作用兩個方面增加問題的復雜程度.在一次講授球面擺習題時，有學生問及若將球面換成橢圓曲面，小球振動規律又會如何變化？小球在兩種不同曲面上擺動時最大區別在于“擺長”是否恒定，擺長指的是球心距離橢圓中心距離.在橢圓曲面上擺動時擺長符合一定規律變化，因此小球擺動屬于一種參數振動，不滿足簡諧振動的條件.

參數振動求解較為復雜，一般的力學教材沒有作深入介紹，但國內外許多學者們對不同類型的參數振動研究比較感興趣.文獻[1]將量子高級絕熱近似方法用來求解變擺長單擺的振動，研究表明，當擺長緩慢變化時，其振動形式近似為簡諧振動.文獻[2]研究表明，從能量轉換角度，任意振幅的單擺均做周期運動，其振動周期可采用橢圓積分形式表示.由于橢圓積分形式運用起來不夠方便，許多學者在此基礎上構造了不同的近似函數描述單擺在任意角度下的振動周期.文獻[3]構造的近似函數最為簡單，該函數只要擺幅控制在90°以內，單擺周期數值計算結果誤差不到1‰.文獻[4]將橢圓積分形式作泰勒展開，構造了新的單擺周期近似公式，采用該公式計算單擺周期時，只要振幅在114°的范圍內，單擺周期數值計算結果相對誤差小于 1‰.文獻[5]重點研究了阻尼及外界驅動力對單擺周期的影響，并分析了具體的運動形式.文獻[6]通過單擺運動與鉛直運動耦合機械模型，研究了變擺長單擺運動特性，通過高精度數值解法得到了變擺長單擺的近似解.文獻[7]利用龍格-庫達法，結合MATLAB軟件得到了非線性單擺振動的多種圖像.

上述研究文獻內容都屬于參數振動，且都是針對單擺從不同方面作了拓展.參數振動形式是多種多樣的，本文以橢圓曲面擺為研究對象，構建小球擺動的動力學方程，運用數值分析和計算機數值仿真的方法，討論小球在橢圓曲面擺動的運動規律，分析影響小球擺動周期的因素.

1 橢圓曲面擺動力學方程及求解

橢圓曲面擺裝置是將球面擺裝置中的球面改成了橢圓曲面，該模型裝置剖面圖如圖1所示.圖中O為平衡位置，O′為橢圓的幾何中心點，θ為小球擺角；d0為平衡位置到橢圓幾何中心的距離，即橢圓豎直方向半軸長度，d是小球球心與橢圓中心的距離(下稱擺長)；G為小球重量，Fn為小球受到曲面的作用力，f為小球擺動時受到的阻力，在速度不快時，可認為阻力與速度大小成正比且方向相反[8]，比例系數記為γ.

圖1 橢圓曲面擺裝置示意圖

1.1 橢圓曲面擺動力學方程

設小球質量為m，根據小球的受力特點，運用定軸轉動定律，擺動過程中的動力學微分方程滿足

(1)

由橢圓方程可知，擺長d滿足

(2)

式(2)中a為橢圓長半軸長度，b為短半軸長度.

令阻尼系數

聯立式(1)和式(2)，可得到小球擺動的動力學微分方程

(3)

1.2 橢圓曲面擺動力學方程數值計算

根據二階微分方程知識，式(3)屬于二階變系數非齊次微分方程，無法給出確切的解析解.若給定初始條件，式(3)可以采用龍格-庫達法給出方程的計算值，并利用計算機進行數值仿真，省去求解微分方程的復雜過程.該方法基本思想是在已知方程導數和初值初始條件，用區間上的若干點進行加權平均得到平均斜率，并經過多次迭代，獲得方程的數值解[9]，其中四階龍格-庫塔法最為經典，具有精度高、誤差小的優點[10].式(3)四階龍格-庫塔法計算基本步驟如下.首先將式(3)用角速度參數表示降為一階微分方程

(4)

此時θ=f(θ,t)，ω=g(θ,t)，設初始條件θ(t0)=θ0，ω(t0)=ω0，確定時間間隔h=δt，采用迭代法依次求出f(θ,t)和ω=g(θ,t)曲線斜率k.最后采用顯式龍格-庫達法式(5)和式(6)求得相應位置角位移和角速度數值

(5)

(6)

式中k1j和k2j分別代表f(θ,t)和ω=g(θ,t)曲線斜率.為了提高計算精度，可以縮短時間間隔h值.四階龍格-庫達法數值計算量大，在實踐中一般是采用計算機程序語言實現，并繪制f(θ,t) 和ω=g(θ,t)曲線圖像.

2 橢圓曲面擺基本振動規律

設橢圓水平方向長半軸a=2 m，豎直方向短半軸b=1 m，設初始條件ω0=0，θ0=arctan(0.2)，阻尼系數β=0.運用MATLAB繪制液體晃動角速度和角加速度圖像，得到如圖2所示的θ-t圖像、ω-t圖像和θ-ω圖像.

(a) θ-t圖像

(b) ω-t圖像

(c) θ-ω圖像

圖2(a)、(b)、(c)圖像可知，在不考慮阻尼作用時，橢圓曲面擺小球擺動類似簡諧振動，小球運動規律可以用正弦或余弦函數表示，圖像振幅和相位取決于初始條件值.因此無阻尼橢圓曲面擺的解形式可以寫成

(7)

式(7)中，φ為擺幅，T為小球擺動周期，由圖像可知T≈1.91 s.

3 橢圓曲面擺周期影響因素

接下來分別通過改變阻尼系數、橢圓曲面幾何尺寸及振幅大小，討論橢圓曲面擺周期特性.

3.1 阻尼系數對周期的影響

設橢圓長半軸a=2 m，短半軸b=1 m，初始條件ω0=0和θ0=arctan(0.2)，阻尼系數β大小分別為0，2，5.27，10，繪制圖3(a)、(b)橢圓曲面擺的θ-t圖像和θ-ω圖像.

(a) θ-t圖像

(b) θ-ω圖像

圖3表明：

(1)當β=2處于欠阻尼狀態，小球的擺幅隨著時間推移逐漸減小，直到最終停留在平衡位置處，而且振動周期大于同一橢圓曲面擺裝置無阻尼振動周期T；

(3)當β=10時處于過阻尼狀態，小球從開始擺動后慢慢回到平衡位置并停下來.

3.2 橢圓曲面幾何尺寸對周期的影響

討論橢圓曲面幾何尺寸對周期的影響時設β=0，初始條件ω0=0，θ0=arctan(0.2)，分兩種情況：

表1 短半軸固定，長半軸增大

表2 長半軸固定，短半軸增大

最終得到圖4(a)、(b)所示的θ-t圖像.

(a) 短半軸固定，長半軸變化

(b) 長半軸固定，短半軸變化

圖4表明，橢圓的長半軸和短半軸長度增大，小球擺動周期均增長.但短半軸的變化對小球擺動周期影響顯著，而長軸的變化對小球擺動周期影響極小.

3.3 振幅對小球擺動周期的影響

圖5 振幅變化的θ-t圖像

圖5表明，伴隨著擺動振幅增大，小球的擺動周期也變長.

4 結束語

(1)在橢圓曲面擺裝置中，小球擺動過程中擺長按一定規律變化，屬于一種參數振動.

(2)小球擺動的動力學方程為二階變系數齊次微分動力學方程，運用龍格-庫達法可以給出方程進行數值解，結合MATLAB軟件繪制角位移和角速度圖像.

(3)橢圓曲面擺振動與單擺振動類似，在無阻尼狀態時符合正弦變化規律，伴隨著阻尼系數、振幅、橢圓曲面長短軸增大，小球擺動周期也相應增大.