Try   HackMD

三維彈性碰撞

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


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

以下共有兩個程式

  1. 17-1 兩球質量相等 (GlowScript 網站動畫連結
  2. 17-2 可以分別設定兩球質量 (GlowScript 網站動畫連結
Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →
17-1 模擬程式畫面截圖

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →
17-2 模擬程式畫面截圖

理論推導

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

m1
m2
,撞前速度分別為
v1
v2
,請推導兩球撞後的速度公式。已知質量相等時的特例為

v1=v1+(v2v1)(r1r2)|r1r2|2(r1r2)

v2=v2+(v1v2)(r2r1)|r2r1|2(r2r1)

Proof:

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

Δp,則

Δp1=Δpr1r2|r1r2|

Δp2=Δp1=Δpr2r1|r2r1|

兩球撞後動量分別為

p1=p1+Δp1  m1v1=m1v1+Δpr1r2|r1r2|

p2=p2+Δp2  m2v2=m2v2+Δpr2r1|r2r1|

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

12m1v12+12m2v22=12m1v12+12m2v22

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

m1v12+m2v22=m1[v1+Δpm1(r1r2)|r1r2|]2+m2[v2+Δpm2(r2r1)|r2r1|]2

m1v12+m2v22=m1v12+2Δpv1(r1r2)|r1r2|+(Δp)2m1(r1r2)2|r1r2|2+m2v22+2Δpv2(r2r1)|r1r2|+(Δp)2m2(r2r1)2|r1r2|2

由於

(r1r2)2=|r1r2|2,可將上式化簡為

2Δp(v1v2)(r1r2)|r1r2|+(Δp)2(1m1+1m2)=0

2(v1v2)(r1r2)|r1r2|+Δpm1+m2m1m2=0

Δp=2m1m2m1+m2(v2v1)(r1r2)|r1r2|

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

v1=v1+2m2m1+m2(v2v1)(r1r2)|r1r2|2(r1r2)

v2=v2+2m1m1+m2(v1v2)(r2r1)|r2r1|2(r2r1)


程式 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</i><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</i><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)

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

  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 圖中可看出 pxpy 皆守恆。

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →
程式17-1模擬結果:px-t

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →
程式17-1模擬結果:py-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</i><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</i><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 圖中可看出 pxpy 皆守恆。

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →
程式17-2模擬結果:px-t

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →
程式17-2模擬結果:py-t

結語

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


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