一個地球人在臺北: 圓周運動

作者：王一哲
日期：2018/4/4

如果一個小球的質量為 m，速度為 v，轉彎時的曲率半徑為 R，則轉彎需要的向心加速度與向心力分別為

$$ F_c = ma_c = m \cdot \frac{v^2}{R} $$

也就是說，如果想要用 VPython 畫出小球在水平面上做等速率圓周運動，我們需要想辦法計算向心加速度的大小與方向。如果成功地畫出水平面上的等速率圓周運動，也許就可以進一步挑戰鉛直面圓周運動。以下共有3個程式：

只畫出小球的等速率圓周運動。　（GlowScript 網站動畫連結）
畫出小球、繩子、轉軸的等速率圓周運動。　（GlowScript 網站動畫連結）
鉛直面圓周運動，畫出速率、切線加速度、法線加速度與時間的關係圖，計算週期。　（GlowScript 網站動畫連結）

程式 7-1：圓周運動

取得程式碼

"""
 VPython教學: 7-1.圓周運動
 Ver. 1: 2018/2/22
 Ver. 2: 2019/9/7
 作者: 王一哲
"""
from vpython import *

"""
 1. 參數設定, 設定變數及初始值
"""
size = 0.5            # 小球半徑
v0 = 10               # 小球初速
R = 5                 # 圓周運動半徑
L = 4*R               # 地板長度
t = 0                 # 時間
dt = 0.001            # 時間間隔

"""
 2. 畫面設定
"""
scene = canvas(title="Circle", width=800, height=400, x=0, y=0, background=vec(0, 0.6, 0.6))
scene.camera.pos = vec(0, L/2, L/2)
scene.camera.axis = vec(0, -L/2, -L/2)
floor = box(pos=vec(0, -size, 0), size=vec(L, 0.01, L), texture=textures.metal)
ball = sphere(pos=vec(R, size, 0), radius=size, color=color.red, make_trail=True, retain=100, v=vec(0, 0, -v0))
arrow_v = arrow(pos=ball.pos, axis=ball.v, radius=0.2*size, shaftwidth=0.4*size, color=color.green)
arrow_a = arrow(pos=ball.pos, axis=vec(0, 0, 0), radius=0.2*size, shaftwidth=0.4*size, color=color.blue)

"""
 3. 物體運動部分
"""
while(True):
    rate(1000)
    axis = vec(0, size, 0) - ball.pos
    ball.a = (ball.v.mag2 / R) * axis.norm()
    ball.v += ball.a*dt
    ball.pos += ball.v*dt
    arrow_v.pos = ball.pos
    arrow_v.axis = ball.v
    arrow_a.pos = ball.pos
    arrow_a.axis = ball.a
    t += dt

水平面圓周運動

參數設定

在此定義的變數有 size、v0、R、L、t、dt，用途都已經寫在該行的註解中。

畫面設定

這裡用到了一個新的功能，camera.pos 用來觀察者所在的位置，camera.axis 則是設定觀察者看畫面的方向。

scene.camera.pos = vector(0, L/2, L/2)
scene.camera.axis = vector(0, -L/2, -L/2)

小球是在 xz 平面上運動，出發點在畫面正右方距離 R 處，初速度方向朝 -z 軸、量值為 v0。
arrow_v 和 arrow_a 是用來表示小球速度和加速度的箭頭。

物體運動

用axis = vector(0, size, 0) - ball.pos找出小球相對於轉軸的位置向量。
用ball.a = (ball.v.mag2 / R) * axis.norm()計算小球的向心加速度。其中 mag2 是用來計算向量的量值平方，假設向量的名稱為 A ，語法有以下兩種
```
A.mag2 = mag2(A)
```
其中norm()是用來計算單位向量，假設向量的名稱為 A ，語法有以下兩種
```
A.norm() = norm(A)
```
最後更新小球的速度、位置，更新箭頭的起點位置、方向及長度，更新時間。

程式 7-2.圓周運動, 加上轉軸及繩子

取得程式碼

"""
 VPython教學: 7-2.圓周運動, 加上轉軸及繩子
 Ver. 1: 2018/2/22
 Ver. 2: 2019/9/7
 作者: 王一哲
"""
from vpython import *

"""
 1. 參數設定, 設定變數及初始值
"""
size = 0.5            # 小球半徑
v0 = 10               # 小球初速
R = 5                 # 圓周運動半徑
L = 4*R               # 地板長度
t = 0                 # 時間
dt = 0.001            # 時間間隔

"""
 2. 畫面設定
"""
scene = canvas(title="Circle with Rope", width=800, height=400, x=0, y=0, background=vec(0, 0.6, 0.6))
scene.camera.pos = vec(0, L/2, L/2)
scene.camera.axis = vec(0, -L/2, -L/2)
floor = box(pos=vec(0, -size, 0), size=vec(L, 0.01, L), texture=textures.metal)
ball = sphere(pos=vec(R, 0, 0), radius=size, color=color.red, make_trail=True, retain=100, v=vec(0, 0, -v0))
center = cylinder(pos=vec(0, -size, 0), axis=vec(0, 2*size, 0), radius=0.1*size, color=color.white)
rope = cylinder(pos=vec(0, 0, 0), axis=ball.pos, radius=0.1*size, color=color.yellow)
arrow_v = arrow(pos=ball.pos, axis=ball.v, radius=0.2*size, shaftwidth=0.4*size, color=color.green)
arrow_a = arrow(pos=ball.pos, axis=vec(0, 0, 0), radius=0.2*size, shaftwidth=0.4*size, color=color.blue)

"""
 3. 物體運動部分
"""
while True:
    rate(1000)
    axis = ball.pos
    ball.a = -(ball.v.mag2 / R) * axis.norm()
    ball.v += ball.a*dt
    ball.pos += ball.v*dt
    rope.axis = axis
    arrow_v.pos = ball.pos
    arrow_v.axis = ball.v
    arrow_a.pos = ball.pos
    arrow_a.axis = ball.a
    t += dt

水平圓周運動（加上繩子和轉軸）

水平圓周運動（繩子被表示加速度的箭頭擋住）

程式 7-2 和 7-1 非常像，只是增加了繩子和轉軸，新增的程式碼為27、28、41行

center = cylinder(pos = vector(0, -size, 0), axis = vector(0, 2*size, 0), radius = 0.1*size, color = color.white)
rope = cylinder(pos = vector(0, 0, 0), axis = ball.pos, radius = 0.1*size, color = color.yellow)
rope.axis = axis

分別用來畫出轉軸、畫出繩子、更新繩子的長度及方向。

程式7-3.鉛直面圓周運動

取得程式碼

"""
 VPython教學: 7-3.鉛直面圓周運動
 Ver. 1: 2018/3/30
 Ver. 2: 2019/9/7
 作者: 王一哲
"""
from vpython import *

"""
 1. 參數設定, 設定變數及初始值
"""
size = 0.5            # 小球半徑
R = 5                 # 圓周運動半徑
g = 9.8               # 重力加速度 9.8 m/s^2
v0 = 1*sqrt(g*R)      # 小球初速, 1 ~ 7 sqrt(g*R)
ratio = 0.1           # 速度, 加速度箭頭長度與實際的比例
i = 0                 # 小球回到出發點的次數
t = 0                 # 時間
dt = 0.0001           # 時間間隔, 取 0.0001 以降低誤差

"""
 2. 畫面設定
"""
scene = canvas(title="Vertical Circle", width=600, height=600, x=0, y=0, background=color.black)
ball = sphere(pos=vec(0, R, 0), radius=size, color=color.red, make_trail=True, retain=100, v=vec(-v0, 0, 0))
center = cylinder(pos=vec(0, 0, -size), axis=vec(0, 0, 2*size), radius=0.1*size, color=color.white)
rope = cylinder(pos=vec(0, 0, 0), axis=ball.pos, radius=0.1*size, color=color.yellow)
arrow_v = arrow(pos=ball.pos, axis=vec(0, 0, 0), radius=0.2*size, shaftwidth=0.4*size, color=color.green)
arrow_a = arrow(pos=ball.pos, axis=vec(0, 0, 0), radius=0.2*size, shaftwidth=0.4*size, color=color.blue)
gd = graph(title="plot", width=600, height=450, x=0, y=600, xtitle="<i>t</i>(s)",
           ytitle="green: <i>v</i>(m/s), red: <i>a<sub>t</sub></i>(m/s<sup>2</sup>), blue: <i>a<sub>n</sub></i>(m/s<sup>2</sup>)")
vt_plot = gcurve(graph=gd, color=color.green)
at_plot = gcurve(graph=gd, color=color.red)
an_plot = gcurve(graph=gd, color=color.blue)

"""
 3. 自訂函式, findan 計算法線加速度, findat 計算切線加速度
"""
def findan(v, pos):
    an = -v.mag2 / R * pos.norm()
    return an

def findat(pos):
    x = pos.x
    y = pos.y
    r = sqrt(x**2 + y**2)
    sintheta = abs(x)/r
    costheta = abs(y)/r
    absat = g*sintheta
    aty = -absat*sintheta
    if((x <= 0 and y <= 0) or (x >=0 and y>= 0)):
        atx = +absat*costheta
    elif((x <= 0 and y >= 0) or (x >= 0 and y <= 0)):
        atx = -absat*costheta
    at = vec(atx, aty, 0)
    return at

"""
 4. 物體運動部分, 小球回到出發點 5 次停止運作
"""
while i < 5:
# 由於 dt 較小，每秒計算 5000 次使動畫速度加快
    rate(5000)
# xp 是小球原來的位置, xc 是小球現在的位置, 用來判斷小球是否回到出發點    
# 計算小球 an, at, 更新加速度, 速度, 位置
    xp = ball.pos.x
    an = findan(ball.v, ball.pos)
    at = findat(ball.pos)
    ball.a = an + at
    ball.v += ball.a*dt
    ball.pos += ball.v*dt
    xc = ball.pos.x
    rope.axis = ball.pos
# 若小球回到出發點, 將 i 加 1, 印出時間 t, 由於誤差會累積, 取第一次回到出發點的時間作為週期
    if(xp > 0 and xc < 0):
        i += 1
        print(i, t)
# 更新代表速度, 加速度的箭頭
    arrow_v.pos = ball.pos
    arrow_v.axis = ball.v * ratio
    arrow_a.pos = ball.pos
    arrow_a.axis = ball.a * ratio
# 更新 v-t, at-t, an-t 圖
    vt_plot.plot(t, ball.v.mag)
    at_plot.plot(t, at.mag)
    an_plot.plot(t, an.mag)
# 更新時間
    t += dt

鉛直圓周運動

鉛直圓周運動速度、切線加速度、法線加速度與時間關係圖

理論分析

為了使計算過程較為簡單，我將圓心設定在原點，小球由 $(0, R, 0)$ 處以初速度 $(-v_0, 0, 0)$ 出發。加速度可以分為改變速度方向用的法線加速度 $a_n$，也就是向心加速度 $a_c$，以及改變速度量值用的切線加速度 $a_t$。$a_n$ 的計算方式與程式 7-1 相同，不再贅述。$a_t$的來源為重力加速度的切線方向分量，假設小球與鉛垂線的夾角為$\theta$，則

$$ a_t = g \sin \theta $$

$a_t$ 的 $y$ 方向分量方向下向，量值為

$$ a_{t,y} = a_t \sin \theta = g \sin^2 \theta $$

$a_t$ 的 $x$ 方向分量方向與位置有關，在第1和第3象限向右，在第2和第4象限向左，量值為

$$ a{t,x} = a_t \cos \theta = g \sin \theta \cos \theta $$

如果只想找出週期，可以改由力學能守恆計算。假設小球在最高點的速度量值為

$$ v_0 = n \sqrt{gR} $$

再由力學能守恆可寫出任意點與最高點的關係式 [5]

$$ \frac{1}{2}mv^2 + mgR \cos \theta = \frac{1}{2}mn^2 gR + mgR $$

$$ v = R \frac{d \theta}{dt} = \sqrt{(n^2 + 2 - 2\cos \theta)gR} $$

$$ T = \int_0^T dt = 2 \sqrt{\frac{R}{g}} \int_0^{2 \pi} \frac{d \theta}{\sqrt{n^2 + 2 - 2\cos \theta}} $$

參數設定

在此定義的變數有 size、R、g、v0、ratio、i、t、dt，用途都已經寫在該行的註解中。為了減少代入的時間長度造成的誤差，將 dt 的值調整為 0.0001。

畫面設定

小球是在 xy 平面上運動，出發點在畫面正上方距離 R 處，初速度方向朝 -x 軸、量值為 v0，不需要調整觀察者的位置及方向。
轉軸的方向改為指向 +z 軸方向。
arrow_v 和 arrow_a 是用來表示小球速度和加速度的箭頭。
用 graph 開啟繪圖視窗，用 gcurve 畫出 v - t、at - t、an - t 關係圖。

自訂函式

自訂函式 findan，輸入的參數為 v 和 pos，格式皆為向量，用來計算小球的法線加速度$ a_n$。
自訂函式 findat，輸入的參數為 pos，格式為向量，用來計算小球的切線加速度 $a_t$。繩子和鉛垂線的夾角為 $\theta$，先計算 $\sin \theta$ 和 $\cos \theta$ ，由於 $\theta$ 只取銳角，因此兩者皆為正值。再計算 $a_{t,y}$ 和 $a_{t,x}$，並由 pos 判斷 $a_{t,x}$ 的正、負號。

物體運動

由小球的位置判斷是否回到出發點。由於小球是由最高點開始逆時鐘方向旋轉，小球回到出發點的條件為：小球原來的位置 xp 在 +x 區域、現在的位置 xc 在 -x 區域。如果回到出發點，印出經過的時間 t、將次數 i 加 1，當小球回到出發點 5 次時停止動畫。
由於 dt 較小，將每秒計算的次數改為 5000，使動畫速度加快。
使用自訂函式 findan、findat 計算小球的法線加速度及切線加速度。
更新小球的速度、位置，更新箭頭的起點位置、方向及長度，更新 v - t、at - t、an - t 關係圖，更新時間。

數據處理部分

若小球在最高點的速度量值為 $v_0 = n \sqrt{gR}$，將 n 用 1 ~ 7 代入，分別由模擬程式及理論計算找出週期，數據如下

n	v0	T理論值	T模擬值
1	7	2.8841511713977500	2.8838000000016613
2	14	1.8728982530658000	1.8728999999998102
3	21	1.3617223964215400	1.3617999999999999
4	28	1.0602946238731300	1.0602999999998997
5	35	0.8646052540953560	0.8645999999999211
6	42	0.7284267507699600	0.7283999999999361
7	49	0.6286255969937600	0.6285999999999471

理論值與模擬值幾乎一樣。將數據匯入 SciDAVis 畫出 T - v0 關係圖如下

T - v₀ 關係圖

如果想要用 Python 計算數值積分，需要用到另一個套件 SymPy，程式碼如下：

from __future__ import division
from sympy import *

x = symbols('x', communtative = True)

def f(n, x):
    return 1/sqrt(n**2 + 2 - 2*cos(x))

for n in range(1, 8, 1):
    print("n = {:d} \tT = {:.16f}".format(n, integrate(f(n, x), (x, 0, pi)).evalf() * 2 * sqrt(5/9.8)))

結語

模擬鉛直圓周運動時，如果代入的 dt 較長，誤差就會很快累積並放大，接下來的動畫及數據就會與原來預設的條件相差很多，這是需要特別注意的地方。

參考資料

canvas: http://www.glowscript.org/docs/VPythonDocs/canvas.html
box: http://www.glowscript.org/docs/VPythonDocs/box.html
sphere: http://www.glowscript.org/docs/VPythonDocs/sphere.html
vector: http://www.glowscript.org/docs/VPythonDocs/vector.html
徐國誠(2013)。鉛直面圓周運動的週期。翰林自然科學天地，42，10-15。

HackMD 版本連結：https://hackmd.io/@yizhewang/HkyXcLfzm

一個地球人在臺北

熱門文章

2018年4月4日星期三

圓周運動

程式 7-1：圓周運動

參數設定

畫面設定

物體運動

程式 7-2.圓周運動, 加上轉軸及繩子

程式7-3.鉛直面圓周運動

理論分析

參數設定

畫面設定

自訂函式

物體運動

數據處理部分

結語

參考資料

沒有留言:

張貼留言

熱門文章

2018年4月4日 星期三

圓周運動

程式 7-1：圓周運動

參數設定

畫面設定

物體運動

程式 7-2.圓周運動, 加上轉軸及繩子

程式7-3.鉛直面圓周運動

理論分析

參數設定

畫面設定

自訂函式

物體運動

數據處理部分

結語

參考資料

沒有留言:

張貼留言

2018年4月4日星期三