分子間の引力、斥力
分子間の引力、斥力を計算に加えましょう。
ここからの数式展開は長いですが、最終目標はLJポテンシャルによる分子間の引力、斥力を計算し、それによって生じる加速度をpixel(長さ), frame(時間)単位で求めることです。 (少し下の方に青色で示されています。)
分子間の相互作用にはLennard-Jonesポテンシャルを用いることにします。
粒子間に働く力は、ポテンシャルをr で微分して – を付ければ得られます。
計算すると
力f がわかれば、運動方程式から物体の運動を記述できます。
a は加速度なので
v の式に順次代入していきます。
さらに質量は分子量M で、
ε をK単位で
表すことにすると
となります。
単位合わせ
さて、上式を計算すれば分子間に働く力によって分子の運動(速度)がいかに変化するかがわかります。
具体的には、以前重力の効果を生じさせるときに
vy = vy + 0.04;
というように、ステップごとにv を変化させたのと同じように分子間の力によってv を変えていけばよい(加速度を生じさせてやればよい)わけです。
ここで単位が問題となります。通常、上式にはSI単位系に合わせた、時間はs (秒)、距離はm (メートル)の単位の数値を入れていきます。しかし、コンピュータ上で分子を動かすためには、時間はframe(1ステップあたりの時間)、距離はpixel(画面上のピクセル)単位での値が必要です。(分子の速度の単位は m s-1の代わりにpixel frame-1 となる)
まず、前ページと同じように、両者の対応を定義してやります。
tscaleとsscaleは、それぞれ1 frame, 1 pixel が(今回のシミュレーションの世界の中で)何ps、何pmに相当するかを指定するパラメータで、プログラム冒頭で変数として宣言し、変更できるようにします。
frameとpixelはそれぞれ時間と距離の単位になります。上式に単位をつけた数値を代入していき、現れたs と m をframe とpixelに置き換えていってやります。
Lennard-Jonesのパラメータは窒素のものを入れてみます。
d はpixel 単位で表した粒子間の距離(の数値部分)です。
さて、現れたs と m を
を使って、(あたかも代数記号のように)置き換えていきます。
単位を整理すると、下記のようにpixel / frame単位の式が得られます。
(d はpixel 単位で表した粒子中心間の距離)
Processingでは乗数はpow(a, b)を使って計算できます。
例えば1023はpow(10, 23)、37012はpow(370,12)、d7はpow(d,7)です。
プログラムへの組み込み
さて、いよいよ分子間の相互作用をプログラム中に組み込みましょう。
処理はParticleクラス内の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(particle[i].visible>0){ //相手の粒子が存在しているとき float angle = atan2(dy, dx); if (distance < 0.9 * minDist) { //半径の和の0.9倍より距離が小さく if(distance < old_distance[i]){ //かつ距離が縮まっている場合、衝突として処理する //粒子間のなす角を算出 //速度を衝突方向の水平成分と垂直成分に分ける 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); } } //衝突処理終わり else if(distance < 10 * minDist){ //粒子間距離が粒子直径の10倍より近い時 //分子間の引力、斥力を処理 float a; //f に両者間に働く力によって生じる加速度を入れる。ヒント:距離はdistanceという変数に入っている。(pixel単位) a = 0; vx = vx - a * cos(angle); //x成分と vy = vy - a * sin(angle); //y成分に分けて速度に足す particle[i].vx = particle[i].vx + a * cos(angle); //作用反作用で、相手の粒子にも(反対方向に)加速度を生じさせる particle[i].vy = particle[i].vy + a * sin(angle); } //分子間の相互作用の処理終わり } old_distance[i] = distance; //計算した粒子間距離は保存しておく } } }
- 15行目~35行目まで
粒子間距離がσ の0.9倍以下の時
前回やったように、衝突(弾性衝突)として処理します。(後述) - 37行目~48行目まで
ここが粒子間の相互作用を処理する部分となります。- 上で導いた式を使って、42行目でa=0となっている右辺に、速度の増分(加速度; 1 pixel frame-1単位)を記入してください。
(v = v0 + at。1 frameあたりなので、t = 1 frameとなりt はプログラム上は消える。) - 44~47行目では、42行目で定義されたaをx成分、y成分に分け、それまでの速度に加えています。相互作用は2つの粒子の対の両方に働くので、相手の粒子にも(加速度の方向を反対にして)同じ処理をしています。
- この処理は全ての粒子対について行っています。つまり、複数の分子間で生じる相互作用は、単純に「足し合わせ」で表されると仮定しているわけです。
- 上で導いた式を使って、42行目でa=0となっている右辺に、速度の増分(加速度; 1 pixel frame-1単位)を記入してください。
分子間のポテンシャルについて
以下、追加の話です。
ここで用いている分子間のポテンシャルについて少々補足します。
LJポテンシャルはこのような形状です。(青)
横軸はσ、縦軸は-ε で規格化してあります。
赤は斥力項(r の-12乗の項)、緑は引力項(r の-6乗の項)です。
さて、遠い方から見ていくと、上の青線はr = 3 σ くらいまで行くとほぼ0になっていることがわかります。
上のプログラムでは38行目でσ の10倍より近い時だけ計算して、それ以遠については実質0なので、計算をしませんでした。
このような距離を「カットオフ半径」といい、これによって計算時間を短縮することができます。
分子を近づけていくと、r = 1.1 σ のあたり(正確には2の1/6乗≈1.122 σ)で青線は最小値となります。この距離より遠いと分子間には引力が働き(緑の効果が赤より大きくなる)、この距離より近いと分子間には斥力が働きます(赤の効果が緑より大きくなる)。
従って、分子の相対速度が小さい時(= 温度が低い時)、2つの分子はこの距離を中心に振動運動をするはずです。
r = σ は分子間のポテンシャルが0となる位置であり、分子同士が速い速度で衝突すると、分子はこの距離より近づきます。画面上で分子を直径σ で描いているなら、両者は「めり込む」ことがある、ということです。
計算(シミュレーション)を行う際、ここで問題が生じます。r < σ の領域では、少し近づいただけでやたらとポテンシャルが上昇します。
このプログラムは時間を(frame単位で)とびとびで計算しているため、速い相対速度で衝突すると、分子同士がめり込んだ後、ものすごい強い力で両者が反発してしまうのです。
これを避けるため、r = 0.9 σ で弾性衝突をする、という処理(15行目)を加えています。
先週やったような、r = σ で弾性衝突をする、という計算(剛体球モデル、ですね)は、ポテンシャルとして次のような仮定をしていることに相当します。
すなわち、r ≤σ のとき U = ∞, r > σ のとき U = 0 です。
2つの粒子は、どんな速い速度で衝突してもr ≤σ となることはありません。
上のプログラムではr = 0.9 σ で弾性衝突ですから
こんなポテンシャルということになります。
本当は(LJポテンシャルに従うとすると)点線のようなポテンシャルですが、r = 0.9 σ より近い時はLJ相互作用を計算せず、剛体球として処理しているということですね。
時間の刻み幅が分子の運動に対して十分に小さければ、r = 0.9 σ まで近づかず、r = 1 σ ~ 0.9 σ の間で、 LJ相互作用によって跳ね返るはずです。
追加の話のさらに追加になりますが、最初から分子がr = 0.9 σ より近い場合、(分子を生成させるとき、同じ場所でクリック連打などをしたときに起こる)
相互作用は0としていますから、より正確には
というようなポテンシャルを仮定していることになりますね。
(r < 0.9 σ では、分子間距離が近づいている場合剛体球として処理、遠ざかっている場合はF = 0。1)dU / dr) = 0
float elast = 1;
//他の粒子との衝突
if(visible>0){ //自分自身が存在しているとき
for (int i = id + 1; i < particle_number; i++) { //他の粒子との相互作用を計算
//粒子間距離の算出
float dx = others[i].x – x;
float dy = others[i].y – y;
float distance = sqrt(dx*dx + dy*dy); //粒子間距離
float minDist = others[i].r/2 + r/2; //粒子の半径の和
if(others[i].visible>0){ //相手の粒子が存在しているとき
float angle = atan2(dy, dx);
if (distance < 0.9 * minDist) { //半径の和の0.9倍より距離が小さく
if(distance < old_distance[i]){ //かつ距離が縮まっている場合、衝突として処理する
//粒子間のなす角を算出
//速度を衝突方向の水平成分と垂直成分に分ける
float v1h = vx*cos(angle)+vy*sin(angle);
float v1v = vx*sin(angle)-vy*cos(angle);
float v2h = others[i].vx*cos(angle)+others[i].vy*sin(angle);
float v2v = others[i].vx*sin(angle)-others[i].vy*cos(angle);
//垂直成分はそのまま
//水平成分は運動量保存則と弾性係数から新しい値を計算する
float nv1h = (v2h-v1h)*(1+elast)/(m/others[i].m+1)+v1h;
float nv2h = (v1h-v2h)*(1+elast)/(others[i].m/m+1)+v2h;
//速度の水平、垂直成分を画面上のx方向、y方向に書き戻す
vx = nv1h*cos(angle)+v1v*sin(angle);
vy = nv1h*sin(angle)-v1v*cos(angle);
others[i].vx = nv2h*cos(angle)+v2v*sin(angle);
others[i].vy = nv2h*sin(angle)-v2v*cos(angle);
}
} //衝突処理終わり
else if(distance < 10 * minDist){ //粒子間距離が粒子直径の10倍より近い時
//分子間の引力、斥力を処理
float a;
//f に両者間に働く力によって生じる加速度を入れる。ヒント:距離はdistanceという変数に入っている。(pixel単位)
a = kr/pow(distance, 13)-ki/pow(distance, 7);
vx = vx – a * cos(angle); //x成分と
vy = vy – a * sin(angle); //y成分に分けて速度に足す
others[i].vx = others[i].vx + a * cos(angle); //作用反作用で、相手の粒子にも(反対方向に)加速度を生じさせる
others[i].vy = others[i].vy + a * sin(angle);
} //分子間の相互作用の処理終わり
}
old_distance[i] = distance; //計算した粒子間距離は保存しておく
}
}
}
脚注
↑1 | dU / dr) = 0 |