林炯輝，劉玉穎，宋 敏

(1.中國農(nóng)業(yè)大學(xué) 工學(xué)院，北京 100083；2.中國農(nóng)業(yè)大學(xué) 理學(xué)院，北京 100083)

自17世紀(jì)牛頓經(jīng)典力學(xué)建立以來，就一直得到普遍的應(yīng)用，它打開了近代自然科學(xué)的大門，在其基礎(chǔ)上發(fā)展出了彈性力學(xué)、流體力學(xué)、結(jié)構(gòu)力學(xué)等多門學(xué)科. 在牛頓力學(xué)中，如果已知物體的受力和初始狀態(tài)，對于有序系統(tǒng)其未來的運(yùn)動(dòng)狀態(tài)是確定的. 學(xué)者們根據(jù)經(jīng)典力學(xué)的定律和萬有引力定律曾經(jīng)精確地預(yù)言彗星和小行星等的運(yùn)動(dòng)，并且得到了驗(yàn)證.

近年來研究者們在大學(xué)物理教學(xué)中引入Matlab計(jì)算機(jī)數(shù)值模擬，培養(yǎng)學(xué)生用計(jì)算機(jī)進(jìn)行學(xué)習(xí)和研究物理問題的能力[1-4]，將VPython引進(jìn)大學(xué)物理課堂的方式也逐漸興起[5-7]. VPython可以利用簡單的代碼模擬物體在三維空間中的運(yùn)動(dòng)，實(shí)時(shí)觀測物體的速度變化及其運(yùn)動(dòng)軌跡. 質(zhì)點(diǎn)平動(dòng)、剛體定軸轉(zhuǎn)動(dòng)和復(fù)擺振動(dòng)是基本的運(yùn)動(dòng)形式，本文選用VPython軟件對剛體的平動(dòng)、轉(zhuǎn)動(dòng)、振動(dòng)進(jìn)行模擬與探究，利用轉(zhuǎn)動(dòng)定律、角動(dòng)量定理、動(dòng)量定理等基本規(guī)律，從剛體的定軸轉(zhuǎn)動(dòng)出發(fā)，研究其在恒定力矩和變化力矩作用下的運(yùn)動(dòng)狀態(tài)以及其能量特征，進(jìn)而研究剛體的平面平行運(yùn)動(dòng)；通過模擬不同角度復(fù)擺的運(yùn)動(dòng)，觀察復(fù)擺受迫振動(dòng)時(shí)的運(yùn)動(dòng)特征，檢驗(yàn)相關(guān)結(jié)論.

1 Vpython概況

VPython是Python默認(rèn)的3D模塊，有著和Python一樣風(fēng)格的編程方式，可以快速創(chuàng)建三維場景和動(dòng)畫. Vpython可創(chuàng)建多種實(shí)體類型，如球體、圓柱體、長方體、錐體等，通過代碼改變其在空間中的位置即可實(shí)現(xiàn)運(yùn)動(dòng)模擬. 利用VPython模擬剛體轉(zhuǎn)動(dòng)的代碼框架為：① 定義常數(shù)(m，g，k等)和運(yùn)動(dòng)初始值(初始角位置，初始角動(dòng)量等)；② 設(shè)置步長Δt，Δt越小模擬越精確；③ 進(jìn)入while循環(huán)，在循環(huán)中根據(jù)當(dāng)前狀態(tài)計(jì)算合外力矩；④ 用該力矩更新系統(tǒng)角動(dòng)量；⑤ 利用角動(dòng)量更新系統(tǒng)角位置. 程序中以差分代替微分，得到角速度、角位置對時(shí)間的迭代式(Euler-Cromer方法)：

ωi+1=ωi+f(θi,ωi,ti)Δt，

θi+1=θi+ωi+1Δt，

ti+1=ti+Δt

Δt越小，模擬的精度越高.

2 剛體轉(zhuǎn)動(dòng)的VPython模擬

2.1 轉(zhuǎn)動(dòng)定律簡介

在大學(xué)物理的學(xué)習(xí)中，我們依據(jù)牛頓第二定律導(dǎo)出了剛體定軸轉(zhuǎn)動(dòng)定律M=Jα，即剛體所受的對于某定軸的合外力矩等于剛體對此定軸的轉(zhuǎn)動(dòng)慣量與剛體在此合外力矩作用下所獲得的角加速度的乘積；力矩的時(shí)間累積效應(yīng)——角動(dòng)量定理dL=Mdt，即剛體角動(dòng)量的增量等于所受合外力矩的角沖量.

2.2 恒定力矩下勻質(zhì)圓柱體的定軸轉(zhuǎn)動(dòng)

有一個(gè)質(zhì)量為2 kg，長度為1 m，半徑為0.1 m的勻質(zhì)圓柱體，轉(zhuǎn)軸為過其質(zhì)心且與圓柱體中心軸線垂直. 勻質(zhì)圓柱體最初處于靜止?fàn)顟B(tài)，對其施加一恒定力矩，力矩方向沿轉(zhuǎn)軸正方向，大小為2 Nm,對其運(yùn)動(dòng)進(jìn)行分析.

對于勻質(zhì)圓柱體，轉(zhuǎn)軸(圖中所示z軸)位于圓柱體中心且垂直于中心軸線，建立如圖1所示的坐標(biāo)系.

圖1 勻質(zhì)圓柱體，其轉(zhuǎn)軸為z軸

(1)

我們在VPython中對該勻質(zhì)圓柱體的轉(zhuǎn)動(dòng)進(jìn)行模擬，利用VPython顯示勻質(zhì)圓柱體尺寸、顏色以及轉(zhuǎn)軸位置(如圖2所示)，將棒的運(yùn)動(dòng)可視化.

圖2 勻質(zhì)圓柱體繞定軸轉(zhuǎn)動(dòng)及對應(yīng)坐標(biāo)系，轉(zhuǎn)軸沿z軸正向

在VPython模擬中，可以直觀地觀察它的運(yùn)動(dòng)過程，程序中使用的物理公式及對應(yīng)的程序循環(huán)語句如下.

whilet<2:

rate(100)

torque=vector(0,0,2) #M=(0,0,2)

L=L+torque*dt#dL=M×dt

dtheta=omega*dt#dθ=|ω|×dt

theta=theta+dtheta

rod.rotate(angle=dtheta,axis=(0,0,1),origin= axle.pos)

t=t+dt

VPython模擬結(jié)果表明：t=2 s時(shí)，角位置θ=23.3044661359 rad，角速度ω=23.3021359223 rad/s，ω/θ=0.999900，與理論計(jì)算結(jié)果符合得很好.

本文還探究了勻質(zhì)圓柱體轉(zhuǎn)動(dòng)過程中角速度、角位置隨時(shí)間的變化關(guān)系：由于恒定力矩作用，角加速度為一恒量，角速度隨時(shí)間線性增大，角位置為時(shí)間的二次函數(shù)，與時(shí)間t的平方成正比(如圖3所示).

圖3 角速度、角位置與時(shí)間關(guān)系圖像

2.3 正余弦變化的力矩下勻質(zhì)圓柱體的定軸轉(zhuǎn)動(dòng)

(t=0,ω=0)

(2)

-0.699cos(5t)+0.699,(t=0,θ=0)

(3)

進(jìn)行VPython模擬，將力矩定義式更改為torque=3cos(5t).模擬結(jié)果表明：角速度ω、角位置θ隨時(shí)間t均做周期性變化，角速度ω、角位置θ隨時(shí)間的變化曲線如圖4所示. 理論與模擬結(jié)果均表明角速度、角位置的變化周期均等于力矩變化周期；且在運(yùn)動(dòng)的一個(gè)周期內(nèi)，圓柱體先加速后減速(圖5a)，當(dāng)角速度減為0時(shí)，轉(zhuǎn)角達(dá)到最大(圖5b).隨后勻質(zhì)圓柱體反向加速(圖5c)到最大值再減速，減速至角速度為零時(shí)，此時(shí)θ=0，棒回到初始位置(圖5d).

圖4 角速度、角位置隨時(shí)間變化關(guān)系曲線

圖5 勻質(zhì)圓柱體轉(zhuǎn)動(dòng)的3D場景

力矩M和轉(zhuǎn)動(dòng)動(dòng)能E隨時(shí)間變化曲線如圖6所示，勻質(zhì)圓柱體動(dòng)能變化周期等于力矩周期的一半.動(dòng)能為

力矩大時(shí)，轉(zhuǎn)動(dòng)動(dòng)能不一定大.通過動(dòng)能與力矩的函數(shù)關(guān)系圖像(圖7)可知：當(dāng)力矩為極值±3Nm時(shí)，動(dòng)能最小為0，當(dāng)力矩為0時(shí)，動(dòng)能為最大值1.048 J.

圖6 動(dòng)能、力矩隨時(shí)間變化關(guān)系曲線

圖7 動(dòng)能與力矩關(guān)系圖像

2.4 平動(dòng)加轉(zhuǎn)動(dòng)

現(xiàn)在在轉(zhuǎn)動(dòng)軸上施加一恒力，勻質(zhì)圓柱體的運(yùn)動(dòng)變?yōu)橘|(zhì)心的平動(dòng)和繞過質(zhì)心軸的轉(zhuǎn)動(dòng).由動(dòng)量定理P2-P1=FΔt，得到動(dòng)量迭代式P=P+Fdt，由P=mv，求得任意時(shí)刻的瞬時(shí)平動(dòng)速度v和dt內(nèi)的位移ds=vdt.在while循環(huán)中更新轉(zhuǎn)軸的位置(設(shè)力沿x軸正方向，故只需改變勻質(zhì)圓柱體質(zhì)心x坐標(biāo)).

whilet<2:

P=P+F*dt

#dPtranslate=F*dt

v=P/M

rod.pos=rod.pos+v*dt

#dx=vdt

axle.pos.x=axle.pos.x+mag(v*dt)

在相同的力矩和大小不同的外力作用下，系統(tǒng)的運(yùn)動(dòng)軌跡不同，如圖8所示.力矩均為M=3cos5t，外力大小分別為0.1 N、0.5 N、1 N沿x軸正向(方向向右)，模擬結(jié)果表明：勻質(zhì)圓柱體端點(diǎn)的軌跡呈鋸齒狀，且隨著外力的增大，“鋸齒”越來越疏松(如圖8所示).

圖8 相同的力矩和不同外力作用下勻質(zhì)圓柱體的運(yùn)動(dòng)軌跡，力矩均為M=3cos5t

3 勻質(zhì)棒的擺動(dòng)(大角度、小角度)的VPython模擬

3.1 重力矩下的棒的擺動(dòng)(大角度、小角度)

在大學(xué)物理振動(dòng)學(xué)章節(jié)中，我們了解了簡諧振動(dòng)是最簡單、最基本的振動(dòng)，小角度復(fù)擺即是簡諧運(yùn)動(dòng)的一個(gè)特例. 在研究復(fù)擺的運(yùn)動(dòng)周期時(shí)，令物體作小振幅擺動(dòng)(擺角<5°). 根據(jù)轉(zhuǎn)動(dòng)定律得

(3)

式中l(wèi)為復(fù)擺質(zhì)心到轉(zhuǎn)軸的距離，當(dāng)θ較小時(shí)，可近似看作sinθ≈θ，式(3)變?yōu)?/p>

解微分方程得

θ=Acos(ωt+ψ)

振動(dòng)的圓頻率和周期分別為[8]

(4)

(5)

通過積分理論計(jì)算，得大擺角運(yùn)動(dòng)時(shí)復(fù)擺的周期為[10]

(6)

可以看出任意角度下擺的運(yùn)動(dòng)方程與運(yùn)動(dòng)周期的理論推導(dǎo)式都與小角度近似(sinθ≈θ)存在差異，θ0為大角度時(shí)，擺的運(yùn)動(dòng)已非簡諧運(yùn)動(dòng).在VPython中，對一根懸掛在低摩擦軸上的桿的運(yùn)動(dòng)(忽略阻力)進(jìn)行建模(圖9)，設(shè)桿的質(zhì)量為m，桿長為L，并設(shè)立一組遞增的初始偏置角度.模擬結(jié)果如圖10所示，虛線為小角度近似得到的余弦曲線:

θ=θ0cos(ω0t),

實(shí)線為擺桿實(shí)際角位移曲線.當(dāng)θ0為小角度時(shí)(<10°)，可見擺桿實(shí)際角位移圖線與簡諧運(yùn)動(dòng)角位移理論計(jì)算結(jié)果曲線相重合，擺桿的運(yùn)動(dòng)為簡諧運(yùn)動(dòng)；當(dāng)θ0>10°時(shí)，隨著初始角的增大，實(shí)際角位移曲線與簡諧運(yùn)動(dòng)角位移理論計(jì)算結(jié)果曲線偏差越大，擺桿不再作簡諧運(yùn)動(dòng)(如圖10(c)-(f)所示)；而當(dāng)θ0=180°時(shí)，擺桿達(dá)到往復(fù)擺動(dòng)與圓周運(yùn)動(dòng)臨界點(diǎn).

圖9 擺桿Vpython建模實(shí)體

圖10 不同初始偏轉(zhuǎn)角度下擺桿實(shí)際角位移曲線與簡諧運(yùn)動(dòng)角位移理論曲線對比圖

令縱坐標(biāo)表示角速度，橫坐標(biāo)表示角位置，圓棒的運(yùn)動(dòng)狀態(tài)對應(yīng)坐標(biāo)軸上的點(diǎn)，這樣的圖形叫做相圖.相圖沒有直觀地顯示角速度、角位置與時(shí)間的關(guān)系，卻能反映出運(yùn)動(dòng)的全局性.相圖法是非線性力學(xué)中最基本的研究方法[11].將不同初始條件下擺桿的運(yùn)動(dòng)狀態(tài)用ω-θ相圖表示如圖11所示.

圖11 不同初始條件下擺桿受重力矩的運(yùn)動(dòng)狀態(tài)ω-θ相圖

由圖11可知，當(dāng)θ0較小時(shí)，相圖軌線為圓形閉合曲線，表明運(yùn)動(dòng)為簡諧振動(dòng)；θ0越大，軌線越“方”.

3.2 考慮阻尼和周期性外力矩時(shí)擺運(yùn)動(dòng)特性的研究

在經(jīng)典力學(xué)中，已知物體初始狀態(tài)和物體運(yùn)動(dòng)過程中受力情況，可以通過牛頓運(yùn)動(dòng)定律和動(dòng)量定理預(yù)測物體的運(yùn)動(dòng)；但自然界中并非所有物體運(yùn)動(dòng)都是確定且唯一的，即使給定了初始條件，物體依然可能有多條運(yùn)動(dòng)軌跡，非線性系統(tǒng)就是典型例子.當(dāng)給任意角度復(fù)擺一個(gè)周期性外力矩并考慮阻力，使其受迫振動(dòng)，當(dāng)驅(qū)動(dòng)力矩達(dá)到一定值時(shí)擺的運(yùn)動(dòng)狀態(tài)無法預(yù)料，且對初始條件具有敏感性，這在物理上被稱為混沌.混沌的一種著名表述是“蝴蝶效應(yīng)”，初始條件的極小偏差，都將可能會(huì)引起結(jié)果的極大差異.

受迫振動(dòng)微分方程為

(7)

C為阻力系數(shù)，M為力矩最大值，ωp為力矩周期變化的角頻率，上式化簡為

(8)

圖12 截取不同初始條件下擺桿受迫阻尼振動(dòng)的ω -θ相圖

由于疊加原理、傅立葉變換等工具不適于求解非線性方程，而計(jì)算機(jī)性能的提高使得研究者一般采用計(jì)算機(jī)進(jìn)行迭代運(yùn)算對其進(jìn)行數(shù)值求解.通過有阻尼受迫振動(dòng)的模擬，發(fā)現(xiàn)在一定范圍內(nèi)即使初始條件f變化很小，卻會(huì)引起復(fù)擺的運(yùn)動(dòng)形式有很大的不同.例如當(dāng)f=1.505時(shí)，相軌跡為無規(guī)則的發(fā)散曲線；f=1.508時(shí)，相軌跡為收斂的三個(gè)極限環(huán).與當(dāng)年洛倫茲在利用計(jì)算機(jī)研究地球大氣時(shí)提出的“蝴蝶效應(yīng)”相類似，對于十分接近的兩個(gè)f的值，初始時(shí)相軌跡重疊，當(dāng)時(shí)間積累到一定程度后，相軌跡會(huì)出現(xiàn)較大偏差：設(shè)定模擬精度dt=0.001 s，f初始值相差0.00001時(shí)，程序運(yùn)行200 s后，兩個(gè)復(fù)擺的運(yùn)動(dòng)將截然不同.隨著f值的變化，相軌跡運(yùn)動(dòng)呈現(xiàn)出由收斂、有序走向混沌的狀態(tài).混沌狀態(tài)下相點(diǎn)軌跡時(shí)而做周期運(yùn)動(dòng)，時(shí)而做非周期運(yùn)動(dòng)；時(shí)而向外延伸，時(shí)而在某一位置纏繞.對于某一時(shí)刻而言復(fù)擺運(yùn)動(dòng)具有不確定性(如圖11(e)—(g)).如果不斷增大f值，研究f>1.7以后的狀態(tài)(不再贅述)，收斂的極限環(huán)與混沌相圖將繼續(xù)交替出現(xiàn)，說明復(fù)擺的周期性運(yùn)動(dòng)與混沌狀態(tài)交替出現(xiàn).

簡單的復(fù)擺系統(tǒng)實(shí)際上的運(yùn)動(dòng)卻“不簡單”.混沌的發(fā)現(xiàn)對經(jīng)典力學(xué)造成了沖擊，同時(shí)也令其不斷完善和發(fā)展. 隨著人們對混沌理論的深入了解，將會(huì)給我們生活帶來有利的應(yīng)用. Vpython模擬可以檢驗(yàn)非線性系統(tǒng)相關(guān)結(jié)論：1)疊加原理不再成立；2)初始條件的不同，會(huì)導(dǎo)致很大差異的運(yùn)動(dòng)形式；3)可能出現(xiàn)完全隨機(jī)的混沌狀態(tài).

4 結(jié)論

本文利用VPython對剛體運(yùn)動(dòng)和剛體擺動(dòng)進(jìn)行運(yùn)動(dòng)模擬和數(shù)值模擬. 利用轉(zhuǎn)動(dòng)定律和角動(dòng)量定理，對恒定力矩和周期性變化力矩下剛體的定軸轉(zhuǎn)動(dòng)運(yùn)動(dòng)狀態(tài)進(jìn)行了探究，并進(jìn)一步觀察了剛體做平面平行運(yùn)動(dòng)的軌跡. 探究了大、小角度下復(fù)擺的運(yùn)動(dòng)方程和周期差別. 有趣的是，通過對受迫阻尼振動(dòng)的模擬，檢驗(yàn)了非線性系統(tǒng)相關(guān)結(jié)論.

VPython不僅將課本知識形象化、可視化，加深了理解，還能快速解決學(xué)生學(xué)習(xí)中產(chǎn)生的疑問，激發(fā)學(xué)生對問題更深入地探索與思考. 在計(jì)算機(jī)技術(shù)迅猛發(fā)展，各種3D建模、數(shù)值分析軟件“百花齊放”的今天，培養(yǎng)學(xué)生能夠熟練運(yùn)用功能強(qiáng)大的軟件輔助學(xué)習(xí)和研究,對深化大學(xué)物理課堂教學(xué)改革具有一定啟示.