三維彈性碰撞

# 三維彈性碰撞 > 作者：王一哲 > 日期：2018/5/6 這個程式主要是參考臺大物理系石明豐教授的講義〈[VPhysics大一課程：碰撞](https://drive.google.com/file/d/0B7pfgC0XDvVHWnhCTHJKQjQ4eG8/view)〉，但是將其中的程式碼改寫成 python 3.X 版的格式。在寫好這個程式之後，可以將它用來模擬理想氣體分子之間的碰撞，作出分子數量 - 速率分布圖，但由於這個程式較為複雜，請參考 VPython 範例程式 "[A hard-sphere gas](http://www.glowscript.org/#/user/GlowScriptDemos/folder/Examples/program/HardSphereGas-VPython)"。以下共有兩個程式 1. 17-1 兩球質量相等　（[GlowScript 網站動畫連結](http://www.glowscript.org/#/user/yizhe/folder/Public/program/17-13Dcollision)） 2. 17-2 可以分別設定兩球質量　（[GlowScript 網站動畫連結](http://www.glowscript.org/#/user/yizhe/folder/Public/program/17-23Dcollision)） <img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://imgur.com/cjvP7d3.png"> <div style="text-align:center">17-1 模擬程式畫面截圖</div> <img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://imgur.com/vSR8t3y.png"> <div style="text-align:center">17-2 模擬程式畫面截圖</div> ## 理論推導假設空間中有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 （[取得程式碼](https://github.com/YiZheWangTw/VPythonTutorial/blob/master/17.%E4%B8%89%E7%B6%AD%E5%BD%88%E6%80%A7%E7%A2%B0%E6%92%9E/17-1_3D_collision.py)） ```python= """ 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="px-t plot", width=600, height=450, x=0, y=600, xtitle="t (s)", ytitle="blue: px1, red: px2, green: px (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="py-t plot", width=600, height=450, x=0, y=1050, xtitle="t (s)", ytitle="blue: py1, red: py2, green: py (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 計算撞後速度，函式中的內容 ```python 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) ``` 就是將理論推導中的公式用程式碼寫出來，在此應用了兩種向量的計算 1. **dot(a, b)**：將向量 a、b 取內積 2. **mag2(a) = a.mag2**：計算向量 a 量值的平方 ### 物體運動部分 1. 更新小球位置。 2. 繪製小球 v-t 圖。 3. 若 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 圖中可看出 px、py 皆守恆。 <img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://imgur.com/b2GeWLa.png"> <div style="text-align:center">程式17-1模擬結果：px-t 圖</div> <img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://imgur.com/qi0NBrn.png"> <div style="text-align:center">程式17-1模擬結果：py-t 圖</div> ## 程式 17-2.三維彈性碰撞, m1 ≠ m2 （[取得程式碼](https://github.com/YiZheWangTw/VPythonTutorial/blob/master/17.%E4%B8%89%E7%B6%AD%E5%BD%88%E6%80%A7%E7%A2%B0%E6%92%9E/17-2_3D_collision.py)） ```python= """ 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 (m1 ≠ m2)", 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="px-t plot", width=600, height=450, x=0, y=600, xtitle="t (s)", ytitle="blue: px1, red: px2, green: px (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="py-t plot", width=600, height=450, x=0, y=1050, xtitle="t (s)", ytitle="blue: py1, red: py2, green: py (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 圖中可看出 px、py 皆守恆。 <img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://imgur.com/1H8nQah.png"> <div style="text-align:center">程式17-2模擬結果：px-t 圖</div> <img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://imgur.com/qBScZWa.png"> <div style="text-align:center">程式17-2模擬結果：py-t 圖</div> ## 結語接下來可以將這個程式應用於理想氣體速率分布的模擬，在 VPython 範例 "[A hard-sphere gas](http://www.glowscript.org/#/user/GlowScriptDemos/folder/Examples/program/HardSphereGas-VPython/edit)" 當中沒有採用這條公式，所以程式碼會比較複雜一點，有興趣的同學可以研究看看。 ## VPython官方說明書 1. **canvas**: http://www.glowscript.org/docs/VPythonDocs/canvas.html 2. **sphere**: http://www.glowscript.org/docs/VPythonDocs/sphere.html 3. **graph**: http://www.glowscript.org/docs/VPythonDocs/graph.html --- ###### tags:`VPython`

Read more

斜向抛射

Python 串列 (List)

使用 Microsoft Excel 繪製帶有誤差槓的 XY 散布圖

終端速度