一個地球人在臺北: 簡諧運動

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

在水平光滑桌面上有一個木塊質量為 m，用一條彈性係數為 k 的彈簧連接到左側的牆壁上，若將木塊向右拉一段距離 R 再由靜止釋放，木塊所受合力與加速度的關係為

$$ F = -kx = ma ~\Rightarrow -kx = m \frac{d^2 x}{dt^2} $$

此時木塊的運動方式稱為簡諧運動(simple harmonic motion, S.H.M.)，由上式可以解出

$$ x(t) = R \cos (\omega t + \phi) $$

$$ v(t) = -\omega R \sin (\omega t + \phi) $$

$$ a(t) = -\omega^2 R \cos (\omega t + \phi) $$

上式中 $\omega$ 為角頻率

$$ \omega = \sqrt{\frac{k}{m}} $$

週期為

$$ T = 2\pi \sqrt{\frac{m}{k}} $$

理論上我們只要在 VPython 中設定好木塊所受的彈簧回復力與木塊離開平衡點位置的關係，應該就能畫出簡諧運動的過程與週期。以下共有2個程式：

理想的簡諧運動。　（GlowScript 網站動畫連結）
考慮阻尼(damping)的簡諧運動。　（GlowScript 網站動畫連結）

簡諧運動畫面截圖

程式 8-1：簡諧運動

取得程式碼

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

"""
 1. 參數設定, 設定變數及初始值
"""
m = 4               # 木塊質量 4 kg
size = 1            # 木塊邊長 1 m
R = 5               # 振幅 5 m
k = 1               # 彈性常數 1 N/m
L0 = R + size       # 彈簧原長
i = 0               # 木塊回到初位置的次數
t = 0               # 時間
dt = 0.001          # 時間間隔

"""
 2. 畫面設定
"""
# 產生動畫視窗、地板、木塊、彈簧
scene = canvas(title="Simple Harmonic Motion", width=800, height=400, x=0, y=0, background=vec(0, 0.6, 0.6))
floor = box(pos=vec(0, -(size+0.1)/2, 0), size=vec(2*L0, 0.1, R), texture=textures.metal)
wall = box(pos=vec(-L0, 0, 0), size=vec(0.1, size, R), texture=textures.metal)
block = box(pos=vec(R+size/2, 0, 0), size=vec(size, size, size), texture=textures.wood, v=vec(0, 0, 0))
spring = helix(pos=vec(-L0, 0, 0), radius=0.3*size, thickness=0.05*size, color=color.yellow)
spring.axis = block.pos - spring.pos - vec(size/2, 0, 0)
# 產生表示速度、加速度的箭頭
arrow_v = arrow(pos=block.pos + vec(0, size, 0), axis=vec(0, 0, 0), shaftwidth=0.3*size, color=color.green)
arrow_a = arrow(pos=block.pos + vec(0, 2*size, 0), axis=vec(0, 0, 0), shaftwidth=0.3*size, color=color.magenta)
# 繪圖部分
gd = graph(title="plot", width=600, height=450, x=0, y=400, xtitle="<i>t</i>(s)", 
           ytitle="blue: <i>x</i>(m), green: <i>v</i>(m/s), magenta: <i>a</i>(m/s<sup>2</sup>)")
xt = gcurve(graph=gd, color=color.blue)
vt = gcurve(graph=gd, color=color.green)
at = gcurve(graph=gd, color=color.magenta)

"""
 3. 物體運動部分, 重覆3個週期
"""
while(i < 3):
    rate(1000)
# 計算彈簧長度、伸長量、回復力
    spring.axis = block.pos - spring.pos - vec(0.5*size, 0, 0)
    F = -k * (spring.axis - vec(L0, 0, 0))
# 計算木塊加速度, 更新速度、位置
    block.a = F/m
    block.v += block.a*dt
    block.pos += block.v*dt
# 更新代表速度、加速度的箭頭位置、方向、長度
    arrow_v.pos = block.pos + vec(0, size, 0)
    arrow_a.pos = block.pos + vec(0, 2*size, 0)
    arrow_v.axis = block.v
    arrow_a.axis = block.a
# 畫出 x-t, v-t, a-t 圖
    xt.plot(pos=(t, block.pos.x - size/2))
    vt.plot(pos=(t, block.v.x))
    at.plot(pos=(t, block.a.x))
# 判斷木塊是否回到出發點, 計算回到出發點的次數
    if(block.pos.x > R + size/2 and block.v.x > 0):
        print(i, t)
        i += 1
# 更新時間
    t += dt

參數設定

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

畫面設定

由於我希望木塊的中心位於 x 軸上，因此將地板的位置向下移動到 y = -(size+0.1)/2 處。
為了讓木塊在 $-R \leq x \leq R$ 之間來回運動，因此將牆壁向左移動到 x = -L0 處。
由於木塊的位置是以中心點為準，而且彈簧連接在木塊左，因此需要將木塊放在彈簧右端點的右側 size/2 處。
為了畫出彈簧需要用到一個新的物件 helix（螺旋線），可以調整的選項主要有
- 位置(pos)，以螺旋線的端點為基準。
- 半徑(radius)
- 軸(axis)，以位置為起點，格式為向量。
- 厚度(thickness)，線條本身的半徑，預設值為 radius/20。
- 顏色(color)

物體運動

先用 spring.axis = block.pos - spring.pos - vector(size/2, 0, 0) 計算並更新彈簧的軸方向及長度，再用 F = -k * (spring.axis - vector(L0, 0, 0)) 計算彈簧的伸長量、回復力。
判斷木塊是否回到出發點的部分有很多種寫法，我是用 block.pos.x > R + size/2 且木塊速度方向向右來判斷，若計算小球回到出發點時將次數 i 加 1 並印出時間。

模擬結果

由下圖可以看出 x、v、a 與時間的關係符合理論預測，的確是 cos、-sin、-cos 的型式。週期模擬值為 12.565，理論值為 12.566，相當接近。

m = 4、k = 1、R = 5 的簡諧運動 x - t、v - t、a - t 關係圖

程式 8-2：簡諧運動, 有阻尼

取得程式碼

"""
 VPython教學: 8-2.簡諧運動, 有阻尼
 Ver. 1: 2018/2/25
 Ver. 2: 2019/9/8
 作者: 王一哲
"""
from vpython import *

"""
 1. 參數設定, 設定變數及初始值
"""
m = 4               # 木塊質量 4 kg
size = 1            # 木塊邊長 1 m
R = 5               # 振幅 5 m
k = 1               # 彈性常數 1 N/m
L0 = R + size       # 彈簧原長
b = 0.3               # 阻尼 f = -bv, overdamped: b^2 > 4mk, critical damping: b^2 = 4mk, underdamped: b^2 < 4mk
T = 2*pi*sqrt(m/k)  # 週期理論值
i = 0               # 木塊運動經過的週期次數
t = 0               # 時間
dt = 0.001          # 時間間隔

"""
 2. 畫面設定
"""
# 產生動畫視窗、地板、木塊、彈簧
scene = canvas(title="Simple Harmonic Motion", width=800, height=400, x=0, y=0, background=vec(0, 0.6, 0.6))
floor = box(pos=vec(0, -(size+0.1)/2, 0), size=vec(2*L0, 0.1, R), texture=textures.metal)
wall = box(pos=vec(-L0, 0, 0), size=vec(0.1, size, R), texture=textures.metal)
block = box(pos=vec(R+size/2, 0, 0), size=vec(size, size, size), texture=textures.wood, v=vec(0, 0, 0))
spring = helix(pos=vec(-L0, 0, 0), radius=0.3*size, thickness=0.05*size, color=color.yellow)
spring.axis = block.pos - spring.pos - vec(size/2, 0, 0)
# 產生表示速度、加速度的箭頭
arrow_v = arrow(pos=block.pos + vec(0, size, 0), axis=vec(0, 0, 0), shaftwidth=0.3*size, color=color.green)
arrow_a = arrow(pos=block.pos + vec(0, 2*size, 0), axis=vec(0, 0, 0), shaftwidth=0.3*size, color=color.magenta)
# 繪圖部分
gd = graph(title="plot", width=600, height=450, x=0, y=400, xtitle="<i>t</i>(s)",
           ytitle="blue: <i>x</i>(m), green: <i>v</i>(m/s), magenta: <i>a</i>(m/s<sup>2</sup>)")
xt = gcurve(graph=gd, color=color.blue)
vt = gcurve(graph=gd, color=color.green)
at = gcurve(graph=gd, color=color.magenta)

"""
 3. 物體運動部分, 重覆5個週期
"""
vp = block.v.x
while(i < 5 and t < 5*T):
    rate(1000)
# 計算彈簧長度、伸長量、回復力
    spring.axis = block.pos - spring.pos - vec(0.5*size, 0, 0)
    F = -k * (spring.axis - vec(L0, 0, 0)) - b * block.v
# 計算木塊加速度, 更新速度、位置
    block.a = F/m
    block.v += block.a*dt
    block.pos += block.v*dt
# 更新代表速度、加速度的箭頭位置、方向、長度
    arrow_v.pos = block.pos + vec(0, size, 0)
    arrow_a.pos = block.pos + vec(0, 2*size, 0)
    arrow_v.axis = block.v
    arrow_a.axis = block.a
# 畫出 x-t, v-t, a-t 圖
    xt.plot(pos = (t, block.pos.x - size/2))
    vt.plot(pos = (t, block.v.x))
    at.plot(pos = (t, block.a.x))
# 判斷木塊是否經過一個週期
    vc = block.v.x
    if(vp > 0 and vc < 0):
        i += 1
        print(i, t)
    vp = vc
# 更新時間
    t += dt

理論計算

假設木塊所受阻尼 $f = -bv$，由牛頓第二定律可得

$$ ma + bv + kx = 0 $$

$$ m \frac{d^2 x}{dt^2} + b \frac{dx}{dt} + kx = 0 $$

共有 3 種情形：

$b^2 > 4mk ~~~~~$過阻尼 (overdamped)
$b^2 = 4mk ~~~~~$臨界阻尼 (critical damping)
$b^2 < 4mk ~~~~~$次阻尼 (underdamped)

程式設計部分

程式 8-2 與 8-1 很像，只有 2 個不同之處。

在計算木塊所受合力時需要改成 F = -k * (spring.axis - vec(L0, 0, 0)) - b * block.v
由於木塊不會回到出發點，在判斷木塊是否經過一個週期時，改成用速度來判斷，如果木塊的速度原來向右、後來向左，代表木塊經過一個週期。

vc = block.v.x
if(vp > 0 and vc < 0):
    i += 1
    print(i, t)
vp = vc

b = 6，過阻尼

b = 4，臨界阻尼

b = 1，次阻尼

b = 0.3，次阻尼

參考資料

canvas: http://www.glowscript.org/docs/VPythonDocs/canvas.html
box: http://www.glowscript.org/docs/VPythonDocs/box.html
helix: http://www.glowscript.org/docs/VPythonDocs/helix.html
arrow: http://www.glowscript.org/docs/VPythonDocs/arrow.html
graph: http://www.glowscript.org/docs/VPythonDocs/graph.html
Damped Harmonic Oscillator: http://hyperphysics.phy-astr.gsu.edu/hbase/oscda.html

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

一個地球人在臺北

2018年4月5日星期四

簡諧運動

程式 8-1：簡諧運動

參數設定

畫面設定

物體運動

模擬結果

程式 8-2：簡諧運動, 有阻尼

理論計算

程式設計部分

參考資料

沒有留言:

張貼留言

2018年4月5日 星期四

簡諧運動

程式 8-1：簡諧運動

參數設定

畫面設定

物體運動

模擬結果

程式 8-2：簡諧運動, 有阻尼

理論計算

程式設計部分

參考資料

沒有留言:

張貼留言

2018年4月5日星期四