1、剛體運(yùn)動(dòng)三維可視化模擬和復(fù)擺相關(guān)問題探討摘要:本文應(yīng)用剛體定軸轉(zhuǎn)動(dòng)定律和角動(dòng)量定理，結(jié)合VPython軟件，首先模擬剛體的定軸轉(zhuǎn)動(dòng)，例如勻質(zhì)圓柱體在不同 力矩作用下的轉(zhuǎn)動(dòng)，直觀地顯示剛體運(yùn)動(dòng)軌跡以及角位移、角速度等物理量隨時(shí)間的變化關(guān)系，以及重力矩作用下的剛體能 量守恒等，使剛體運(yùn)動(dòng)過程實(shí)現(xiàn)可視化.有趣的是,本文還研究了剛體的平面平行運(yùn)動(dòng)（平動(dòng)加轉(zhuǎn)動(dòng)）以及剛體的振動(dòng),模擬了 大角度下剛體擺動(dòng)時(shí)出現(xiàn)的混沌狀態(tài).可視化模擬深化了對(duì)課本知識(shí)的理解.關(guān)鍵詞:力矩;VPython;可視化;轉(zhuǎn)動(dòng)定律;角動(dòng)量定理自17世紀(jì)牛頓經(jīng)典力學(xué)建立以來，就一直得到 普遍的應(yīng)用，它打開了近代自然科學(xué)的大門,在其基 礎(chǔ)上

2、發(fā)展出了彈性力學(xué)、流體力學(xué)、結(jié)構(gòu)力學(xué)等多門 學(xué)科.在牛頓力學(xué)中，如果已知物體的受力和初始 狀態(tài)，對(duì)于有序系統(tǒng)其未來的運(yùn)動(dòng)狀態(tài)是確定的. 學(xué)者們根據(jù)經(jīng)典力學(xué)的定律和萬有引力定律曾經(jīng)精 確地預(yù)言彗星和小行星等的運(yùn)動(dòng)，并且得到了驗(yàn)證.近年來研究者們?cè)诖髮W(xué)物理教學(xué)中引入Matlab 計(jì)算機(jī)數(shù)值模擬，培養(yǎng)學(xué)生用計(jì)算機(jī)進(jìn)行學(xué)習(xí)和研 究物理問題的能力3，將VPython引進(jìn)大學(xué)物理課 堂的方式也逐漸興起VPython可以利用簡(jiǎn)單的 代碼模擬物體在三維空間中的運(yùn)動(dòng)，實(shí)時(shí)觀測(cè)物體 的速度變化及其運(yùn)動(dòng)軌跡.質(zhì)點(diǎn)平動(dòng)、剛體定軸轉(zhuǎn) 動(dòng)和復(fù)擺振動(dòng)是基本的運(yùn)動(dòng)形式，本文選用 VPython軟件對(duì)剛體的平動(dòng)、轉(zhuǎn)動(dòng)、振動(dòng)進(jìn)行

3、模擬與 探究,利用轉(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模塊,有著和Py- thon 一樣風(fēng)格的編程方式，可以快速創(chuàng)建三維場(chǎng)景 和動(dòng)畫.Vpython可創(chuàng)建多種實(shí)體類型，如球體、圓 柱體、長(zhǎng)方體、錐體等,通過代碼改變其在空間中的 位置即可實(shí)現(xiàn)運(yùn)動(dòng)模擬.利用VPython模擬剛體轉(zhuǎn) 動(dòng)的代碼框架為:定義常數(shù)（&,g,（等）和運(yùn)動(dòng)初 始值（初始角

4、位置，初始角動(dòng)量等）;設(shè)置步長(zhǎng) 越小模擬越精確；進(jìn)入while循環(huán),在循環(huán) 中根據(jù)當(dāng)前狀態(tài)計(jì)算合外力矩;用該力矩更新系 統(tǒng)角動(dòng)量;利用角動(dòng)量更新系統(tǒng)角位置.程序中 以差分代替微分,得到角速度、角位置對(duì)時(shí)間的迭代 式（Euler-Cromer 方法）：!+*= 偵仇,!* ,t.） #）,知=t+#t#t越小，模擬的精度越高.2剛體轉(zhuǎn)動(dòng)的VPython模擬2.1轉(zhuǎn)動(dòng)定律簡(jiǎn)介在大學(xué)物理的學(xué)習(xí)中，我們依據(jù)牛頓第二定律 導(dǎo)出了剛體定軸轉(zhuǎn)動(dòng)定律M = 即剛體所受的對(duì) 于某定軸的合外力矩等于剛體對(duì)此定軸的轉(zhuǎn)動(dòng)慣量 與剛體在此合外力矩作用下所獲得的角加速度的乘 積;力矩的時(shí)間累積效應(yīng)一角動(dòng)量定理d = Md

5、t, 即剛體角動(dòng)量的增量等于所受合外力矩的角沖量.2.2恒定力矩下勻質(zhì)圓柱體的定軸轉(zhuǎn)動(dòng)有一個(gè)質(zhì)量為2 kg,長(zhǎng)度為1 m,半徑為0.1 m 的勻質(zhì)圓柱體，轉(zhuǎn)軸為過其質(zhì)心且與圓柱體中心軸 線垂直.勻質(zhì)圓柱體最初處于靜止?fàn)顟B(tài)，對(duì)其施加 一恒定力矩，力矩方向沿轉(zhuǎn)軸正方向，大小為2 Nm, 對(duì)其運(yùn)動(dòng)進(jìn)行分析.對(duì)于勻質(zhì)圓柱體，轉(zhuǎn)軸(圖中所示-軸)位于圓 柱體中心且垂直于中心軸線，建立如圖1所示的坐 標(biāo)系.$圖1勻質(zhì)圓柱體，其轉(zhuǎn)軸為-軸將長(zhǎng)為.，半徑為r,質(zhì)量為&,密度為p的勻質(zhì) 圓柱體分割為無數(shù)個(gè)垂直于軸、厚為心的薄圓2 1盤.已知薄圓盤對(duì)其直徑軸的轉(zhuǎn)動(dòng)慣量J,=.由 平行軸定理，距離原點(diǎn)處的薄圓盤對(duì)-

6、軸轉(zhuǎn)動(dòng)慣2 1量為dJ=、+2d&，積分得該圓柱體對(duì)-軸的轉(zhuǎn) 動(dòng)慣量為.2J- / $d/ Rpr2(: + 2) d2/#2 L+LCpr D土 $) =4 12/&r2 &L21412(1)當(dāng)沿轉(zhuǎn)軸正向的恒定力矩作用于勻質(zhì)圓柱體 時(shí)，圓柱體從靜止開始轉(zhuǎn)動(dòng).依據(jù)剛體的角動(dòng)量定 理，對(duì)于)時(shí)刻圓柱角動(dòng)量有廣L = Mt = J 在一段很小的時(shí)間間隔#t內(nèi),其運(yùn)動(dòng)可看作勻速轉(zhuǎn) 動(dòng),角位移故t時(shí)刻圓柱體角位置(1)令 #t =dt,得角位移 / $ ! dt / $ dt.將具體數(shù)值代入該公式，理論計(jì)算得t = 2 s時(shí),！23.30 rad、,%23.30 rad.S我們?cè)赩Python中對(duì)該勻

7、質(zhì)圓柱體的轉(zhuǎn)動(dòng)進(jìn)行 模擬，利用YPythoJ顯示勻質(zhì)圓柱體尺寸、顏色以及 轉(zhuǎn)軸位置(如圖2所示)，將棒的運(yùn)動(dòng)可視化.圖2勻質(zhì)圓柱體繞定軸轉(zhuǎn)動(dòng)及對(duì)應(yīng)坐標(biāo)系,轉(zhuǎn)軸沿-軸正向 在VPython模擬中，可以直觀地觀察它的運(yùn)動(dòng) 過程，程序中使用的物理公式及對(duì)應(yīng)的程序循環(huán)語 句如下.while t2:rate( 100)torque = vector( 0,0,2)L = L+torque* dt#M = ( #M = ( 0,0,2) idL = Mxdt,l L l# l! l =#d = l! lxdtomega = mag( L) /theta = theta+dthetarod.rotate(

8、angle = dtheta, axis = ( 0,0,1) , origin = axle,pos)t = t+dtVPython模擬結(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所示). 2.3正余弦變化的力矩下勻質(zhì)圓柱體的定軸轉(zhuǎn)動(dòng)對(duì)于m = 2 kg,L = 1 m,r

9、= 0.1 m的勻質(zhì)圓柱體, 當(dāng)轉(zhuǎn)軸位于質(zhì)心且與中心軸線垂直時(shí)，由2.1可得22mr mL其轉(zhuǎn)動(dòng)慣量J=亍%0.171667 kgm2,探究力矩 呈正余弦規(guī)律變化時(shí)($net = 3cos( 5t)，勻質(zhì)圓柱體 的運(yùn)動(dòng)狀態(tài)又如何？其運(yùn)動(dòng)周期等于力矩變化周期 2s嗎？轉(zhuǎn)動(dòng)動(dòng)能隨時(shí)間如何變化？所受力矩與其圖3角速度、角位置與時(shí)間關(guān)系圖像 轉(zhuǎn)動(dòng)動(dòng)能有何關(guān)系？力矩取極值時(shí)，圓柱體轉(zhuǎn)動(dòng)動(dòng) 能大小如何？dt%t 根據(jù)理論推導(dǎo)二#net二在)時(shí)亥叭d)積分可得角速度和角位置分別為1 !二一I 3cos( 5) d) % 3.495sin( 5),Jo()/ 0,!/ 0)( 2)= d)/ $ 3.495s

10、in( 5) d)二-0.699cos( 5) + 0.699,() / 0,/ 0) ( 3) 進(jìn)行VPython模擬，將力矩定義式更改為torque = 3cos( 5)模擬結(jié)果表明：角速度！、角位置隨時(shí)間 )均做周期性變化，角速度！、角位置隨時(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)系曲線力矩3和轉(zhuǎn)動(dòng)動(dòng)能E隨時(shí)

11、間變化曲線如圖6圖5勻質(zhì)圓柱體轉(zhuǎn)動(dòng)的3D場(chǎng)景所示圖5勻質(zhì)圓柱體轉(zhuǎn)動(dòng)的3D場(chǎng)景12132M $0 7d)% 0.524( 1 - cos10)力矩大時(shí)，轉(zhuǎn)動(dòng)動(dòng)能不一定大.通過動(dòng)能與力矩的函 數(shù)關(guān)系圖像(圖7)可知：當(dāng)力矩為極值3Nm時(shí)% 0.524( 1 - cos10)圖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)量定 理P*=得到動(dòng)量迭代式9 = P+Fdt,由9 =求得任意時(shí)刻的瞬時(shí)平動(dòng)速度；和dt內(nèi)的位移 ds = vdt.在while循環(huán)中更新轉(zhuǎn)軸的位置（設(shè)力沿！ 軸正方向，故只需改變勻質(zhì)圓柱體質(zhì)心！坐標(biāo)）.wh

12、ile t2:P = P+：* dt#dPgnslate 二d；=P/MP#% =Mrod.pos = rod.pos+;* dt#dx = ;dtaxle.pos.! = axle.pos.x+mag（ ;* dt）在相同的力矩和大小不同的外力作用下，系統(tǒng) 的運(yùn)動(dòng)軌跡不同，如圖8所示.力矩均為M =3cos5t, 外力大小分別為0.1 N、0.5 N、1 N沿！軸正向（方向 向右），模擬結(jié)果表明:勻質(zhì)圓柱體端點(diǎn)的軌跡呈鋸 齒狀,且隨著外力的增大，“鋸齒”越來越疏松（如圖圖8相同的力矩和不同外力作用下勻質(zhì)圓柱體的 運(yùn)動(dòng)軌跡，力矩均為M = 3cos5t3 勻質(zhì)棒的擺動(dòng)（大角度、小角度）的VPy

13、thon 模擬3.1重力矩下的棒的擺動(dòng)（大角度、小角度）在大學(xué)物理振動(dòng)學(xué)章節(jié)中，我們了解了簡(jiǎn)諧振 動(dòng)是最簡(jiǎn)單、最基本的振動(dòng)，小角度復(fù)擺即是簡(jiǎn)諧運(yùn) 動(dòng)的一個(gè)特例.在研究復(fù)擺的運(yùn)動(dòng)周期時(shí),令物體 作小振幅擺動(dòng)（擺角5。）.根據(jù)轉(zhuǎn)動(dòng)定律得d2 一&g=sin = J-（ 3）式中=為復(fù)擺質(zhì)心到轉(zhuǎn)軸的距離，當(dāng)較小時(shí)，可近 似看作sin%,式（3）變?yōu)? = 0,解微分方程得=7cos（ !t+%）振動(dòng)的圓頻率和周期分別為T =2mgl(4T =2mgl對(duì)于大角度下復(fù)擺的擺動(dòng)情況,大多數(shù)物理教 材沒有詳細(xì)提及.我們常常會(huì)思考：當(dāng)擺動(dòng)角5。時(shí)，物體還做簡(jiǎn)諧振動(dòng)嗎？它的運(yùn)動(dòng)周期又變?yōu)槎?少呢？因?yàn)楫?dāng)較大時(shí)，

14、擺動(dòng)微分方程理+!）sin=0為一個(gè)非線性微分方程.對(duì)于非線性系統(tǒng),非線性 運(yùn)動(dòng)方程不再滿足疊加原理，難以求出精確解析解; 初始條件不同會(huì)導(dǎo)致不同的運(yùn)動(dòng)形式;還可能會(huì)出 現(xiàn)隨機(jī)的混沌運(yùn)動(dòng).有研究表明可以應(yīng)用迭代法求 得方程三次迭代近似解為（o為初始偏置角）3=osin!+1120 ( 1+b540 ) sin3!+4Bo?o0sin5!(5)通過積分理論計(jì)算，得大擺角運(yùn)動(dòng)時(shí)復(fù)擺的周 期為回Qlsin2 -sin2。=2（B）=cos( !t),/3 - 9.8/2 rad/s%3.834 rad/s可以看出任意角度下擺的運(yùn)動(dòng)方程與運(yùn)動(dòng)周期 的理論推導(dǎo)式都與小角度近似（sin%）存在差異, 0為

15、大角度時(shí)，擺的運(yùn)動(dòng)已非簡(jiǎn)諧運(yùn)動(dòng).在VPython 中，對(duì)一根懸掛在低摩擦軸上的桿的運(yùn)動(dòng)（忽略阻 力）進(jìn)行建模（圖9=cos( !t),/3 - 9.8/2 rad/s%3.834 rad/s實(shí)線為擺桿實(shí)際角位移曲線.當(dāng)為小角度時(shí) (10?時(shí)，隨著初始角的增大，實(shí)際角位移曲 線與簡(jiǎn)諧運(yùn)動(dòng)角位移理論計(jì)算結(jié)果曲線偏差越大, 擺桿不再作簡(jiǎn)諧運(yùn)動(dòng)(如圖10( c) -( f)所示);而當(dāng) o = 18O。時(shí)，擺桿達(dá)到往復(fù)擺動(dòng)與圓周運(yùn)動(dòng)臨界點(diǎn).圖9擺桿Vpython建模實(shí)體(e)0o=12O(b)eqo(d)%=90。(f)0o=17O(a)R=5。(c)%=20。圖10不同初始偏轉(zhuǎn)角度下擺桿實(shí)際角位移曲

16、線與簡(jiǎn)諧 運(yùn)動(dòng)角位移理論曲線對(duì)比圖由圖11可知，當(dāng)較小時(shí)，相圖軌線為圓形閉 合曲線，表明運(yùn)動(dòng)為簡(jiǎn)諧振動(dòng)；o越大，軌線越 “方” .3.2 考慮阻尼和周期性外力矩時(shí)擺運(yùn)動(dòng)特性的 研究在經(jīng)典力學(xué)中，(e)0o=12O(b)eqo(d)%=90。(f)0o=17O(a)R=5。(c)%=20。圖10不同初始偏轉(zhuǎn)角度下擺桿實(shí)際角位移曲線與簡(jiǎn)諧 運(yùn)動(dòng)角位移理論曲線對(duì)比圖由圖11可知，當(dāng)較小時(shí)，相圖軌線為圓形閉 合曲線，表明運(yùn)動(dòng)為簡(jiǎn)諧振動(dòng)；o越大，軌線越 “方” .3.2 考慮阻尼和周期性外力矩時(shí)擺運(yùn)動(dòng)特性的 研究在經(jīng)典力學(xué)中，已知物體初始狀態(tài)和物體運(yùn)動(dòng) 過程中受力情況，可以通過牛頓運(yùn)動(dòng)定律和動(dòng)量定 理預(yù)

17、測(cè)物體的運(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ù)料，且對(duì)初始條件具有敏感性，這在物理上 被稱為混沌混沌的一種著名表述是“蝴蝶效應(yīng)”， 初始條件的極小偏差，都將可能會(huì)引起結(jié)果的極大 差異.受迫振動(dòng)微分方程為-&g=sin -!+3cos!p)(7)c為阻力系數(shù),3為力矩最大值,！，為力矩周期變化 的角頻率，上式化簡(jiǎn)為一Bsin+Cos !pt( 8)7 = -，B = g，/ =-在 VPython 中進(jìn)行模

18、擬，取7 =圖11不同初始條件下擺桿受重力矩的運(yùn)動(dòng)狀態(tài)！-相圖(a)戶0.5時(shí)，軌跡穩(wěn)定且收斂0.5,B =l,%=2/3,初始角(a)戶0.5時(shí)，軌跡穩(wěn)定且收斂0.5,B =l,%=2/3,初始角 = 1.6 rad % 91.673。，/從0.5開始取一系列值，觀察運(yùn)動(dòng)相圖如圖12所示.(d)戶1.115時(shí)，出現(xiàn)兩個(gè)疊加極限環(huán)(e)戶1.2時(shí)，軌跡混亂無序(f)戶1.38,相軌跡具有周期性截取不同初始條件下擺桿受迫阻尼振動(dòng)的！ 0相圖由于疊加原理、傅立葉變換等工具不適于求解 非線性方程，而計(jì)算機(jī)性能的提高使得研究者一般 采用計(jì)算機(jī)進(jìn)行迭代運(yùn)算對(duì)其進(jìn)行數(shù)值求解.通過 有阻尼受迫振動(dòng)的模擬，發(fā)現(xiàn)在一定范圍內(nèi)即使初 始條件/變化很小,卻會(huì)引起復(fù)擺的運(yùn)動(dòng)形式有很 大的不同.例如當(dāng)/ =1.505時(shí),相軌跡為無規(guī)則的發(fā) 散曲線;/ =1.508時(shí),相軌跡為收斂的三個(gè)極限環(huán).與 當(dāng)年洛倫茲在利用計(jì)算機(jī)研究地球大氣時(shí)提出的 “蝴蝶效應(yīng)”相類似,對(duì)于十分接近的兩個(gè)/的值，初 始時(shí)相軌跡重疊，當(dāng)時(shí)間積累到一定程度后，相軌跡 會(huì)出現(xiàn)較大偏差:設(shè)定模擬精度d) = 0.001 s,f初始 值相差0.00001時(shí)，程序運(yùn)行200 s后，兩個(gè)復(fù)擺的 運(yùn)動(dòng)將截然不同.隨著/值的變化，相軌跡