Lecture 3 動力學實例模擬－圓周運動、振動與擺動

建國高中特色選修課程 - 物理現象的程式設計與模擬
作者：曾靖夫
日期：2018/7/5

一、等速率圓周運動的力學分析 - 微分法

上一堂課的最後，學到了如何利用三角函數在程式中畫出圓周運動，並且紀錄了三個數據點，我們簡單稱為「三點記錄法」。在物理學中，只要三個點我們就可計算出平均速度與平均加速度，如果時間間隔取的夠短，就可以擬真呈現出瞬時速度與瞬時加速度。以下我們簡單利用一維等加速度直線運動進行概念說明：

上方左圖為一物體在直線上做等加速度運動以相同時間間隔記錄下來的三個點，假設經過這三點的瞬間速度分別為 $v_1$、$v_2$、$v_3$ 標示在旁邊，上方右圖為物體運動過程的 $v-t$ 圖。由位置1到位置2的平均速度為：

而我們又可以利用上方前半段和後半段的速度差計算出平均加速度為：

注意這裡雖然 $\bar{v}_{12}$、$\bar{v}_{23}$ 是在講平均速度，但是在 $v-t$ 中你可以看到 $\bar{v}_{12}$ 指的也是質點1到質點2時間中點的瞬時速度，同理 $\bar{v}_{23}$ 是質點2到質點3時間中點的瞬時速度。

以上雖然只是在分析直線運動，但是其觀念運用到立體的空間向量也依然適用，以下我們就簡單分析如何利用圓周運動過程中持續在紀錄的三個數據點計算出每一時刻的速度、加速度及作用力。下圖中我們畫出被紀錄的三個位置，其位置向量分別為 $\vec{r}_{1}$、$\vec{r}_{2}$、$\vec{r}_{3}$，則：

前前時刻到前一時刻的平均速度：

Example 1 : 等速率圓周運動的力學分析

from vpython import *

size = 0.1
theta = 0.0
R = 1.0
omega = 2*pi

scene = canvas(width=500, height=500, center=vector(0,0,0), range=1.5*R, background=vector(148.0/255,228.0/255,204.0/255))
ball = sphere(radius=size, color=color.blue, make_trail=True, interval=1, retain=1000)
ball.pos = vector(R,0,0)

v_vector = arrow(color=color.green)  #畫描述v的箭頭
v_text = label(box = False, opacity = 0, height = 25, color=color.green, text = 'v')  #產生文字標籤'v'

dt = 0.001
t = 0.0

pre_theta = theta

while True:
    rate(1/dt)
    t = t + dt

    pre_pre_theta = pre_theta    #前前時刻角度
    pre_theta = theta            #前一時刻角度
    theta = theta + omega*dt     #現在時刻角度

    pre_pre_ball_pos = vector(R*cos(pre_pre_theta),R*sin(pre_pre_theta),0) #球前前時刻的位置
    pre_ball_pos = vector(R*cos(pre_theta),R*sin(pre_theta),0)             #球前一時刻的位置
    now_ball_pos = vector(R*cos(theta),R*sin(theta),0)                     #球現在時刻的位置
    ball.pos = pre_ball_pos    #將球物件畫在前一時刻的位置（ps.上面三個可任挑，dt很小肉眼看不出來差別）
    
    ball_v_12 = (pre_ball_pos - pre_pre_ball_pos)/dt   #前半程平均速度
    ball_v_23 = (now_ball_pos - pre_ball_pos)/dt       #後半程平均速度
    ball.v = ball_v_12         #將球的速度指定為前半程平均速度（ps.這兩個平均速度也可任挑，理由同上）
    
    v_vector.pos = ball.pos    #將速度箭頭的位置放在球的位置
    if mag(ball.v) > 0: v_vector.axis = ball.v/mag(ball.v)*R/2    #將速度箭頭的軸方向指向速度方向，並將箭頭的長度設定為半徑的一半
    
    v_text.pos = v_vector.pos + v_vector.axis*1.2 #文字標籤'v'的位置

最後你會發現微分法的優點是，只要我知道物體隨時間怎麼動，就可以反推出它每一個時刻的速度與受力，透過很表面的現象就能反推背後的力學機制。即使我不知道物體隨時間的運動函數，也可以直接用攝影機錄下來，再利用影像分析軟體，例如 tracker 截下每一個時間點的物體所在位置座標，再將數據匯入後一樣可以算出物體的受力情況。

課堂作業 3-1

請用微分法完成圓周運動速度與加速度向量模擬，執行結果如下附圖，綠色箭頭是速度向量、紅色箭頭是加速度向量，同時在箭頭尖端處以 label 指令產生文字描述物理量名稱。最後對照此程式模擬出來的加速度量值 mag(ball.a) 與向心加速度公式 $a_C = R\omega^{2} = \dfrac{v^{2}}{R}$ 計算出來的結果是吻合的。

二、彈力 - 虎克定律在程式中的寫法

在國中同學們學到的虎克定律為彈力大小與彈簧的伸長量成正比，但現在還必須要注意彈力方向，而且這裡指的是彈簧對物體的作用力。虎克定律向量形式為：$\vec{F}_{s}=-k\vec{x}$，這裡的 $k$ 稱為彈力常數，一般彈簧在彈性限度內 $k$ 為一定值，其SI單位為$\rm{N/m}$。

在程式中彈力的寫法可以引用python內建「函數」的寫法，其指令如下：

def functionname( parameters ):	    #Python function官方指令
    "function_docstring"			
    function_suite
    return [expression]

彈力函數：

def SpringForce(r, L):                    #建立彈力函數
    return - k*(mag(r)-L) * r / mag(r)    #回傳虎克定律計算的彈力向量

接著我們將模擬一個物體懸掛在一鉛錘懸掛的彈簧上振動，這時候物體每瞬間同時都受到兩個力量作用－重力與彈力，而這兩力的合力使物體在彈簧上來回振動。你可以任意修改彈力常數，改變物體的質量與釋放點，從中觀察各種參數造成的影響，看看是否吻合你的物理直覺。以下請參考彈簧振動的範例程式：

Example 2 : 鉛直彈簧振動模擬

from vpython import *
g = 9.8                 #重力加速度 9.8 m/s^2
size = 0.05             #球半徑 0.05 m            
L = 0.5                 #彈簧原長 0.5m
k = 10                  #彈簧力常數 10 N/m
m = 0.1                 #球質量 0.1 kg
Fg = m*vector(0,-g,0)   #球所受重力向量

def SpringForce(r,L):
    return -k*(mag(r)-L)*r/mag(r)

scene = canvas(width=600, height=600, center=vector(0, -L*0.8, 0), range=1.2*L)#設定畫面
ceiling = box(length=0.4, height=0.005, width=0.4, opacity = 0.2)#畫天花板
ball = sphere(radius = size,  color=color.yellow, make_trail = True, retain = 1000, interval=1)#畫球
spring = helix(radius=0.02, thickness =0.01)#畫彈簧
	
ball.pos = vector(0, -L, 0)     #球在t=0時的位置
ball.v = vector(0, 0, 0)        #球初速
spring.pos = vector(0, 0, 0)    #彈簧頭端的位置

dt = 0.001

while True:
    rate(1/dt)
    spring.axis = ball.pos - spring.pos           #彈簧的軸方向，由頭端指向尾端的距離向量
    ball.a = (Fg + SpringForce(spring.axis,L))/m  #球的加速度由重力和彈力所造成
    ball.v += ball.a*dt      #球的速度
    ball.pos += ball.v*dt    #球的位置

進階作業 3-1

請分析彈簧系統的力學能，力學能＝動能+彈力位能+重力位能，各種位能的零位面請自訂，請利用graph指令畫圖驗證彈簧系統在沒有任何阻力作用下，力學能守恆。同時也請在執行視窗中將彈簧所受的作用力及其合力與物體運動速度以箭頭標示，並在箭頭尖端處用label表示物理量的代號，就像是在電腦裡直接畫力圖一樣，執行畫面如下所示。

三、單擺運動

在我們了解了彈力之後，也許你可以猜得到，如果一個彈簧的彈力常數超大，其實就可以看成一根不可伸縮的棒子，如果在桿子的尾端黏上一個球並由靜止狀態釋放，就會形成單擺運動。以下我們將以此概念模擬單擺：

如上圖，假設初始擺角為 $\theta$，棒子長度為 $L$，棒子頭端固定在原點 $\left ( 0,0,0 \right )$，尾端座標為 $\left ( L {\rm sin}\theta, -L {\rm cos}\theta,0 \right )$，在釋放後的每一瞬間，球只受到重力和棒子作用在球上的力。請參考以下範例程式：

Example 3 : 單擺運動模擬

from vpython import *
g = 9.8                 #重力加速度 9.8 m/s^2
size = 0.05             #球半徑 0.05 m            
L = 0.5                 #彈簧原長 0.5m
k = 100000              #彈簧力常數 100000 N/m
m = 0.1                 #球質量 0.1 kg
theta = 30 * pi/180
Fg = m*vector(0,-g,0)   #球所受重力向量

def SpringForce(r,L):
    return -k*(mag(r)-L)*r/mag(r)

scene = canvas(width=600, height=600, center=vector(0, -L*0.8, 0), range=1.2*L)#設定畫面
ceiling = box(length=0.4, height=0.005, width=0.4, opacity = 0.2)#畫天花板
ball = sphere(radius = size,  color=color.yellow, make_trail = True, retain = 1000, interval=1)#畫球
rod = cylinder(radius=size/10)#畫棒子
	
ball.pos = vector(L*sin(theta), -L*cos(theta), 0)   #球的初始位置
ball.v = vector(0, 0, 0)                            #球初速
rod.pos = vector(0, 0, 0)                           #棒子頭端的位置

dt = 0.001    #時間間隔
t = 0.0       #初始時間

pre_x = ball.pos.x         #三點記錄法，初始設定

while True:
    rate(1/dt)
    rod.axis = ball.pos - rod.pos                #桿子的軸方向：由桿子頭端指向尾端的向量
    ball.a = (Fg + SpringForce(rod.axis,L))/m    #牛頓第二定律：加速度＝合力/質量

    pre_pre_x = pre_x      #三點記錄法，前前時刻x座標
    pre_x = ball.pos.x     #三點記錄法，前一時刻x座標

    ball.v += ball.a*dt    #速度
    ball.pos += ball.v*dt  #三點記錄法，現在時刻x座標
    t += dt

    if pre_x > pre_pre_x and pre_x > ball.pos.x:    #計算單擺週期
        print ('simulated period = ', t, ', theoretical period = ', 2*pi*(L/g)**0.5 )    #印出單擺週期
        t = 0    #印出單擺週期

執行畫面如下：

在程式中我們使用三點記錄法來計算週期，正好可以驗證小角度單擺週期的理論公式 $T=2\pi\sqrt{\dfrac{L}{g}}$，是否與模擬出來的週期吻合。跑了一段時間後，執行視窗記錄下來的模擬週期與週期理論值的對照如下：

你會發現這兩組數據雖然接近但還是差了一些些！你有猜到為什麼會這樣嗎？

最後我們將重力、彈力、合力、速度這些向量用 arrow 指令標示上去，但你會發現箭頭會瘋狂的抖動！原因是我們這裡用了非常高剛性的彈簧，因此會造成高頻率的微小幅度振動，不畫箭頭看不出來，畫了箭頭後整個振動的效果就相當明顯。

以上是由靜止狀態釋放的單擺運動，如果一開始就給擺錘初速度的話，而且速度量值剛好為 $v_{0}=\sqrt{gL \ {\rm sin} \theta \ {\rm tan} \theta}$，方向沿 $z$ 方向，則擺錘將在同一水平面上形成錐動擺運動。這裡我們直接給出錐動擺條件的答案，這在高二的物理課會證明，同學們其實可以自己隨便設定初速，試試看會有什麼運動型態產生。下面是標準的錐動擺運動示意圖：

課堂作業 3-2

請由單擺運動程式碼更改擺錘初速度，使之形成錐動擺。並在執行視窗印出模擬週期與週期理論值 $T=2\pi\sqrt{\dfrac{L {\rm cos} \theta}{g}}$，驗證兩者數值相同。同時也將各作用力與速度利用箭頭標示在執行畫面上，執行畫面如下圖所示。另外你也可以更改彈力常數，觀察其運動型態有什麼不同。

進階作業 3-2

請模擬雙擺運動，將單擺程式再加入一個擺掛在下端，使之形成雙擺，並以力學能守恆驗證模擬的正確性。執行畫面如下圖所示：

右上的數據圖中，綠色線為綠球的力學能，紅色線為紅球的力學能，藍色線為總力學能，你會發現藍色線數值幾乎不變呈水平狀態。

四、Future more …

2016 IYPT #2

2018 IYPT #11

作業繳交：

平台連結：…

此平台由南港高中-高慧君老師設計建置