作者:王一哲
日期:2018/4/5
將小球用一條理想的繩子懸掛在天花板底下,若繩子與鉛垂線的夾角(擺角)為
若定義轉動慣量 (moment of inertia)
由於力矩
若
可以解出
上式中
以下共有2個程式:
"""
VPython教學: 9-1.單擺
Ver. 1: 2018/2/25
Ver. 2: 2019/9/8
Ver. 3: 2023/4/27 懸掛點改成原點,所有物件向上平移
作者: 王一哲
"""
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
alpha = 0 # 角加速度, 初始值為 0
omega = 0 # 角速度, 初始值為 0
i = 0 # 小球經過週期次數
t = 0 # 時間
dt = 0.001 # 時間間隔
print("週期理論值 T = {:.12f} s".format(T))
"""
2. 畫面設定
"""
# 產生動畫視窗、天花板、小球、繩子
scene = canvas(title="Pendulum", width=600, height=600, x=0, y=0, background=vec(0, 0.6, 0.6), center=vec(0, -L/2, 0), range=L)
roof = box(pos=vec(0, 0.05, 0), size=vec(L, 0.1, 0.5*L), color=color.blue)
ball = sphere(pos=vec(L*sin(theta0), -L*cos(theta0), 0), radius=size, color=color.red,
make_trail=True, retain=100, v=vec(0, 0, 0))
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), 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
alpha = cross(r, vec(0, -m*g, 0)).z/(m*L**2)
#alpha = -m*g*ball.pos.x/(m*L**2)
omega += alpha*dt
theta += omega*dt
# 更新小球的位置、速度, 繩子的軸方向及長度
ball.pos = vec(L*sin(theta), -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("第{:d}次擺動經過的時間 t = {:.12f} s".format(i, t))
omega_p = omega_c
# 更新時間
t += dt
在此定義的變數有 m、size、L、theta0、theta、g、T、alpha、omega、i、t、dt,用途都已經寫在該行的註解中。
由於我們將小球的懸掛點改到原點,將天花板、小球、繩子都向上平移,因此計算 r 等於小球的位置。
r = ball.pos
計算角加速度的方法有兩種,如果要用數學上的外積 (cross product)先計算力矩,由於
alpha = cross(r, vec(0, -m*g, 0)).z/(m*L**2)
第二種寫法,用力乘以力臂計算重力產生的力矩,再計算角加速度
alpha = -m*g*ball.pos.x / (m*L**2)
再計算角速度、擺角。
更新小球的位置、速度, 繩子的軸方向及長度。
更新代表速度的箭頭位置、方向、長度。
畫出
當小球經過一個週期時,角度速會由逆時鐘方向(正值)變為順時鐘方向(負值),藉此判斷小球是否經過一個週期。
最後記得要更新時間。
當
"""
VPython教學: 9-2.單擺, 有空氣阻力
Ver. 1: 2018/2/25
Ver. 2: 2019/9/8
Ver. 3: 2023/4/27 懸掛點改成原點,所有物件向上平移
作者: 王一哲
"""
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,
alpha = 0 # 角加速度, 初始值為 0
omega = 0 # 角速度, 初始值為 0
i = 0 # 小球經過週期次數
t = 0 # 時間
dt = 0.001 # 時間間隔
print("週期理論值 T = {:.12f} s".format(T))
"""
2. 畫面設定
"""
# 產生動畫視窗、天花板、小球、繩子
scene = canvas(title="Pendulum with Air Drag", width=600, height=600, x=0, y=0, background=vec(0, 0.6, 0.6), center=vec(0, -L/2, 0), range=L)
roof = box(pos=vec(0, 0.05, 0), size=vec(L, 0.1, 0.5*L), color=color.blue)
ball = sphere(pos=vec(L*sin(theta0), -L*cos(theta0), 0), radius=size, color=color.red,
make_trail=True, retain=100, v=vec(0, 0, 0))
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), 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
F = vec(0, -m*g, 0) - b*ball.v
alpha = cross(r, F).z/(m*L*L)
omega += alpha*dt
theta += omega*dt
# 更新小球的位置、速度, 繩子的軸方向及長度
ball.pos = vec(L*sin(theta), -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("第{:d}次擺動經過的時間 t = {:.12f} s".format(i, t))
omega_p = omega_c
# 更新時間
t += dt
程式 9-2 與 9-1 非常類似,不同之處在於
第 18 行,新增變數 b = 0.1。
第 58 行,計算小球所受合力時要加上空氣阻力
F = vec(0, -m*g, 0) - b * ball.v
取
"""
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
在前一篇文章〈簡諧運動〉畫出了木塊、彈簧系統的運動情形,在程式 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 行,改成畫位置、速度、加速度與時間關係圖。
取 𝛳0 = 30° ,k = 10 N/m,L0 = 4 m,小球由靜止釋放後的位置、速度、加速度與時間關係圖如下,我實在看不出有任何的規律性。有時候小球會超過出發時的高度,這是有可能發生的,因為小球出發時除了有重力位能之外還有小球和彈簧之間的彈性位能,如果小球正好在較比較靠近懸掛點處向上移動,就有可能把大部分的彈性位能換成重力位能。
即使系統比較複雜,但只要知道如何描述物體的受力情形,就可以用 VPython 把運動過程畫出來。
VPython