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

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

這個程式主要是參考臺大物理系石明豐教授的講義〈VPhysics大一課程：碰撞〉，但是將其中的程式碼改寫成 python 3.X 版的格式。在寫好這個程式之後，可以將它用來模擬理想氣體分子之間的碰撞，作出分子數量 - 速率分布圖，但由於這個程式較為複雜，請參考 VPython 範例程式 "A hard-sphere gas"。

以下共有兩個程式

17-1 兩球質量相等　（GlowScript 網站動畫連結）
17-2 可以分別設定兩球質量　（GlowScript 網站動畫連結）

17-1 模擬程式畫面截圖

17-2 模擬程式畫面截圖

理論推導

假設空間中有2個大小相等的彈性球體，其質量分別為

$m_1$ 、

$m_2$ ，撞前速度分別為

$v_1$ 、

$v_2$ ，請推導兩球撞後的速度公式。已知質量相等時的特例為

$\mathbf{v_1'} = \mathbf{v_1} + \frac{(\mathbf{v_2} - \mathbf{v_1}) \cdot (\mathbf{r_1} - \mathbf{r_2})}{|\mathbf{r_1} - \mathbf{r_2}|^2}(\mathbf{r_1} - \mathbf{r_2})$

$\mathbf{v_2'} = \mathbf{v_2} + \frac{(\mathbf{v_1} - \mathbf{v_2}) \cdot (\mathbf{r_2} - \mathbf{r_1})}{|\mathbf{r_2} - \mathbf{r_1}|^2}(\mathbf{r_2} - \mathbf{r_1})$

Proof:

假設兩球相撞時的動量變化量值為

$\Delta p$ ，則

$\mathbf{\Delta p_1} = \Delta p \frac{\mathbf{r_1} - \mathbf{r_2}}{|\mathbf{r_1} - \mathbf{r_2}|}$

$\mathbf{\Delta p_2} = -\mathbf{\Delta p_1} = \Delta p \frac{\mathbf{r_2} - \mathbf{r_1}}{|\mathbf{r_2} - \mathbf{r_1}|}$

兩球撞後動量分別為

$\mathbf{p_1'} = \mathbf{p_1} + \mathbf{\Delta p_1} ~\Rightarrow~ m_1 \mathbf{v_1'} = m_1 \mathbf{v_1} + \Delta p \frac{\mathbf{r_1} - \mathbf{r_2}}{|\mathbf{r_1} - \mathbf{r_2}|}$

$\mathbf{p_2'} = \mathbf{p_2} + \mathbf{\Delta p_2} ~\Rightarrow~ m_2 \mathbf{v_2'} = m_2 \mathbf{v_2} + \Delta p \frac{\mathbf{r_2} - \mathbf{r_1}}{|\mathbf{r_2} - \mathbf{r_1}|}$

由於兩球間為彈性碰撞，碰撞前後動能沒有損失，因此

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

將 v1' 、v2' 代入上式並同乘以2

$m_1 v_1^2 + m_2 v_2^2 = m_1 \left[ \mathbf{v_1} + \frac{\Delta p}{m_1} \frac{(\mathbf{r_1} - \mathbf{r_2})}{|\mathbf{r_1} - \mathbf{r_2}|}\right ]^2 + m_2 \left[ \mathbf{v_2} + \frac{\Delta p}{m_2} \frac{(\mathbf{r_2} - \mathbf{r_1})}{|\mathbf{r_2} - \mathbf{r_1}|}\right ]^2$

$m_1 v_1^2 + m_2 v_2^2 = m_1 v_1^2 + 2 \Delta p \frac{\mathbf{v_1} \cdot (\mathbf{r_1} - \mathbf{r_2})}{|\mathbf{r_1} - \mathbf{r_2}|} + \frac{(\Delta p)^2}{m_1} \frac{(\mathbf{r_1} - \mathbf{r_2})^2}{|\mathbf{r_1} - \mathbf{r_2}|^2} \\+ m_2 v_2^2 + 2 \Delta p \frac{\mathbf{v_2} \cdot (\mathbf{r_2} - \mathbf{r_1})}{|\mathbf{r_1} - \mathbf{r_2}|} + \frac{(\Delta p)^2}{m_2} \frac{(\mathbf{r_2} - \mathbf{r_1})^2}{|\mathbf{r_1} - \mathbf{r_2}|^2}$

由於

$(\mathbf{r_1} - \mathbf{r_2})^2 = |\mathbf{r_1} - \mathbf{r_2}|^2$ ，可將上式化簡為

$2 \Delta p \frac{(\mathbf{v_1} - \mathbf{v_2}) \cdot (\mathbf{r_1} - \mathbf{r_2})}{|\mathbf{r_1} - \mathbf{r_2}|}+(\Delta p)^2 \left( \frac{1}{m_1} + \frac{1}{m_2} \right) = 0$

$\frac{2(\mathbf{v_1} - \mathbf{v_2}) \cdot (\mathbf{r_1} - \mathbf{r_2})}{|\mathbf{r_1} - \mathbf{r_2}|} + \Delta p \frac{m_1 + m_2}{m_1 m_2} = 0$

$\Delta p = \frac{2 m_1 m_2}{m_1 + m_2} \frac{(\mathbf{v_2} - \mathbf{v_1}) \cdot (\mathbf{r_1} - \mathbf{r_2})}{|\mathbf{r_1} - \mathbf{r_2}|}$

將代入最上面2式可得撞後速度

$\mathbf{v_1'} = \mathbf{v_1} + \frac{2m_2}{m_1 + m_2}\frac{(\mathbf{v_2} - \mathbf{v_1}) \cdot (\mathbf{r_1} - \mathbf{r_2})}{|\mathbf{r_1} - \mathbf{r_2}|^2}(\mathbf{r_1} - \mathbf{r_2})$

$\mathbf{v_2'} = \mathbf{v_2} + \frac{2m_1}{m_1 + m_2}\frac{(\mathbf{v_1} - \mathbf{v_2}) \cdot (\mathbf{r_2} - \mathbf{r_1})}{|\mathbf{r_2} - \mathbf{r_1}|^2}(\mathbf{r_2} - \mathbf{r_1})$

程式 17-1.三維彈性碰撞, m1 = m2

取得程式碼

"""
 VPython教學: 17-1.三維彈性碰撞, m1=m2
 Ver. 1: 2018/3/4
 Ver. 2: 2019/9/14 改成畫 px-t 及 py-t 圖
 作者: 王一哲
"""
from vpython import *

"""
 1. 參數設定, 設定變數及初始值
"""
m1, r1, c1, v1 = 1, 1, color.blue, vec(8, 0, 0)     # 小球1質量, 半徑, 顏色, 初速
m2, r2, c2, v2 = 1, 1, color.red, vec(0, 0, 0)      # 小球2質量, 半徑, 顏色, 初速
L, t, dt = 10, 0, 0.001    # 畫面邊長, 時間, 時間間隔

"""
 2. 畫面設定
"""
# 產生動畫視窗及小球
scene = canvas(title="3D Collision", width=600, height=600, x=0, y=0, background=vec(0, 0.6, 0.6), range=L)
b1 = sphere(pos=vec(-L+r1, r1*(2/3), 0), m=m1, v=v1, radius=r1, color=c1, make_trail=True)
b2 = sphere(pos=vec(0, 0, 0), m=m2, v=v2, radius=r2, color=c2, make_trail=True)
# px-t plot
gd1 = graph(title="<i>p<sub>x</sub></i>-<i>t</i> plot", width=600, height=450, x=0, y=600, xtitle="<i>t</i> (s)",
            ytitle="blue: <i>p</i><sub><i>x</i>1</sub>, red: <i>p</i><sub><i>x</i>2</sub>, green: <i>p<sub>x</sub></i> (kg m/s)")
px1 = gcurve(graph=gd1, color=c1)
px2 = gcurve(graph=gd1, color=c2)
px = gcurve(graph=gd1, color=color.green)
# py-t plot
gd2 = graph(title="<i>p<sub>y</sub></i>-<i>t</i> plot", width=600, height=450, x=0, y=1050, xtitle="<i>t</i> (s)",
            ytitle="blue: <i>p</i><sub><i>y</i>1</sub>, red: <i>p</i><sub><i>y</i>2</sub>, green: <i>p<sub>y</sub></i> (kg m/s)")
py1 = gcurve(graph=gd2, color=c1)
py2 = gcurve(graph=gd2, color=c2)
py = gcurve(graph=gd2, color=color.green)
# 計算撞後速度的函式
def af_col_v(v1, v2, x1, x2):
    v1_prime = v1 + dot((v2 - v1), (x1 - x2)) / mag2(x1 - x2) * (x1 - x2)
    v2_prime = v2 + dot((v1 - v2), (x2 - x1)) / mag2(x2 - x1) * (x2 - x1)
    return (v1_prime, v2_prime)

"""
 3. 物體運動部分, 小球到達畫面邊緣時停止運作
"""
# 印出撞前動能
K1 = 0.5*b1.m*b1.v.mag2
K2 = 0.5*b2.m*b2.v.mag2
print("K1 =", K1, "K2 =", K2, "K =", K1 + K2)

while(abs(b1.pos.x) < L and abs(b1.pos.y) < L and abs(b2.pos.x) < L and abs(b2.pos.y) < L):
    rate(500)
# 更新小球位置
    b1.pos += b1.v*dt
    b2.pos += b2.v*dt
# 繪製小球 p-t 圖
    px1.plot(pos=(t, b1.m*b1.v.x))
    px2.plot(pos=(t, b2.m*b2.v.x))
    px.plot(pos=(t, b1.m*b1.v.x + b2.m*b2.v.x))
    py1.plot(pos=(t, b1.m*b1.v.y))
    py2.plot(pos=(t, b2.m*b2.v.y))
    py.plot(pos=(t, b1.m*b1.v.y + b2.m*b2.v.y))
# 若 b1、b2 相撞則計算撞後速度並重新指定給 v1, v2
    if(mag(b1.pos - b2.pos) <= r1 + r2 and dot((b1.pos - b2.pos), (b1.v - b2.v)) <=0):
        b1.v, b2.v = af_col_v(b1.v, b2.v, b1.pos, b2.pos)
        cm = sphere(pos=(b1.pos + b2.pos)/2, radius=r1/5, color=color.yellow)
# 更新時間
    t += dt

# 印出撞後動能
K1 = 0.5*b1.m*b1.v.mag2
K2 = 0.5*b2.m*b2.v.mag2
print("K1 =", K1, "K2 =", K2, "K =", K1 + K2)

參數設定

在此設定變數為小球的半徑、質量、顏色、初速度，畫面邊長、時間、時間間隔，對應的變數名稱請參考程式碼。

畫面設定

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

自訂函式

自訂函式 af_col_v 計算撞後速度，函式中的內容

v1_prime = v1 + dot((v2 - v1), (x1 - x2)) / mag2(x1 - x2) * (x1 - x2)
v2_prime = v2 + dot((v1 - v2), (x2 - x1)) / mag2(x2 - x1) * (x2 - x1)

就是將理論推導中的公式用程式碼寫出來，在此應用了兩種向量的計算

dot(a, b)：將向量 a、b 取內積
mag2(a) = a.mag2：計算向量 a 量值的平方

物體運動部分

更新小球位置。
繪製小球 v-t 圖。
若 mag(b1.pos - b2.pos) <= r1 + r2 且 dot((b1.pos - b2.pos), (b1.v - b2.v)) <=0 代表 b1、b2 相撞，將 b1.v、b2.v、b1.pos、b2.pos 代入自訂函式 af_col_v 中計算撞後速度，再重新指定給 b1.v、b2.v。

模擬結果

假設 m1 = m2 = 1, v1 = (8, 0, 0), v2 = (0, 0, 0)，撞前動能 K1 = 32.0, K2 = 0.0, K = 32.0，撞前動能 K1 = 3.574460479870966, K2 = 28.425539520129036, K = 32.0，由 p-t 圖中可看出 p_x、p_y 皆守恆。

程式17-1模擬結果：p_x-t 圖

程式17-1模擬結果：p_y-t 圖

程式 17-2.三維彈性碰撞, m1 ≠ m2

取得程式碼

"""
 VPython教學: 17-2.三維彈性碰撞, m1 ≠ m2
 Ver. 1: 2018/3/4
 Ver. 2: 2019/9/14 改成畫 px-t 及 py-t 圖
 作者: 王一哲
"""
from vpython import *

"""
 1. 參數設定, 設定變數及初始值
"""
m1, r1, c1, v1 = 1, 1, color.blue, vec(8, 0, 0)     # 小球1質量, 半徑, 顏色, 初速
m2, r2, c2, v2 = 2, 1, color.red, vec(0, 0, 0)      # 小球2質量, 半徑, 顏色, 初速
L, t, dt = 10, 0, 0.001    # 畫面邊長, 時間, 時間間隔

"""
 2. 畫面設定
"""
# 產生動畫視窗及小球
scene = canvas(title="3D Collision (<i>m</i><sub>1</sub> ≠ <i>m</i><sub>2</sub>)", width=600, height=600,
               x=0, y=0, background=vec(0, 0.6, 0.6), range=L)
b1 = sphere(pos=vec(-L+r1, r1*(2/3), 0), m=m1, v=v1, radius=r1, color=c1, make_trail=True)
b2 = sphere(pos=vec(0, 0, 0), m=m2, v=v2, radius=r2, color=c2, make_trail=True)
# px-t plot
gd1 = graph(title="<i>p<sub>x</sub></i>-<i>t</i> plot", width=600, height=450, x=0, y=600, xtitle="<i>t</i> (s)",
            ytitle="blue: <i>p</i><sub><i>x</i>1</sub>, red: <i>p</i><sub><i>x</i>2</sub>, green: <i>p<sub>x</sub></i> (kg m/s)")
px1 = gcurve(graph=gd1, color=c1)
px2 = gcurve(graph=gd1, color=c2)
px = gcurve(graph=gd1, color=color.green)
# py-t plot
gd2 = graph(title="<i>p<sub>y</sub></i>-<i>t</i> plot", width=600, height=450, x=0, y=1050, xtitle="<i>t</i> (s)",
            ytitle="blue: <i>p</i><sub><i>y</i>1</sub>, red: <i>p</i><sub><i>y</i>2</sub>, green: <i>p<sub>y</sub></i> (kg m/s)")
py1 = gcurve(graph=gd2, color=c1)
py2 = gcurve(graph=gd2, color=c2)
py = gcurve(graph=gd2, color=color.green)
# 計算撞後速度的函式
def af_col_v(m1, m2, v1, v2, x1, x2):
    v1_prime = v1 + (2*m2)/(m1 + m2) * dot((v2 - v1), (x1 - x2)) / mag2(x1 - x2) * (x1 - x2) 
    v2_prime = v2 + (2*m1)/(m1 + m2) * dot((v1 - v2), (x2 - x1)) / mag2(x2 - x1) * (x2 - x1)
    return (v1_prime, v2_prime)

"""
 3. 物體運動部分, 小球到達畫面邊緣時停止運作
"""
# 印出撞前動能
K1 = 0.5*b1.m*b1.v.mag2
K2 = 0.5*b2.m*b2.v.mag2
print("K1 =", K1, "K2 =", K2, "K =", K1 + K2)

while(abs(b1.pos.x) < L and abs(b1.pos.y) < L and abs(b2.pos.x) < L and abs(b2.pos.y) < L):
    rate(500)
# 更新小球位置
    b1.pos += b1.v*dt
    b2.pos += b2.v*dt
# 繪製小球 p-t 圖
    px1.plot(pos=(t, b1.m*b1.v.x))
    px2.plot(pos=(t, b2.m*b2.v.x))
    px.plot(pos=(t, b1.m*b1.v.x + b2.m*b2.v.x))
    py1.plot(pos=(t, b1.m*b1.v.y))
    py2.plot(pos=(t, b2.m*b2.v.y))
    py.plot(pos=(t, b1.m*b1.v.y + b2.m*b2.v.y))
# 若 b1、b2 相撞則計算撞後速度並重新指定給 v1, v2
    if(mag(b1.pos - b2.pos) <= r1 + r2 and dot((b1.pos - b2.pos), (b1.v - b2.v)) <=0):
        b1.v, b2.v = af_col_v(b1.m, b2.m, b1.v, b2.v, b1.pos, b2.pos)
        cm = sphere(pos=(b1.pos + b2.pos)/2, radius=r1/5, color=color.yellow)
# 更新時間
    t += dt

# 印出撞後動能
K1 = 0.5*b1.m*b1.v.mag2
K2 = 0.5*b2.m*b2.v.mag2
print("K1 =", K1, "K2 =", K2, "K =", K1 + K2)

程式設計部分

程式 17-2 和 17-1 幾乎相同，只是將計算撞後速度的公式改為有質量的版本。

模擬結果

假設 m1 1, m2 = 2, v1 = (8, 0, 0), v2 = (0, 0, 0)，撞前動能 K1 = 32.0, K2 = 0.0, K = 32.0，撞前動能 K1 = 6.7328537598853035, K2 = 25.2671462401147, K = 32.0，由 p-t 圖中可看出 p_x、p_y 皆守恆。

程式17-2模擬結果：p_x-t 圖

程式17-2模擬結果：p_y-t 圖

結語

接下來可以將這個程式應用於理想氣體速率分布的模擬，在 VPython 範例 "A hard-sphere gas" 當中沒有採用這條公式，所以程式碼會比較複雜一點，有興趣的同學可以研究看看。

VPython官方說明書

HackMD 版本連結：https://hackmd.io/@yizhewang/H1-59qBMm

一個地球人在臺北

熱門文章

2018年5月6日星期日

三維彈性碰撞

理論推導

Proof:

程式 17-1.三維彈性碰撞, m1 = m2

參數設定

畫面設定

自訂函式

物體運動部分

模擬結果

程式 17-2.三維彈性碰撞, m1 ≠ m2

程式設計部分

模擬結果

結語

VPython官方說明書

沒有留言:

張貼留言

熱門文章

2018年5月6日 星期日

三維彈性碰撞

理論推導

Proof:

程式 17-1.三維彈性碰撞, m1 = m2

參數設定

畫面設定

自訂函式

物體運動部分

模擬結果

程式 17-2.三維彈性碰撞, m1 ≠ m2

程式設計部分

模擬結果

結語

VPython官方說明書

沒有留言:

張貼留言

2018年5月6日星期日