簡諧運動

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

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

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

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

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

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

$$a(t)=-{\omega }^{2}R\cos (\omega t+\varphi )$$

上式中

$\omega$ 為角頻率

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

週期為

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

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

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

簡諧運動畫面截圖

程式 8-1：簡諧運動 （取得程式碼）

"""
 VPython教學: 8-1.簡諧運動
 Ver. 1: 2018/2/25
 Ver. 2: 2019/9/8
 Ver. 3: 2022/4/21 加上週期理論值, 修改 print 的格式
 作者: 王一哲
"""
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          # 時間間隔
print("週期理論值為 {:f} s".format(2*pi*sqrt(m/k)))

"""
 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:
        i += 1
        print("第 {:d} 次回到右端點，經過時間為 {:f} s".format(i, t))
# 更新時間
    t += dt

參數設定

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

畫面設定

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

物體運動

1. 先用 spring.axis = block.pos - spring.pos - vector(size/2, 0, 0) 計算並更新彈簧的軸方向及長度，再用 F = -k * (spring.axis - vector(L0, 0, 0)) 計算彈簧的伸長量、回復力。
2. 判斷木塊是否回到出發點的部分有很多種寫法，我是用 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
 Ver. 3: 2022/4/21 加上週期理論值, 修改 print 的格式
 作者: 王一哲
"""
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.1*sqrt(4*m*k) # 阻尼 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          # 時間間隔
print("週期理論值為 {:f} s".format(T))

"""
 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("完成第 {:d} 次振盪，經過時間為 {:f} s".format(i, t))
    vp = vc
# 更新時間
    t += dt

理論計算 [6]

假設木塊所受阻尼

$f=-bv$，由牛頓第二定律可得

$$ma+bv+kx=0$$

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

共有 3 種情形：

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

程式設計部分

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

1. 在計算木塊所受合力時需要改成 F = -k * (spring.axis - vec(L0, 0, 0)) - b * block.v

2. 由於木塊不會回到出發點，在判斷木塊是否經過一個週期時，改成用速度來判斷，如果木塊的速度原來向右、後來向左，代表木塊經過一個週期。

   ​​​​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，次阻尼

參考資料

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

---

tags:VPython