一個地球人在臺北: 相疊木塊

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

這是會出現在基礎物理2B下的例題，出現的章節甚至橫跨動量守恆、功與能量、碰撞等三章，題目敘述如下：

木塊A質量為 $2m$、木塊B質量為 $m$，兩個木塊相疊，桌面與木塊A之間沒有摩擦力，A的初速度為 $v$。兩個木塊之間的動摩擦係數為 $\mu_k$，若A的長度夠長，A、B最後能以相同的速度前進，請問：(1) A、B的末速為何？ (2) 經過多久之後兩者速度相同。

兩木塊相疊示意圖

題目有兩種版本，第一種是像上圖的版本，木塊B的初位置在木塊A的最右側且原為靜止，木塊A的初速度向右；第二種則是木塊B的初位置在木塊A的最左側且初速度向右，木塊A原為靜止。這次的目標是將運動過程中木塊A、B的速度 - 時間關係圖、能量 - 時間關係圖畫出來。

模擬程式畫面截圖

程式 13-1.相疊的木塊完全非彈性碰撞

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

"""
 VPython教學: 13-1.相疊的木塊完全非彈性碰撞
 Ver. 1: 2017/6/16
 Ver. 2: 2018/1/4   修改為 Python 3.X 版程式碼
 Ver. 3: 2018/2/26   增加 v-t 圖
 Ver. 4: 2019/9/8
 作者: 王一哲
"""

from vpython import *

"""
 1. 參數設定, 設定變數及初始值
    (1) m1 = 0.2, v1 = 0.0, m2 = 0.1, v2 = 2.0
    (2) m1 = 0.2, v1 = 1.0, m2 = 0.1, v2 = 0.0
    (3) m1 = 0.2, v1 = 0.0, m2 = 0.1, v2 = 2.5 => b2 會飛出去
    (4) m1 = 0.2, v1 = 2.0, m2 = 0.1, v2 = 0.0 => b1 到畫面右側時還沒有達到等速度
"""
d1 , h1, w1 = 1.8, 0.2, 0.2             # 下方木塊1的長度 = 1.8m, 高度 = 0.2m, 寬度 = 0.2 m
d2 , h2, w2 = 0.2, 0.2, 0.2             # 上方木塊2的長度 = 0.2m, 高度 = 0.2m, 寬度 = 0.2 m
m1, v1, c1 = 0.2, 0.0, color.red        # 下方木塊1的質量 = 0.2 kg, 初速 = 0.0 m/s, 紅色
m2, v2, c2 = 0.1, 2.0, color.green      # 上方木塊2的質量 = 0.1 kg, 初速 = 2.0 m/s, 綠色
xmax, xmin = 2.0, -2.0                  # x 軸範圍
g = 9.8                                 # 重力加速度 = 9.8 m/s^2
mu = 0.1                                # 動摩擦係數
dt = 0.0005                      # 畫面更新的時間間隔，單位為s, 原為0.001但不夠準確, 故改為0.0005
t = 0                                  # 模擬所經過的時間 ，單位為s，初始值為0
bx = 0                               # 計算 b2 初位置用的變數
i, te = 0, -1                           # 記錄 b1、b2 達到等速度經過時間用的變數

"""
 2. 畫面設定
"""
# 產生動畫視窗
scene = canvas(title="Two Blocks", width=800, height=300, center=vec(0, 0.4, 0), background=vec(0, 0.6, 0.6))
# 產生地板
floor = box(pos=vec(0, -0.5*h2, 0), size=vec(xmax - xmin, 0.05, 0.8), color=color.blue)
# 產生下方木塊 b1, 位於畫面最左側 xmin + d1/2, 初速度 v1
b1 = box(pos=vec(xmin + d1/2, 0, 0), size=vec(d1, h1, w1), color = c1, v=vec(v1, 0, 0))
# 產生上方木塊 b2, 若 v2 >= v1 則位於畫面最左側 xmin + d2/2; 若 v2 < v1, 則位於 b1 最右側 xmin + d1 - d2/2, 初速度 v2
if(v2 >= v1): bx = xmin + 0.5*d2
else: bx = xmin + d1 - 0.5*d2
b2 = box(pos=vec(bx, h1, 0), size=vec(d2, h2, w2), color=c2, v=vec(v2, 0, 0))
# 繪圖部分
gd = 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>, blue: <i>E</i> (J)")
kt1 = gcurve(graph=gd, color=c1)
kt2 = gcurve(graph=gd, color=c2)
et = gcurve(graph=gd, color=color.blue)
gd2 = graph(title="<i>v</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)")
vt1 = gcurve(graph=gd2, color=c1)
vt2 = gcurve(graph=gd2, color=c2)

"""
 3. 物體運動部分, 重複執行直到: (1) b1 抵達邊緣時, (2) b2 抵達 b1 邊緣時
"""
while((b1.pos.x <= xmax - d1/2) and (b2.pos.x + d2/2 <= b1.pos.x + d1/2 + 0.001)):
#動畫執行頻率, 每秒500次
    rate(500)
#由 b1, b2 速度判斷動摩擦力量值及加速度方向, 
    if(b2.v.x > b1.v.x):
        force = mu * m2 * g
        b1.a = vec(force / m1, 0, 0)
        b2.a = vec(-force / m2, 0, 0)
    elif(b2.v.x < b1.v.x):
        force = mu * m2 * g
        b1.a = vec(-force / m1, 0, 0)
        b2.a = vec(force / m2, 0, 0)
    else:
        force = 0
        b1.a = vec(0, 0, 0)
        b2.a = vec(0, 0, 0)
# 記錄 b1、b2 達到等速度需要的時間
    if(abs(b2.v.x - b1.v.x) < 0.0005 and i == 0):
        te = t
        i = 1     
#計算 b1, b2 的速度, 位置
    b1.v += b1.a * dt
    b2.v += b2.a * dt
    b1.pos += b1.v * dt
    b2.pos += b2.v * dt
# 計算 b1, b2 動能及系統力學能, 畫能量 - 時間關係圖
    k1 = 0.5 * m1 * mag2(b1.v)
    k2 = 0.5 * m2 * mag2(b2.v)
    e = k1 + k2
    kt1.plot(pos = (t, k1))            
    kt2.plot(pos = (t, k2))
    et.plot(pos = (t, e))
# 畫速度 - 時間關係圖
    vt1.plot(pos = (t, b1.v.x))            
    vt2.plot(pos = (t, b2.v.x))
#更新時間
    t += dt
#結束 while 迴圈, 印出末速及 b1、b2 達到等速度需要的時間
print("v1 = ", b1.v.x)
print("v2 = ", b2.v.x)
print("te = ", te)
print("end")

參數設定

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

畫面設定

大部分的功能在之前的動畫中都已經使用過了，應該只需要解釋比較特別的部分：由 v1、v2 決定 b2 的初位置。當 v2 >= v1 時，b2 的初位置要在 b1 的左側，否則在執行動畫之後 b2 會立刻飛出 b1；反之，當 v2 < v1 時，b2 的初位置要在 b1 的右側，因此需要以下兩行程式碼計算 b2 的初位置：

if(v2 >= v1): bx = xmin + d2/2

else: bx = xmin + d1 - d2/2

物體運動

為了使動畫重覆執行，直到以下 2 種情況發生為止：
1. b1 抵達畫面邊緣時
2. b2 抵達 b1 邊緣時
因此 while 迴圈中的條件設定為
```
((b1.pos.x <= xmax - d1/2) and (b2.pos.x + d2/2 <= b1.pos.x + d1/2 + 0.001)
```
第二個條件中加上 0.001 是為了避免當 v1 > v2 時，動畫一開始執行就停止運作。
b1、b2 之間的動摩擦力方向與當時的速度有關，共有 3 種狀況：
1. 當 b2.v.x > b1.v.x 時，動摩擦力對 b1 的作用方向向右、對 b2 的作用方向向左。
2. 當 b2.v.x < b1.v.x 時，動摩擦力對 b1 的作用方向向左、對 b2 的作用方向向右。
3. 當 b2.v.x = b1.v.x 時，動摩擦力為 0。
因此需要第 60 ~ 71 行的程式碼，利用 if … elif … else 決定動摩擦力的量值及方向。其實 else 的部分可以省略，因為在程式中 b1.v.x 和 b2.v.x 都是浮點數，幾乎不可能相等。
為了記錄 b1、b2 達到等速度經過的時間，需要第 73 ~ 75 行程式碼。由於 b1.v.x 和 b2.v.x 幾乎不可能相等，因此條件設定為abs(b2.v.x - b1.v.x) < 0.0005。由於只需要記錄 b1.v.x 和 b2.v.x 幾乎相等的第一個時刻，因此條件中再加上 i == 0，並在執行 if 中的程式碼時設定 i = 1，就不會再執行第二次。
計算木塊的加速度，更新速度、位置，計算木塊的動能，畫速度 - 時間關係圖及能量 - 時間關係圖，更新時間。
最後印出木塊末速及 b1、b2 達到等速度需要的時間。

模擬結果

b1、b2 之間的動摩擦力與加速度量值分別為

$$ f_k = \mu_k N = \mu_k m_2 g~~~~~a_1 = \frac{\mu_k m_2 g}{m_1}~~~~~a_2 = \mu_k g $$

可以由上式及動量守恆計算木塊的末速及達成等速度需要的時間。我測試了 4 種不同的條件：

m1 = 0.2, v1 = 0.0, m2 = 0.1, v2 = 2.0
理論值：v1’ = v2’ = 2/3, te ≈ 1.360544
模擬結果：v1’ = 0.6668900000000263, v2’ = 0.6662199999997995, te = 1.3604999999999057

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

條件1的模擬結果：能量 - 時間關係圖

m1 = 0.2, v1 = 1.0, m2 = 0.1, v2 = 0.0
理論值：v1’ = v2’ = 2/3, te ≈ 0.680272
模擬結果：v1’ = 0.666554999999931, v2’ = 0.6668899999999899, te = 0.6799999999999805