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

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

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

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

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

{\bar{v}}_{12} = \frac{x_{2} - x_{1}}{d t}

$\bar{v}_{12} = \dfrac{x_2-x_1}{dt}$ 位置2到位置3的平均速度為：

{\bar{v}}_{23} = \frac{x_{3} - x_{2}}{d t}

$\bar{v}_{23} = \dfrac{x_3-x_2}{dt}$

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

\bar{a} = \frac{v_{23} - v_{12}}{d t} = \frac{(x_{3} - x_{2}) - (x_{2} - x_{1})}{d t^{2}}

$\bar{a} = \dfrac{v_{23}-v_{12}}{dt} = \dfrac{\left ( x_3-x_2 \right ) - \left ( x_2-x_1 \right )}{dt^{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}$ ，則：

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

{\vec{v}}_{12} = \frac{Δ {\vec{r}}_{21}}{d t} = \frac{{\vec{r}}_{2} - {\vec{r}}_{1}}{d t}

$\vec{v}_{12} = \dfrac{\Delta \vec{r}_{21}}{dt} = \dfrac{\vec{r}_2-\vec{r}_1}{dt}$ 前一時刻到現在時刻的平均速度：

{\vec{v}}_{23} = \frac{Δ {\vec{r}}_{32}}{d t} = \frac{{\vec{r}}_{3} - {\vec{r}}_{2}}{d t}

$\vec{v}_{23} = \dfrac{\Delta \vec{r}_{32}}{dt} = \dfrac{\vec{r}_3-\vec{r}_2}{dt}$ 如果程式中模擬的時間

d t

$dt$ 足夠小，則

{\vec{v}}_{12}

$\vec{v}_{12}$ 和

{\vec{v}}_{23}

$\vec{v}_{23}$ 都可以視為現在時刻的瞬時速度。而這段期間的平均加速度為：

\vec{a} = \frac{{\vec{v}}_{23} - {\vec{v}}_{12}}{d t} = \frac{({\vec{r}}_{3} - {\vec{r}}_{2}) - ({\vec{r}}_{2} - {\vec{r}}_{1})}{d t^{2}}

$\vec{a} = \dfrac{\vec{v}_{23}-\vec{v}_{12}}{dt} = \dfrac{\left ( \vec{r}_3-\vec{r}_2 \right ) - \left ( \vec{r}_2-\vec{r}_1 \right )}{dt^{2}}$ 也因為

d t

$dt$ 極小，我們只要將這個式子直接寫到程式中就可以近似模擬出每瞬間的加速度，這套方法我們稱為「微分法」。而且你會發現這個加速度非常特別，其方向永遠指向圓心，稱為「向心加速度」！以下我們示範如何計算瞬時速度，請參考範例程式：

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'的位置

用 label 標示物理量名稱的設定技巧：

先在 while 迴圈上面設定 label 指令：
```
v_text = label(box = False, opacity = 0, height = 25, color=color.green, text = 'v')
```
接著在 while 迴圈裡面設定 label 的位置：
```
v_text.pos = v_vector.pos + v_vector.axis*1.2
```
label 指令中的 box 是用來設定文字是否需要外框，False 為無框、True 為有框；opacity 是設定文字標籤的透明度，數值設定在 0~1 之間，0 表示完全透明、1 表示完全不透明；height為文字大小；text 為顯示在執行畫面上的文字，需用單引號設定為字串。

習慣上直接將文字設定在箭頭的尖端處，畫面效果就相當精美。上面的 1.2 數字是為了讓文字離箭頭尖端處有一點點距離，否則文字會和箭頭的尖端重疊。

最後你會發現微分法的優點是，只要我知道物體隨時間怎麼動，就可以反推出它每一個時刻的速度與受力，透過很表面的現象就能反推背後的力學機制。即使我不知道物體隨時間的運動函數，也可以直接用攝影機錄下來，再利用影像分析軟體，例如 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)    #回傳虎克定律計算的彈力向量

用函數撰寫虎克定律說明：
.

r向量：由彈簧的頭端指向尾端的距離向量。
mag( r )：mag用來計算向量的量值，可透過mag( r )算出 $\vec{r}$ 的量值。
r / mag( r )：這是一個「單位向量」，方向沿著 $\vec{r}$ 方向，長度為 1 的向量。在這裡只用來指出彈力方向，不會影響彈力量值。

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

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表示物理量的代號，就像是在電腦裡直接畫力圖一樣，執行畫面如下所示。

提示：

均勻重力場的重力位能公式為 $U_{g}=mgh$ ，位能零位面請自訂，寫法為：
```
Ug = m * g * ball.pos.y。
```
彈力位能公式為 $\dfrac{1}{2}kx^{2}$ ，寫法為：
```
Us = 0.5 * k * (mag(spring.axis) - L)**2
```
質點的動能 $\dfrac{1}{2}mv^{2}$ ，寫法為：
```
K = 0.5 * m * mag(ball.v)**2
```

arrow 設定技巧：

先在 while 迴圈上面設定 arrow 指令繪圖：

v_vector = arrow(shaftwidth = 0.02, color=color.green)
Ftot_vector = arrow(shaftwidth = 0.02, color=color.red)
mg_vector = arrow(shaftwidth = 0.02, color=color.yellow)
Fs_vector = arrow(shaftwidth = 0.02, color=color.white)
X = 0.2

這裡的 X 是用來縮放調整箭頭的長度，以較美觀的比例呈現圖形。

接著請在 while 迴圈裡面設定箭頭的位置與軸方向：

#設定箭頭的位置，要將圖形中會重疊的箭頭微微錯開，會比較好觀察
v_vector.pos = ball.pos + vector(2.5*size,0,0)      #將速度箭頭設定在球的右方2.5size處
Ftot_vector.pos = ball.pos + vector(-2.5*size,0,0)  #將合力箭頭設定在球的左方2.5size處
mg_vector.pos = ball.pos    #將重力箭頭設定在球上
Fs_vector.pos = ball.pos    #將彈力箭頭設定在球上

#設定箭頭的軸向量
v_vector.axis = ball.v * X                                #設定速度箭頭，將軸方向設定為速度向量
Ftot_vector.axis = (Fg + SpringForce(spring.axis,L)) * X  #設定合力箭頭，將軸方向設定為合力向量
mg_vector.axis = Fg * X                                   #設定重力箭頭，將軸方向設定為重力向量
Fs_vector.axis = SpringForce(spring.axis,L) * X           #設定彈力箭頭，將軸方向設定為彈力向量

箭頭的長度的縮放，請務必留意相同單位的物理量縮放的比例要一致，例如：這裡的重力、彈力、合力單位都是牛頓，縮放比例必須一致，這樣在三角形法、平行四邊形法的向量合成比例才會是正確的。

label 設定技巧：

先在 while 迴圈上面設定 label 指令繪圖：

v_text = label(box = False, opacity = 0, height = 25, color=color.green, text = 'v')
Fs_text = label(box = False, opacity = 0, height = 25, color=color.white, text = 'Fs')
mg_text = label(box = False, opacity = 0, height = 25, color=color.yellow, text = 'Fg')
Ftot_text = label(box = False, opacity = 0, height = 25, color=color.red, text = 'F_tot')

接著在 while 迴圈裡面設定文字的位置：

v_text.pos = v_vector.pos + v_vector.axis*1.2
Ftot_text.pos = Ftot_vector.pos + Ftot_vector.axis*1.2
mg_text.pos = mg_vector.pos + mg_vector.axis*1.2
Fs_text.pos = Fs_vector.pos + Fs_vector.axis*1.2

三、單擺運動

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

如上圖，假設初始擺角為 $\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 指令標示上去，但你會發現箭頭會瘋狂的抖動！原因是我們這裡用了非常高剛性的彈簧，因此會造成高頻率的微小幅度振動，不畫箭頭看不出來，畫了箭頭後整個振動的效果就相當明顯。

高剛性彈簧的高頻振動處理：

這裡我們沿著彈簧的軸方向設計一個「避震器」，讓彈簧沿軸方向的振動迅速被消弭掉，整體的受力情況就會逐漸趨於穩定狀態。有關避震器的設定如下：
```
def SpringDamp(v, r):  #避震器
    cos_theta = dot(v,r)/(mag(v)*mag(r))                    #用向量內積找v和r夾角的餘弦函數
    r_unit_vector = norm(r)                                 #沿彈簧軸方向的單位向量
    projection_vector = mag(v) * cos_theta * r_unit_vector  #計算v在r方向的分量
    spring_damp = - damp * projection_vector                #沿彈簧軸方向的阻力
    return spring_damp
```
這裡暫時不詳細解釋這段程式碼的細節，同學只要記得計算物體合力的時候，要把避震器加上去就好了。

避震器的實際結構大致如下圖：主要有兩部份，一是彈簧、二是阻尼器，避震的功能主要是阻尼器的貢獻，阻尼器壓縮時有阻力、伸長時也有阻力，故能迅速消除彈簧的振動。

圖片來源：http://www.carnews.com/article/info/d9927c78-4b11-11e8-8ee2-42010af00004/

以上是由靜止狀態釋放的單擺運動，如果一開始就給擺錘初速度的話，而且速度量值剛好為 $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

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

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

提示：假設兩根棒子一樣長，且初始狀態為水平，由靜止狀態釋放

一開始紅球的位置在ball_1.pos = vector(L,0,0)，
綠球的位置在ball_2.pos = vector(2*L,0,0)。
內棒的初始位置為rod_1.pos = vector(0,0,0)，
外棒的初始位置為rod_2.pos = vector(L,0,0)。

以上設定好初始條件之後，在迴圈裡計算每一刻兩根棒子的軸方向

while True:
...
    rod_2.pos = ball_1.pos                 #外棒的位子在紅球處
    rod_1.axis = ball_1.pos                #內棒的軸方向由原點指向紅球
    rod_2.axis = ball_2.pos - ball_1.pos   #外棒的軸方向由紅球指向綠球
    ...

棒子的軸方向決定之後，就可以代入彈力的函數中，直接算出這兩跟棒子的彈力，而紅球和綠球即是受彈力和重力而運動。

紅球所受的合力為：

F1 = vector(0, -m*g, 0) + SpringForce(rod_1.axis) - SpringForce(rod_2.axis)

綠球所受的合力為：

F2 = vector(0, -m*g, 0) + SpringForce(rod_2.axis)

知道這兩個球所受的作用力，即可找出每瞬間的加速度，加速度使速度改變，速度使位置改變，即可完成整個動力學模擬。
```
ball_1.v += F1/m*dt
ball_1.pos += ball_1.v*dt
ball_2.v += F2/m*dt
ball_2.pos += ball_2.v*dt
```