一個地球人在臺北: 一維彈性碰撞

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

一維彈性碰撞是基礎物理2B下最後一章的內容。假設在水平光滑桌面上有兩個小球，質量分別為 $m_1$ 及 $m_2$，速度分別為 $v_1$ 及 $v_2$，由於兩個小球的碰撞過程不受外力，系統動量守恆；若碰撞過程沒有能量損失，則碰撞前、後兩小球的總動能相等；由以上兩個條件可以導出小球碰撞前後的速度關係式為

$$ v_1' = \frac{m_1 - m_2}{m_1 + m_2} v_1 + \frac{2m_2}{m_1 + m_2} v_2 $$

$$ v_2' = \frac{2m_1}{m_1 + m_2} v_1 + \frac{m_2 - m_1}{m_1 + m_2} v_2 $$

一維彈性碰撞的過程

以下共有兩個程式，第一個是直接代撞後速度公式，第二個則是在木塊間加上理想彈簧作為緩衝，畫出完整的碰撞過程。

程式15-1畫面截圖

程式15-2畫面截圖

程式 15-1.一維彈性碰撞公式

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

"""
 VPython教學: 15-1.一維彈性碰撞公式
 Ver. 1: 2018/2/28
 Ver. 2: 2019/9/13
 作者: 王一哲
"""

from vpython import *

"""
 1. 參數設定, 設定變數及初始值
    (1) m1=0.5, m2=0.1, v1=1.0, v2=0.0
    (2) m1=0.2, m2=0.4, v1=1.0, v2=-2.0
    (3) m1=0.4, m2=0.2, v1=2.5, v2=-0.5
"""
d1, m1, v1, c1 = 0.2, 0.5, 1.0, color.red        # 木塊1的寬度=0.1 m, 質量=0.5 kg, 初速, 紅色
d2, m2, v2, c2 = 0.2, 0.1, 0.0, color.green      # 木塊2的寬度=0.1 m, 質量=0.1 kg, 初速, 綠色
xmax, xmin = 2.0, -2.0                           # x 軸範圍
t, dt = 0, 0.0005                                # 時間, 畫面更新的時間間隔，單位為s, 原為0.001但能量不夠準確, 故改為0.0005

"""
 2. 畫面設定
"""
# 產生動畫視窗
scene = canvas(title="1 Dimension Collision", width=800, height=300, center=vec(0, 0.4, 0),
               background=vec(0, 0.6, 0.6))
# 產生地板
floor = box(pos=vec(0, -d1/2.0, 0), size=vec((xmax - xmin), 0.05, 0.8), color=color.blue)
# 產生左側木塊 b1,  右側木塊 b2 並設定初速度
b1 = box(pos=vec(xmin + 0.5*d1, 0, 0), size=vec(d1, d1, d1), color=c1, m=m1, v=vec(v1, 0, 0))
b2 = box(pos=vec(0, 0, 0), size=vec(d2, d2, d2), color=c2, m=m2, v=vec(v2, 0, 0))
# 繪圖部分
gd = graph(title="<i>v</i>-<i>t</i> plot", x=0, y=300, width=600, height=450,
           xtitle="<i>t</i> (s)", ytitle="red: <i>v</i><sub>1</sub>, green: <i>v</i><sub>2</sub> (m/s)")
vt1 = gcurve(graph=gd, color=c1)
vt2 = gcurve(graph=gd, color=c2)

# 自訂函式，一維彈性碰撞速度公式
def af_col_v(m1, m2, v1, v2):
    v1_prime = (m1-m2)/(m1+m2)*v1 + (2*m2)/(m1+m2)*v2
    v2_prime = (2*m1)/(m1+m2)*v1 + (m2-m1)/(m1+m2)*v2
    return (v1_prime, v2_prime)

"""
 3. 物體運動部分, 重複執行直到木塊抵達邊緣時
"""
# 印出木塊撞前速度
print("m1 =", b1.m, "m2 =", b2.m)
print(b1.v.x, b2.v.x)

while ((b2.pos.x <= xmax - d2/2) and (b1.pos.x >= xmin + d1/2)):
    rate(1000)
# 計算木塊間的距離, 若發生碰撞則計算撞後速度
    dx = b2.pos.x - b1.pos.x
    if(dx <= (d1 + d2)/2):
        b1.v.x, b2.v.x = af_col_v(b1.m, b2.m, b1.v.x, b2.v.x)
# 更新木塊的位置
    b1.pos += b1.v * dt
    b2.pos += b2.v * dt
# 畫 v-t 圖
    vt1.plot(pos=(t, b1.v.x))
    vt2.plot(pos=(t, b2.v.x))
# 更新時間
    t += dt

# 印出木塊撞後速度
print(b1.v.x, b2.v.x)

參數設定

在此設定變數為木塊的寬度、質量、初速、顏色，x軸的範圍、時間、時間間隔，其中時間間隔 dt 設定為 0.0005，這是因為設定為 0.001 時計算木塊能量的誤差較大，故選擇較小的數值。

畫面設定

產生動畫視窗、地板、木塊、繪圖視窗的程式碼在之前動畫當中已經出現很多次，這裡就不再贅述。

自訂函式

程式碼第 36 ~ 39 行設定了一個名為 af_col_v 的函式，輸入的參數為 (m1, m2, v1, v2)，將以上參數代入一維彈性碰撞速度公式，回傳撞後速度 v1_prime, v2_prime。我們在程式碼第 58 行呼叫這個函式

b1.v.x, b2.v.x = af_col_v(b1.m, b2.m, b1.v.x, b2.v.x)

由於回傳值有兩個，程式會將第1個回傳值指定給 b1.v.x，將第2個回傳值指定給 b2.v.x。

物體運動部分

當木塊抵達邊緣時停止動畫，因此 while 迴圈的條件設定為
```
(b2.pos.x <= xmax - d2/2) and (b1.pos.x >= xmin + d1/2)
```
若兩個木塊中心的距離小於 (d1 + d2) / 2，表示木塊之間發生碰撞，呼叫自訂函式 af_col_v 計算撞後速度。
更新木塊的位置，畫速度 - 時間關係圖。

模擬結果

我測試了 3 種不同的條件：

m1 = 0.5, m2 = 0.1, v1 = 1.0, v2 = 0.0
m1 = 0.2, m2 = 0.4, v1 = 1.0, v2 = -2.0
m1 = 0.4, m2 = 0.2, v1 = 2.5, v2 = -0.5

這是將數值直接代入一維彈性碰撞速度公式得到的結果，主要是作為下一個程式的對照組。

條件1的模擬結果：速度 - 時間關係圖

條件2的模擬結果：速度 - 時間關係圖

條件3的模擬結果：速度 - 時間關係圖

程式 15-2.木塊彈簧系統彈性碰撞

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

"""
 VPython教學: 15-2.木塊彈簧系統彈性碰撞
 Ver. 1: 2017/6/15
 Ver. 2: 2017/12/31 改為新版VPython語法
 Ver. 3: 2018/2/28 加上 v-t, a-t 圖
 Ver. 4: 2019/9/13
 作者: 王一哲
"""

from vpython import *

"""
 1. 參數設定, 設定變數及初始值
    (1) m1 = 0.5, m2 = 0.1, v1 = 1.0, v2 = 0.0, k = 2.0
    (2) m1 = 0.2, m2 = 0.4, v1 = 1.0, v2 = -2.0, k = 5.0
    (3) m1 = 0.4, m2 = 0.2, v1 = 2.5, v2 = -0.5
"""
d1, m1, v1, c1 = 0.2, 0.4, 2.5, color.red         # 木塊1的寬度 = 0.1 m, 質量 = 0.5 kg, 初速, 紅色
d2, m2, v2, c2 = 0.2, 0.2, -0.5, color.green       # 木塊2的寬度 = 0.1 m, 質量 = 0.1 kg, 初速, 綠色
L0, k = 0.5, 5.0                                  # 彈簧的原長 = 0.5 m, 彈性常數 = 2.0 N/m
xmax, xmin = 2.0, -2.0                            # x 軸範圍
t, dt = 0, 0.0005                                 # 時間, 畫面更新的時間間隔，單位為s, 原為0.001但能量不夠準確, 故改為0.00005

"""
 2. 畫面設定
"""
# 產生動畫視窗
scene = canvas(title="1 Dimension Collision", width=800, height=300, center=vec(0, 0.4, 0),
               background=vec(0, 0.6, 0.6))
# 產生地板
floor = box(pos=vec(0, -d1/2.0, 0), size=vec((xmax - xmin), 0.05, 0.8), color=color.blue)
# 產生左側木塊 b1,  右側木塊 b2 並設定初速度
b1 = box(pos=vec(-L0 - 1, 0, 0), size=vec(d1, d1, d1), color=c1, m=m1, v=vec(v1, 0, 0))
b2 = box(pos=vec(0, 0, 0), size=vec(d2, d2, d2), color=c2, m=m2, v=vec(v2, 0, 0))
# 產生彈簧, 起點為(-0.5*d2, 0, 0), 方向為(-L0, 0, 0)
spring = helix(pos=b2.pos + vec(-0.5*d2, 0, 0), axis=vec(-L0, 0, 0), radius=0.05, thickness=0.03)
# 繪圖部分
gd1 = graph(title="<i>E</i>-<i>t</i> plot", x=0, y=300, width=600, height=450, xtitle="<i>t</i> (s)", 
            ytitle="red: <i>K</i><sub>1</sub>, green: <i>K</i><sub>2</sub>, orange: <i>U</i>, blue: <i>E</i> (J)")
kt1 = gcurve(graph=gd1, color=c1)
kt2 = gcurve(graph=gd1, color=c2)
ut = gcurve(graph=gd1, color=color.orange)
et = gcurve(graph=gd1, color=color.blue)
gd2 = graph(title="<i>v</i>-<i>t</i> and <i>a</i>-<i>t</i> plot", x=0, y=750, width=600, height=450,
            xtitle="<i>t</i> (s)", ytitle="red: <i>v</i><sub>1</sub>, green: <i>v</i><sub>2</sub> (m/s); orange: <i>a</i><sub>1</sub>, blue: <i>a</i><sub>2</sub> (m/s<sup>2</sup>)")
vt1 = gcurve(graph=gd2, color=c1)
vt2 = gcurve(graph=gd2, color=c2)
at1 = gcurve(graph=gd2, color=color.orange)
at2 = gcurve(graph=gd2, color=color.blue)

"""
 3. 物體運動部分, 重複執行直到木塊抵達邊緣時
    由虎克定律求彈簧的回復力 force = -k*delta_x
    為了使 force 變為向量, 需要乘以 spring.axis 的單位向量
    mag(a) = a.mag => 計算向量 a 的量值
    a / mag(a) = a.norm() => 回傳向量 a 的單位向量
"""
# 印出木塊撞前速度
print("m1 =", b1.m, "m2 =", b2.m)
print(b1.v.x, b2.v.x)

while ((b2.pos.x <= xmax - d2/2) and (b1.pos.x >= xmin + d1/2)):
    rate(1000)
# 計算木塊間的距離, 更新彈簧起點位置
    dx = b2.pos.x - b1.pos.x - 0.5*d1 - 0.5*d2
    spring.pos = b2.pos + vec(-0.5*d2, 0, 0)
# 若木塊間的距離大於等於彈簧原長, 彈簧未被壓縮, 回復力 = 0, 木塊加速度 = 0
# 若木塊間的距離小於彈簧原長, 彈簧被壓縮, 計算彈簧回復力, 木塊加速度
    if(dx >= L0):
        spring.axis = vec(-L0, 0, 0)
        dL = 0
        b1.a = vec(0, 0, 0)
        b2.a = vec(0, 0, 0)
    else:
        spring.axis = vec(-dx, 0, 0)
        dL = L0 - dx
        force = vec(-k*dL, 0, 0)
        b1.a = force/b1.m
        b2.a = -force/b2.m
# 更新木塊速度、位置
    b1.v += b1.a * dt
    b2.v += b2.a * dt
    b1.pos += b1.v * dt
    b2.pos += b2.v * dt
# 計算木塊動能、系統彈性位能、力學能並作圖
    k1 = 0.5 * m1 * b1.v.mag2
    k2 = 0.5 * m2 * b2.v.mag2
    u = 0.5 * k * dL**2
    e = k1 + k2 + u
    kt1.plot(pos=(t, k1))
    kt2.plot(pos=(t, k2))
    ut.plot(pos=(t, u))
    et.plot(pos=(t, e))
# 畫v-t圖及a-t圖
    vt1.plot(pos=(t, b1.v.x))
    vt2.plot(pos=(t, b2.v.x))
    at1.plot(pos=(t, b1.a.x))
    at2.plot(pos=(t, b2.a.x))
# 更新時間
    t += dt

# 印出木塊撞後速度
print(b1.v.x, b2.v.x)

程式設計部分

由於程式 15-2 和 15-1 幾乎相同，只是加上一條理想彈簧作為緩衝，因此以下只說明兩者的不同之處。

第 17 行，設定彈簧的原長 L0 = 0.5、彈性常數 k = 5.0。
第 33 行，用 helix 物件產生彈簧，彈簧的起點位置為 b2 的左側，彈簧的軸為 vec(-L0, 0, 0)。
第 41 ~ 46 行，增加一個繪圖視窗，用來畫能量 - 時間關係圖。
第 58 行，計算兩個木塊之間的距離 dx = b2.pos.x - b1.pos.x - d1/2 - d2/2。第 59 行，更新彈簧的起點位置。
第 62 ~ 72 行：若 dx >= L0，代表彈簧沒被壓縮，將彈簧的軸設定為 vec(-L0, 0, 0)，彈簧的回復力及 b1、b2 的加速度皆為0；反之，代表彈簧被壓縮，將彈簧的軸設定為 vec(-dx, 0, 0)，計算彈簧的回復力及 b1、b2 的加速度。
第 79 ~ 86 行，計算木塊動能、系統彈性位能、力學能並作圖。