粒子同士の相互作用(2)

分子間の引力、斥力

分子間の引力、斥力を計算に加えましょう。

ここからの数式展開は長いですが、最終目標は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 となる)

まず、前ページと同じように、両者の対応を定義してやります。

t_scaleとs_scaleは、それぞれ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)を使って計算できます。
例えば10²³はpow(10, 23)、370¹²はpow(370,12)、d⁷は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 = v₀ + at。1 frameあたりなので、t = 1 frameとなりt はプログラム上は消える。)
  • 44~47行目では、42行目で定義されたaをx成分、y成分に分け、それまでの速度に加えています。相互作用は2つの粒子の対の両方に働くので、相手の粒子にも(加速度の方向を反対にして)同じ処理をしています。
  • この処理は全ての粒子対について行っています。つまり、複数の分子間で生じる相互作用は、単純に「足し合わせ」で表されると仮定しているわけです。

分子間のポテンシャルについて

以下、追加の話です。
ここで用いている分子間のポテンシャルについて少々補足します。

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。¹⁾dU / dr) = 0