一個地球人在臺北: 單擺

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

將小球用一條理想的繩子懸掛在天花板底下，若繩子與鉛垂線的夾角（擺角）為 $\theta_0$，小球由靜止釋放。當擺角為 $\theta$ 時，小球受到重力及繩子張力的作用，相對於懸掛點產生的力矩為

$$ \vec \tau = \vec r \times \vec F ~\Rightarrow~ \tau = -rF \sin \theta $$

$$ \tau = I \alpha ~\Rightarrow~ -rF \sin \theta = mr^2 \frac{d^2 \theta}{dt^2} $$

由於上式中 $\tau$ 與 $\theta$ 的方向相反，需要加上負號。若 $r = L$，$F = mg$，當 $\theta < 5^{\circ}$ 時， $\sin \theta \approx \theta$，可以解出 $$ \theta (t) = \theta_0 \cos \left(\sqrt{\frac{g}{L}} t \right) $$ 週期為 $$ T = 2 \pi \sqrt{\frac{L}{g}} $$

以下共有2個程式：

理想的單擺，改變起始的擺角計算運動過程及週期。　（GlowScript 網站動畫連結）
考慮空氣阻力的單擺。　（GlowScript 網站動畫連結）

單擺畫面截圖

程式 9-1：單擺

取得程式碼

"""
 VPython教學: 9-1.單擺
 Ver. 1: 2018/2/25
 Ver. 2: 2019/9/8
 作者: 王一哲
"""
from vpython import *

"""
 1. 參數設定, 設定變數及初始值
"""
m = 1               # 小球質量
size = 0.2          # 小球半徑
L = 5               # 擺長
theta0 = radians(30)# 起始擺角, 用 radians 將單位換成 rad
theta = theta0      # 擺角
g = 9.8             # 重力加速度
T = 2*pi*sqrt(L/g)  # 單擺週期理論值, L=5, g=9.8, T=4.48798950512828
alpha = 0           # 角加速度, 初始值為 0
omega = 0           # 角速度, 初始值為 0
i = 0               # 小球經過週期次數
t = 0               # 時間
dt = 0.001          # 時間間隔

"""
 2. 畫面設定
"""
# 產生動畫視窗、天花板、小球、繩子
scene = canvas(title="Pendulum", width=600, height=600, x=0, y=0, background=vec(0, 0.6, 0.6))
roof = box(pos=vec(0, L/2 + 0.05, 0), size=vec(L, 0.1, 0.5*L), color=color.blue)
ball = sphere(pos=vec(L*sin(theta0), L/2 - L*cos(theta0), 0), radius=size, color=color.red,
              make_trail=True, retain=100, v=vec(0, 0, 0))
rope = cylinder(pos=vec(0, L/2, 0), axis=ball.pos - vec(0, L/2, 0), radius=0.1*size, color=color.yellow)
# 產生表示速度的箭頭
arrow_v = arrow(pos=ball.pos, axis=vec(0, 0, 0), shaftwidth=0.3*size, color=color.green)
arrow_vx = arrow(pos=ball.pos, axis=vec(0, 0, 0), shaftwidth=0.3*size, color=color.magenta)
arrow_vy = arrow(pos=ball.pos, axis=vec(0, 0, 0), shaftwidth=0.3*size, color=color.orange)
# 繪圖部分
gd = graph(title="plot", width=600, height=450, x=0, y=600, xtitle="<i>t</i> (s)", 
           ytitle="blue: theta (rad), green: omega (rad/s), red: alpha (rad/s<sup>2</sup>)")
theta_t = gcurve(graph=gd, color=color.blue)
omega_t = gcurve(graph=gd, color=color.green)
alpha_t = gcurve(graph=gd, color=color.red)

"""
 3. 物體運動部分, 重覆5個週期
"""
omega_p = omega

while(i < 5):
    rate(1000)
# 計算小球所受力矩、角加速度、角速度、擺角
    r = ball.pos - vec(0, L/2, 0)
    F = vec(0, -m*g, 0)
    torque = cross(r, F)
    alpha = torque.z/(m*L*L)
    #alpha = -m*g*ball.pos.x/(m*L*L)
    omega += alpha*dt
    theta += omega*dt
# 更新小球的位置、速度, 繩子的軸方向及長度
    ball.pos = vec(L*sin(theta), L/2 - L*cos(theta), 0)
    v = L*omega
    vx = v*cos(theta)
    vy = v*sin(theta)
    rope.axis = r
# 更新代表速度的箭頭位置、方向、長度
    arrow_v.pos = ball.pos
    arrow_vx.pos = ball.pos
    arrow_vy.pos = ball.pos
    arrow_v.axis = vec(vx, vy, 0)
    arrow_vx.axis = vec(vx, 0, 0)
    arrow_vy.axis = vec(0, vy, 0)
# 畫出 theta-t, omega-t, alpha-t 圖
    theta_t.plot(pos=(t, theta))
    omega_t.plot(pos=(t, omega))
    alpha_t.plot(pos=(t, alpha))
# 檢驗小球是否經過一個週期
    omega_c = omega
    if(omega_p > 0 and omega_c < 0):
        i += 1
        print(i, t)
    omega_p = omega_c
# 更新時間
    t += dt

參數設定

在此定義的變數有 m、size、L、theta0、theta、g、T、alpha、omega、i、t、dt，用途都已經寫在該行的註解中。

畫面設定

由於我希望物件可以盡量佔滿整個動畫視窗，因此將天花板的位置向上移動到 y = L/2 + 0.05 處，懸掛點不在原點，計算繩子的長度及方向時要減掉天花板的中心位置。
配合懸掛點的位置，小球位於 y = L/2 - L*cos(theta0) 處。
目前只畫出代表 v、vx、vy 的箭頭，可以仿照這個方法畫出代表 a、ax、ay 的箭頭。
繪圖部分只畫出 $\theta - t$、$\omega - t$、$\alpha - t$ 關係圖。

物體運動

先用
```
r = ball.pos - vec(0, L/2, 0)
```
計算並更新繩子的軸方向及長度，由於只有重力會產生力矩，因此只接計算角加速度
```
alpha = -m * g * ball.pos.x / (m * L * L)
```
再計算角速度、擺角。
更新小球的位置、速度, 繩子的軸方向及長度。
更新代表速度的箭頭位置、方向、長度。
畫出 $\theta - t$、$\omega - t$、$\alpha - t$ 關係圖。
當小球經過一個週期時，角度速會由逆時鐘方向（正值）變為順時鐘方向（負值），藉此判斷小球是否經過一個週期。
最後記得要更新時間。

模擬結果

當 $\theta < 5^{\circ}$ 時，週期理論值為 4.48798950512828 s；當 $\theta = 1^{\circ}$ 時週期為 4.4879999999998335 s；當 $\theta = 30^{\circ}$ 時週期為 4.5659999999998595；當 $\theta = 60^{\circ}$ 時週期為 4.815999999999943，當 $\theta$ 越大時週期與理論值的差距越大。從 $\theta - t$、$\omega - t$、$\alpha - t$ 關係圖可以看出，三者大約是 cos、-sin、-cos 的關係，與理論計算相符。

𝛳 = 1°，𝛳 - t 、𝜔 - t 、𝛼 - t 關係圖

𝛳 = 30°，𝛳 - t 、𝜔 - t 、𝛼 - t 關係圖

𝛳 = 60°，𝛳 - t 、𝜔 - t 、𝛼 - t 關係圖

程式 9-2.單擺, 有空氣阻力

取得程式碼

"""
 VPython教學: 9-2.單擺, 有空氣阻力
 Ver. 1: 2018/2/25
 Ver. 2: 2019/9/8
 作者: 王一哲
"""
from vpython import *

"""
 1. 參數設定, 設定變數及初始值
"""
m = 1                    # 小球質量
size = 0.2               # 小球半徑
L = 5                    # 擺長
theta0 = radians(30)     # 起始擺角, 用 radians 將單位換成 rad
theta = theta0           # 擺角
g = 9.8                  # 重力加速度
b = 0.1                  # 空氣阻力 f = -bv
T = 2*pi*sqrt(L/g)       # 單擺週期理論值, L = 5, g = 9.8, T = 4.48798950512828
alpha = 0                # 角加速度, 初始值為 0
omega = 0                # 角速度, 初始值為 0
i = 0                    # 小球經過週期次數
t = 0                    # 時間
dt = 0.001               # 時間間隔

"""
 2. 畫面設定
"""
# 產生動畫視窗、天花板、小球、繩子
scene = canvas(title="Pendulum", width=600, height=600, x=0, y=0, background=vec(0, 0.6, 0.6))
roof = box(pos=vec(0, L/2 + 0.05, 0), size=vec(L, 0.1, 0.5*L), color=color.blue)
ball = sphere(pos=vec(L*sin(theta0), L/2 - L*cos(theta0), 0), radius=size, color=color.red,
              make_trail=True, retain=100, v=vec(0, 0, 0))
rope = cylinder(pos=vec(0, L/2, 0), axis=ball.pos - vec(0, L/2, 0), radius=0.1*size, color=color.yellow)
# 產生表示速度的箭頭
arrow_v = arrow(pos=ball.pos, axis=vec(0, 0, 0), shaftwidth=0.3*size, color=color.green)
arrow_vx = arrow(pos=ball.pos, axis=vec(0, 0, 0), shaftwidth=0.3*size, color=color.magenta)
arrow_vy = arrow(pos=ball.pos, axis=vec(0, 0, 0), shaftwidth=0.3*size, color=color.orange)
# 繪圖部分
gd = graph(title="plot", width=600, height=450, x=0, y=600, xtitle="<i>t</i> (s)", 
           ytitle="blue: theta (rad), green: omega (rad/s), red: alpha (rad/s<sup>2</sup>)")
theta_t = gcurve(graph=gd, color=color.blue)
omega_t = gcurve(graph=gd, color=color.green)
alpha_t = gcurve(graph=gd, color=color.red)

"""
 3. 物體運動部分, 重覆5個週期
"""
omega_p = omega

while(i < 5):
    rate(1000)
# 計算小球所受力矩、角加速度、角速度、角位移
    r = ball.pos - vec(0, L/2, 0)
    F = vec(0, -m*g, 0) - b*ball.v
    torque = cross(r, F)
    alpha = torque.z/(m*L*L)
    omega += alpha*dt
    theta += omega*dt
# 更新小球的位置、速度, 繩子的軸方向及長度
    ball.pos = vec(L*sin(theta), L/2 - L*cos(theta), 0)
    v = L*omega
    vx = v*cos(theta)
    vy = v*sin(theta)
    ball.v = vec(vx, vy, 0)
    rope.axis = r
# 更新代表速度的箭頭位置、方向、長度
    arrow_v.pos = ball.pos
    arrow_vx.pos = ball.pos
    arrow_vy.pos = ball.pos
    arrow_v.axis = vec(vx, vy, 0)
    arrow_vx.axis = vec(vx, 0, 0)
    arrow_vy.axis = vec(0, vy, 0)
# 畫出 theta-t, omega-t, alpha-t 圖
    theta_t.plot(pos = (t, theta))
    omega_t.plot(pos = (t, omega))
    alpha_t.plot(pos = (t, alpha))
# 檢驗小球是否經過一個週期
    omega_c = omega
    if(omega_p > 0 and omega_c < 0):
        i += 1
        print(i, t)
    omega_p = omega_c
# 更新時間
    t += dt

程式設計部分

程式 9-2 與 9-1 非常類似，不同之處在於

第 18 行，新增變數 b = 0.1。
第 55 ~ 57 行，先用
```
F = vec(0, -m * g, 0) - b * ball.v
```
計算小球所受合力，再用
```
torque = cross(r, F)
```
計算相對於懸掛點的力矩，最後用
```
alpha = torque.z / (m * L * L)
```
計算角加速度。

模擬結果

取 $\theta_0 = 30^{\circ}$，由以下4張圖可以看出當 b 越大時，𝛳、𝜔、𝛼 隨時間衰減的幅度較大，符合理論預測。

b = 0.1，𝛳 - t 、𝜔 - t 、𝛼 - t 關係圖

b = 0.2，𝛳 - t 、𝜔 - t 、𝛼 - t 關係圖

b = 0.3，𝛳 - t 、𝜔 - t 、𝛼 - t 關係圖

b = 0.4，𝛳 - t 、𝜔 - t 、𝛼 - t 關係圖

程式 9-3.單擺 + 彈簧

取得程式碼
GlowScript 網站動畫連結

"""
 VPython教學: 9-3.單擺 + 彈簧
 Ver. 1: 2018/4/5
 Ver. 2: 2019/9/8
 作者: 王一哲
"""
from vpython import *

"""
 1. 參數設定, 設定變數及初始值
"""
m = 1                    # 小球質量
size = 0.2               # 小球半徑
L = 5                    # 擺長
theta0 = radians(30)     # 起始擺角, 用 radians 將單位換成 rad
theta = theta0           # 擺角
g = 9.8                  # 重力加速度
k = 10                   # 彈性常數 10 N/m
L0 = 4                   # 彈簧原長
T = 2*pi*sqrt(L/g)       # 單擺週期理論值
ratio = 0.2              # 箭頭長度與實際值的比例
t = 0                    # 時間
dt = 0.001               # 時間間隔

"""
 2. 畫面設定
"""
# 產生動畫視窗、天花板、小球、彈簧
scene = canvas(title="Pendulum with Spring", width=600, height=600, x=0, y=0, background=vec(0, 0.6, 0.6))
roof = box(pos = vec(0, L/2 + 0.05, 0), size=vec(L, 0.1, 0.5*L), color=color.blue)
ball = sphere(pos=vec(L*sin(theta0), L/2 - L*cos(theta0), 0), radius=size,
              color=color.red, make_trail=True, retain=100, v=vec(0, 0, 0))
spring = helix(pos=vec(0, L/2, 0), axis=ball.pos - vec(0, L/2, 0), radius=0.6*size,
               thickness=0.3*size, color=color.yellow)
# 產生表示速度的箭頭
arrow_v = arrow(pos=ball.pos, axis=vec(0, 0, 0), shaftwidth=0.3*size, color=color.green)
arrow_vx = arrow(pos=ball.pos, axis=vec(0, 0, 0), shaftwidth=0.3*size, color=color.magenta)
arrow_vy = arrow(pos=ball.pos, axis=vec(0, 0, 0), shaftwidth=0.3*size, color=color.orange)
# 繪圖部分
gd1 = graph(title="position", width=600, height=450, x=0, y=600, xtitle="<i>t</i>(s)", 
            ytitle="blue: <i>r</i>(m), red: <i>x</i>(m), green: <i>y</i>(m)")
r_t = gcurve(graph=gd1, color=color.blue)
x_t = gcurve(graph=gd1, color=color.red)
y_t = gcurve(graph=gd1, color=color.green)
gd2 = graph(title="velocity", width=600, height=450, x=0, y=1050, xtitle="<i>t</i>(s)", 
            ytitle="blue: <i>v</i>(m/s), red: <i>v<sub>x</sub></i>(m/s), green: <i>v<sub>y</sub></i>(m/s)")
v_t = gcurve(graph=gd2, color=color.blue)
vx_t = gcurve(graph=gd2, color=color.red)
vy_t = gcurve(graph=gd2, color=color.green)
gd3 = graph(title="acceleration", width=600, height=450, x=0, y=1500, xtitle="<i>t</i>(s)",
            ytitle = "blue: <i>a</i>(m/s<sup>2</sup>), red: <i>a<sub>x</sub></i>(m/s<sup>2</sup>), green: <i>a<sub>y</sub></i>(m/s<sup>2</sup>)")
a_t = gcurve(graph=gd3, color=color.blue)
ax_t = gcurve(graph=gd3, color=color.red)
ay_t = gcurve(graph=gd3, color=color.green)

"""
 3. 物體運動部分
"""
r = ball.pos - vec(0, L/2, 0)
while(t < 5*T):
    rate(1000)
# 計算小球所受合力、加速度、速度、位置, 更新彈簧長度、方向
    F = vec(0, -m*g, 0) - k*(r.mag - L0)*r.norm()
    ball.a = F/m
    ball.v += ball.a*dt
    ball.pos += ball.v*dt
    r = ball.pos - vec(0, L/2, 0)
    spring.axis = r
# 更新代表速度的箭頭位置、方向、長度
    arrow_v.pos = ball.pos
    arrow_vx.pos = ball.pos
    arrow_vy.pos = ball.pos
    arrow_v.axis = ball.v * ratio
    arrow_vx.axis = vec(ball.v.x, 0, 0) * ratio
    arrow_vy.axis = vec(0, ball.v.y, 0) * ratio
# 畫出位置、速度、加速度與時間關係圖
    r_t.plot(pos = (t, ball.pos.mag))
    x_t.plot(pos = (t, ball.pos.x))
    y_t.plot(pos = (t, ball.pos.y))
    v_t.plot(pos = (t, ball.v.mag))
    vx_t.plot(pos = (t, ball.v.x))
    vy_t.plot(pos = (t, ball.v.y))
    a_t.plot(pos = (t, ball.a.mag))
    ax_t.plot(pos = (t, ball.a.x))
    ay_t.plot(pos = (t, ball.a.y))
# 更新時間
    t += dt

將繩子改為彈簧，𝛳₀ = 30°，不考慮空氣阻力的運動畫面截圖

程式設計部分

在前一篇文章〈簡諧運動〉畫出了木塊、彈簧系統的運動情形，在程式 9-1 又畫了單擺的運動情形，於是我很好奇將單擺系統中的繩子換成彈簧會發生什麼事，就拿程式 9-1 來修改一下，以下是修改的地方：

第 18、19 行，新增彈性常數 k = 10、彈簧原長 L0 = 4。
為了使代表速度的箭頭不會太長，在第 21 行新增箭頭長度與實際值的比例 ratio = 0.2。
第 33 行，將原來產生的繩子改用 helix 產生彈簧。
第 40 ~ 54 行，繪圖部分共畫 3 張圖，將位置、速度、加速度與時間關係圖分開畫。
第 59 行，在 while 迴圈外先用
```
r = ball.pos - vec(0, L/2, 0)
```
計算軸的長度及方向。

第 63 ~ 66 行，先用

F = vec(0, -m * g, 0) - k * (r.mag - L0) * r.norm()

計算小球所受合力，再用

ball.a = F/m
ball.v += ball.a * dt
ball.pos += ball.v * dt

計算小球的加速度、更新速度及位置。

第 74 ~ 75 行，為了避免箭頭太長影響到畫面的比例，在更新箭頭的長度及方向時乘以 ratio，使長度縮短。
第 77 ~ 85 行，改成畫位置、速度、加速度與時間關係圖。

模擬結果

取 𝛳₀ = 30° ，k = 10 N/m，L0 = 4 m，小球由靜止釋放後的位置、速度、加速度與時間關係圖如下，我實在看不出有任何的規律性。有時候小球會超過出發時的高度，這是有可能發生的，因為小球出發時除了有重力位能之外還有小球和彈簧之間的彈性位能，如果小球正好在較比較靠近懸掛點處向上移動，就有可能把大部分的彈性位能換成重力位能。

將繩子改為彈簧，𝛳₀ = 30°，小球的位置 - 時間關係圖

將繩子改為彈簧，𝛳₀ = 30°，小球的速度 - 時間關係圖

將繩子改為彈簧，𝛳₀ = 30°，小球的加速度 - 時間關係圖

結語

即使系統比較複雜，但只要知道如何描述物體的受力情形，就可以用 VPython 把運動過程畫出來。

參考資料

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

一個地球人在臺北

熱門文章

2018年4月5日星期四

單擺

程式 9-1：單擺

參數設定

畫面設定

物體運動

模擬結果

程式 9-2.單擺, 有空氣阻力

程式設計部分

模擬結果

程式 9-3.單擺 + 彈簧

程式設計部分

模擬結果

結語

參考資料

沒有留言:

張貼留言

熱門文章

2018年4月5日 星期四

單擺

程式 9-1：單擺

參數設定

畫面設定

物體運動

模擬結果

程式 9-2.單擺, 有空氣阻力

程式設計部分

模擬結果

程式 9-3.單擺 + 彈簧

程式設計部分

模擬結果

結語

參考資料

沒有留言:

張貼留言

2018年4月5日星期四