碰撞
1. 彈性碰撞公式：
Formula for elastic collision :
\[\vec{p}_{1f}=\vec{p}_1+\vec{p}\]
\[\vec{p}_{2f}=\vec{p}_2-\vec{p}\]
\[\vec{p}=\frac{2(m_1\vec{p}_2-m_2\vec{p}_1) \cdot \hat{r}}{m_1 + m_2} \hat{r}\]
\[\hat{r}=\frac{\vec{r}_2-\vec{r}_1}{|\vec{r}_2-\vec{r}_1|}\]
2. \(m_1=m_2\)：
from visual import *
size = 0.1
scene = display(width=800, height=800, background=(0.5,0.5,0))
b1 = sphere(radius = size, color=color.yellow, make_trail=True)
b2 = sphere(radius = size, color=color.green, make_trail=True)
b1.pos = vector( 0.6 , 0.1, 0); b1.v = vector(-0.2, 0, 0)
b2.pos = vector( 0 , 0 , 0); b2.v = vector(0 , 0, 0)
v1=arrow(pos=b1.pos,axis=b1.v)
KC=0; t=0.; dt = 0.001
while True:
rate(1000)
b1.pos += b1.v*dt
b2.pos += b2.v*dt
if abs(b1.pos-b2.pos) < 2*size and dot(b1.v, b2.v) < = 0 :
v1p=b1.v-(b1.pos-b2.pos)*dot(b1.v-b2.v,b1.pos-b2.pos)/abs(b1.pos-b2.pos)**2
v2p=b2.v-(b2.pos-b1.pos)*dot(b2.v-b1.v,b2.pos-b1.pos)/abs(b2.pos-b1.pos)**2
b1.v=v1p
b2.v=v2p
if(KC==0):
v2p=arrow(); v1p=arrow()
(v1p.pos,v1p.axis)=(b1.pos,b1.v)
(v2p.pos,v2p.axis)=(b2.pos,b2.v)
KC=1
t=t+dt
if(abs(b1.pos-b2.pos)>size*10): break
3. General \(m_1, m_2\)：
from visual import *
def collision(v1, v2, p1, p2):
global m1,m2,t
rh=(p2-p1)/abs(p2-p1)
p=rh*(2.*m1*m2/(m1+m2))*dot(v2-v1,rh)
v1p=v1+p/m1; v2p=v2-p/m2
return v1p, v2p
size=0.1; m1=1.; m2=5.
scene = display(width=800, height=800, background=(0.5,0.5,0))
b1 = sphere(radius = size, color=color.yellow, make_trail=True)
b2 = sphere(radius = size, color=color.green, make_trail=True)
b1.pos=vector(0.6,0.1,0); b1.v=vector(-0.2,0,0)
b2.pos=vector(0,0,0); b2.v=vector(0,0,0)
v1=arrow(pos=b1.pos, axis=b1.v)
KC=0
t=0.
dt = 0.001
while True:
rate(1000)
b1.pos += b1.v * dt; b2.pos += b2.v * dt
if abs(b1.pos - b2.pos) < = 2*size :
b1.v, b2.v = collision(b1.v,b2.v,b1.pos,b2.pos)
if(KC==0):
v1p=arrow(pos=b1.pos, axis=b1.v)
v2p=arrow(pos=b2.pos, axis=b2.v)
KC=1
t=t+dt
if(abs(b1.pos-b2.pos)>size*10): break
在下面的影片中，請注意在螢幕最下方\(m_2\)的質量一直在改變，造成碰撞的形態也跟著改變。
4. 子彈與單擺
在下面的程式模擬當中，我們將看到一個水平射出的子彈以水平方向射出以完全非彈性碰撞的方式，直接射入靜止懸掛的單擺當中。單擺吸收了子彈的動量之後，開始進行鉛垂方向的圓週運動，最後成為週期性的單擺運動。請你參考下面的程式，將未完成的部分填入適當的程式設計以完成這個動畫模擬。
from visual import *
m=1.0; M=5.; g=9.8;k=10.; size = 0.05;L0 = 1.0
theta=radians(0.); omega0=-0; omega=omega0; alpha=-g*sin(theta)/L0
scene = display(width=800,height=800,center=(0,-L0/2,0), background=(0.5,0.5,0))
ceiling = box(length=0.2, height=0.001, width=0.2, color=color.blue)
ball = sphere(radius = size, color=color.blue,trail_type="curve",interval=2, retain=20)
string = cylinder(pos=(0,0,0), axis=(0,-1,0),radius=0.003,color=color.cyan)
bull = cylinder(radius = size/2, pos=(-1,-1,0),axis=(size*1.5,0,0),color=color.red)
bull.v=vector(20,0,0)
ball.pos = vector(L0*sin(theta), -L0*cos(theta), 0)
dt = 0.01
t=0.
while abs(bull.pos-ball.pos)>size*2:
rate(20)
bull.pos.x=bull.pos.x+bull.v.x*dt
t=t+dt
fv=m*bull.v.x/(m+M)
omega=fv/L0
while True:
rate(20)
alpha = -g*sin(theta)/L0
omega += alpha*dt
theta += omega*dt
ball.pos = vector(L0*sin(theta), -L0*cos(theta), 0)
bull.pos = vector(L0*sin(theta)-size, -L0*cos(theta), 0)
string.axis = ball.pos - string.pos
t=t+dt
5. 在空中爆炸的拋體
from visual import *
size = 0.5; Lx=40.; g=9.8; v0=20; th=radians(60.); b=0.0; m1=1.; mb=0.25; Rb=0.3
scene = display(width=800, height=800, center=(Lx/2,10,0),background=(0.5,0.5,0))
scene.forward=vector(0,0.2,-1)
ground=box(pos=(Lx/2.,0.,0.),length=Lx+4,height=0.1,width=1,color=(0,0.5,1))
ball1 = sphere(radius = size, color=color.red,make_trail=True,trail_type="points",
interval=100, retain=200)
ball1.pos = vector(0.0,0.0,0)
ball1.v = vector(v0*cos(th), v0*sin(th),0.)
t_explode=1.2
dt = 0.001
t=0.
while t < t_explode:
rate(200)
t=t+dt
ball1.pos=ball1.pos+ball1.v * dt
print "%9.4f %9.4f %9.4f" %(t,ball1.pos.x,ball1.pos.y)
if (ball1.pos.y < size) and (ball1.v.y < 0.): break
if(ball1.pos.y < 0.): break
ball1.v=ball1.v-(0,g*dt,0)+(-b*ball1.v.x*dt,-b*ball1.v.y*dt,0.)
Nbull=4
r_exp=ball1.pos
v_exp=ball1.v
p=vector(3.0,0,0)
p2=vector(0.0,2,0)
bull=[]
for i in range(Nbull):
bull.append(sphere(radius=Rb,color=(1,1,0),make_trail=True,retain=8000))
bull[i].pos=r_exp
bull[0].v=(mb*v_exp+p)/mb
bull[1].v=(mb*v_exp-p)/mb
bull[2].v=(mb*v_exp+p2)/mb
bull[3].v=(mb*v_exp-p2)/mb
print i,bull[i].pos
while t < 2.5:
rate(200)
t=t+dt
for i in range(Nbull):
bull[i].pos=bull[i].pos+bull[i].v * dt
bull[i].v=bull[i].v-(0,g*dt,0)+(-b*bull[i].v.x*dt,-b*bull[i].v.y*dt,0.)
==================================
HOMEWORK for collision
1. 在下面的連結影片中我們將同時發射9個子彈射向一個硬球，子彈將受到硬球的散射，請編寫一個程式來進行這個實驗的模擬。
同學可利用下面的程式來編寫： b1=target, bull=bullet
from visual import *
Rt=0.2; Rb=0.03; m1=1.; m2=100.; Nbull=#
def collision(v1, v2, p1, p2):
global m1,m2,t
rh=(p2-p1)/abs(p2-p1)
p=rh*(2.*m1*m2/(m1+m2))*dot(v2-v1,rh)
v1p = v1 + p/m1
v2p = v2 - p/m2
return v1p, v2p
scene = display(width=800, height=800, background=(0.5,0.5,0))
b1 = sphere(radius=Rt,color=color.yellow, make_trail=True)
b1.pos=vector(0,0,0)
b1.v=vector(0,0,0)
bull=[]
for i in range(#####):
bull.######(sphere(radius=Rb,color=(1,0,0),make_trail=True))
bull[i].###=vector(1.6,-0.3+0.08*i,0)
#######.v=######(-0.2,0,0)
print i,bull[i].pos
t=0.
dt = 0.01
while True:
rate(100)
for i in range(#####):
bull[i].### += #######.# * ##
if abs(#######.pos - ##.pos)<=Rb+Rt:
#########, #### = collision(#######.v,##.v,#######.pos,##.pos)
t=t+dt
if(t>12): break
2. 在下面的影片連結中有一個正方體的氣體容器，其中有50個單原子氣體分子(例如氦原子)。剛開始每一個原子都有相同的速度大小，但方向是隨機的。經過原子間的彈性碰撞，與容器的彈性碰撞後，氣體分子的速度不停地改變，請參考下面的程式來進行這個實驗的模擬。
from visual import *
from random import random # random module for generate random numbers
from visual.graph import *
def collision(v1, v2, p1, p2):
global m,t
m1=m; m2=m
rh=(p2-p1)/abs(p2-p1)
p=rh*(2.*m1*m2/(m1+m2))*dot(v2-v1,rh)
v1p = v1 + p/m1; v2p = v2 - p/m2
return v1p, v2p
N=50; L=((24.4E-3/(6E23))*N)**(1/3.0)/2; m=4E-3/6E23; size=310E-12
L_size = L-size; k=1.38E-23; T=298.0; vrms = (3*k*T/m)**0.5
t=0; dt=0.5E-13
atoms = []
scene = display(width=800, height=800,background=(0.2,0.2,0))
container = box(length = 2*L, height = 2*L, width = 2*L, opacity=0.2, color = color.yellow )
VA=[]
for i in range(N):
position = vector(-L_size+2*L_size*random(),-L_size+2*L_size*random(),-L_size+2*L_size*random())
if i== N-1:
atom = sphere(pos=position, radius = size, color=color.yellow, make_trail = True, retain=50)
else:
atom = sphere(pos=position, radius = size, color=color.cyan)
ra, rb = pi*random(), 2*pi*random()
atom.m, atom.v = m, vector(vrms*sin(ra)*cos(rb), vrms*sin(ra)*sin(rb), vrms*cos(ra)); atom.va=abs(atom.v)
atoms.append(atom); VA.append(atom.va)
while True:
t += dt
rate(100)
kinE=0.
for i in range(N):
atoms[i].pos += atoms[i].v * dt
atoms[i].va=abs(atoms[i].v)
VA[i]=atoms[i].va
kinE=kinE+atoms[i].va**2
KH=0
for i in range(N-1):
for j in range(i+1,N):
rij=abs(atoms[i].pos-atoms[j].pos)
if(rij