粒子同士の衝突
ここまで、粒子同士はたがいに関与せず、すり抜けてしまっていました。
いよいよ粒子同士の相互作用について考えます。 まずは粒子が剛体球だとして、跳ね返る(弾性衝突)処理を行いましょう。
粒子が衝突しているかどうかは、2つの粒子間の距離がそれぞれの半径の和より小さくなっているかどうかを調べればわかります。 この判定は、すべての粒子の間で行われなくてはなりません。 例えば粒子が100個ある場合には1回のステップごとに
for(int i=0; i<100; i++){ //全ての粒子について
for(int j=i+1; j<100; j++){ //他の粒子との距離を調べる
distance = sqrt(pow(x[i]-x[j], 2)+pow(y[i]-y[j], 2));
if(distance<r){
//衝突処理をここに書く
}
}
}
というような処理が必要になります。
sqrt(a)は a の平方根を返す関数、pow(a, b) は a の b 乗を返す関数。
2重のforループで処理をすることになりますが、2つめのループ の初期値が i+1 になっていることに注意してください。
0番の粒子については1~99番の粒子との距離をチェック、
1番の粒子については2~99番の粒子との距離をチェック、
2番の粒子については3~99番の粒子との距離をチェック、
・
・
50番の粒子については51~99番の粒子との距離をチェック、
・
・
97番の粒子と98, 99番の粒子との距離をチェック、
98番の粒子と99番の粒子との距離をチェック、
して全部の「対」(つい)のチェックが終わることになります。
チェックの回数としては
4950の粒子の対についての距離の計算が必要になります。
粒子間の距離を計算する
下記のソースファイルを実行してください。
int particle_number = 64; //粒子の最大数を指定できるようにした
Particle[] particle = new Particle[particle_number];
int j;
void setup(){
size(400, 400);
background(0);
smooth();
noStroke();
for(int i=0; i < particle_number; i++){
particle[i] = new Particle(i);
}
}
void draw(){
background(0);
for(int i=0; i<particle_number; i++){
particle[i].display();
particle[i].move();
particle[i].reflect();
particle[i].collide(); //collideを追加
}
}
void mousePressed(){ //マウスクリック時
particle[j].move_mouse_position();
j = j + 1;
if(j == particle_number){
j = 0;
}
}
class Particle{
int id; //番号
float x; //x座標
float y; //y座標
float vx; //x方向の速度
float vy; //y方向の速度
float visible; //動き始めてからの時間
float r; //直径
float m; //質量
float[] old_distance = new float[particle_number];
int R;
int G;
int B;
Particle(int iid){ //初期化
id = iid;
visible = 0;
for(int i=0; i < particle_number; i++){
old_distance[i] = width + height;
}
}
void display(){ //粒子の表示
if(visible > 0){
fill(R, G, B, visible);
ellipse( x, y, r, r);
}
}
void move(){ //粒子の移動
x = x + vx;
y = y + vy;
}
void move_mouse_position(){ //粒子を生成
float deg;
deg = random(2*PI);
x = mouseX;
y = mouseY;
vx = cos(deg);
vy = -sin(deg);
r = 40;
m = 1;
if(id==0){ //一つ目は赤く
R = 255;
G = 0;
B = 0;
}
else if(id==1){ //二つ目は青く
R = 0;
G = 0;
B = 255;
}
else{ //その他は白
R = G = B = 255;
}
visible = 180; //半透明にした
}
void reflect(){
//壁との衝突
float elast = 1;
if(x + r/2 > width){
x = width - r/2;
vx = -elast*vx;
}
if(x - r/2 < 0){
x = r/2;
vx = -elast*vx;
}
if(y + r/2 > height){
y = width - r/2;
vy = -elast*vy;
}
if(y - r/2 < 0){
y = r/2;
vy = -elast*vy;
}
}
void collide() {
//他の粒子との衝突
if(visible>0){
for (int i = id + 1; i < particle_number; i++) {
if(particle[i].visible>0){
//粒子間距離の算出
float dx = particle[i].x - x;
float dy = particle[i].y - y;
float distance = sqrt(dx*dx + dy*dy);
if( (id==0)&&(i==1) ){
print("distance " + distance + "\n");
} //if文終わり
old_distance[i] = distance;
} //if文終わり
} //for文終わり
} //if文終わり
} //collide()終わり
}
画面上をクリックすると粒子を生成するようになっています。 粒子は依然としてすり抜けますが、粒子同士の間の距離を計算する機能を加えました。赤玉と青玉の間の距離が、ソースファイルの下の部分に表示されます。
説明
- 粒子間の距離を見る処理、collide()を加えました。(117行目より)
全ての粒子対について(for(int i = id; i < particle_number; i++))
自身の位置(x, y)と他の粒子の位置(particle[i].x, paritlce[i].y)から
粒子間距離(distance)を算出しています。
ソースファイルの下にメッセージを出す命令はprintです。(119行目) - 衝突判定の処理に必要なので、クラスParticleに自分自身の番号(id)を持たせ、前回の距離(old_distance[i])を記録するようにしています。
衝突の処理
2粒子間の衝突処理を行うには、衝突が起きた(両者の距離が半径の和より小さくなった)ときに、それぞれの速度を適切な値に変更することになります。
処理後に必要なのは2つの粒子の速度で、2次元の場合、未知変数は4つです。
vx, vy, pariticle[i].vx, particle[i].vy
4つの値を決定するためには4つの方程式が必要となります。今回は以下の4つの式がそれに相当します。
- エネルギー保存則
2つの粒子の持っている運動エネルギー(mv2/2)の和は衝突前と衝突後で等しい。 - 運動量保存則
2つの粒子の持っている運動量(mv)の和は衝突前と衝突後で等しい。
運動量はベクトルなので、- 粒子1の衝突方向に垂直な向きの運動量は衝突前後で変化しない。(その方向には力が加わらないから)
- 粒子2の衝突方向に垂直な向きの運動量は衝突前後で変化しない。(同)
- 2つの粒子の衝突方向に水平な向きの運動量の和は衝突前後で変化しない。
2つの粒子が衝突前に持っていた運動量、運動エネルギーを計算し、上の4つの条件を用いることで衝突後の速度を得ることができます。
ちょっと大変なので、下記にソースコードを置きました。collide() をこれに入れ替えれば、衝突処理が行われるようになるはずです。
void collide() {
float elast = 1;
//他の粒子との衝突
if(visible>0){
for (int i = id + 1; i < particle_number; i++) {
//粒子間距離の算出
float dx = particle[i].x - x;
float dy = particle[i].y - y;
float distance = sqrt(dx*dx + dy*dy);
float minDist = particle[i].r/2 + r/2;
if (distance < minDist) { //半径の和より距離が小さく
if(particle[i].visible>0){
if(distance < old_distance[i]){ //かつ距離が縮まっている場合、衝突 (←この処理をしないと1回の衝突で2回処理をしてしまうことがある)
//粒子間のなす角を算出
float angle = atan2(dy, dx);
//速度を衝突方向の水平成分と垂直成分に分ける
float v1h = vx*cos(angle)+vy*sin(angle);
float v1v = vx*sin(angle)-vy*cos(angle);
float v2h = particle[i].vx*cos(angle)+particle[i].vy*sin(angle);
float v2v = particle[i].vx*sin(angle)-particle[i].vy*cos(angle);
//垂直成分はそのまま
//水平成分は運動量保存則と弾性係数から新しい値を計算する
float nv1h = (v2h-v1h)*(1+elast)/(m/particle[i].m+1)+v1h;
float nv2h = (v1h-v2h)*(1+elast)/(particle[i].m/m+1)+v2h;
//速度の水平、垂直成分を画面上のx方向、y方向に書き戻す
vx = nv1h*cos(angle)+v1v*sin(angle);
vy = nv1h*sin(angle)-v1v*cos(angle);
particle[i].vx = nv2h*cos(angle)+v2v*sin(angle);
particle[i].vy = nv2h*sin(angle)-v2v*cos(angle);
}
}
}
old_distance[i] = distance; //計算した粒子間距離は保存しておく
}
}
}
1回の衝突で衝突処理が2回行われることを防ぐため、old_distance[i] を使っています。
衝突処理は、距離が両者の半径の和よりも小さく、(13行目)
かつ、距離が前回よりも小さくなっている(粒子が近づいている)(15行目)ときに限り
実行しています。
処理の終了後、distance を old_distance[i] に代入し、次のターンで「前回の距離」として利用できるようにします。(37行目)
なお、visible が 0 以下のときは、その粒子が存在していないものとして、処理を行っていません。
下記は2019/1/11版
if (distance < minDist) { //半径の和より距離が小さく
if(distance < old_distance[inum]){ //かつ距離が縮まっている場合、衝突 (←この処理をしないと1回の衝突で2回処理をしてしまうことがある)
//粒子間のなす角を算出
float angle = atan2(dy, dx);
//速度を衝突方向の水平成分と垂直成分に分ける
float v1h = vx*cos(angle)+vy*sin(angle);
float v1v = vx*sin(angle)-vy*cos(angle);
float v2h = ball[inum].vx*cos(angle)+ball[inum].vy*sin(angle);
float v2v = ball[inum].vx*sin(angle)-ball[inum].vy*cos(angle);
//垂直成分はそのまま
//水平成分は運動量保存則と弾性係数から新しい値を計算する
float nv1h = (v2h-v1h)*(1+elast)/(m/ball[inum].m+1)+v1h;
float nv2h = (v1h-v2h)*(1+elast)/(ball[inum].m/m+1)+v2h;
//速度の水平、垂直成分を画面上のx方向、y方向に書き戻す
vx = nv1h*cos(angle)+v1v*sin(angle);
vy = nv1h*sin(angle)-v1v*cos(angle);
ball[inum].vx = nv2h*cos(angle)+v2v*sin(angle);
ball[inum].vy = nv2h*sin(angle)-v2v*cos(angle);
}
}
old_distance[i] = distance; //計算した粒子間距離は保存しておく
//if(distance<turnmind){
// turnmind=distance;
//}
例題
- 重力効果を有効にしたとき、分子同士がすり抜ける場合と跳ね返る場合の挙動を比較してください。
- この時点で、粒子はMaxwell-Bolzmann 分布で表される速度分布をしているはずです。(2次元だけど)
各粒子を、それぞれの速度に応じた色付けをしてみてください。(速いものは濃い赤、遅いものは青、など) - 壁で反射させる代わりに、周期境界条件を課してください。
- 重力効果を有効にした状態で、床(下の壁)を斜めにできるでしょうか。