仮想世界の温度を制御する

系の温度を制御する

相互作用を導入すると、分子間の距離によって系の温度が変わるようになります。
相互作用なし(弾性衝突のみ)だと、エネルギー保存則によって分子の運動エネルギーの和は一定に保たれるので、系の温度も一定です。
しかし、相互作用を導入すると、運動エネルギーに分子間の相互作用による位置エネルギー(=ポテンシャルエネルギー)が加わるので、運動エネルギーは一定ではなくなります。

(運動エネルギーの和)+(位置エネルギーの和) = 一定

系の温度を制御する、たとえば系の温度を下げてやるには、分子の運動エネルギーを小さく、すなわち分子の速度を下げてやればよいことがわかります。

系に含まれる分子の運動エネルギーから算出される温度が目標の温度の倍であれば、すべての分子の運動エネルギーを半分にしてやればよいのです。

それまでの速度をv 、新しい速度をv‘ とすると、温度を半分(1/2)にするには

 

 \frac{1}{2}mv'^2 = \frac{1}{2} \cdot \frac{1}{2}mv^2 

より

 v' = \sqrt{\frac{1}{2}} \cdot v 

とすればよく、温度をk 倍にするには

 v' = \sqrt{k} \cdot v 

と、各分子の速度を修正すればよいことになります。

これはもちろん近似です。
実際には分子の速度は壁の分子との衝突などによって個々に変化し、
長い時間をかけて一定の分布(その温度でのMaxwell-Bolzmann分布)になってゆくでしょう。
(例えば熱いガスを冷たい容器に入れた場合… 速い速度を持った分子が、ゆっくりと格子の所定位置で振動している壁の分子に衝突し、運動エネルギーを壁の分子に奪われ、比較的ゆっくりとした速度で跳ね返る)

いきなり目標温度に分子の速度を修正するのではなく、比較的ゆっくりと修正することにしましょう。
現在の(系の分子の平均運動エネルギーから算出した)温度をTave、目標温度をTobjとし、1回の処理で1割づつ目標温度に近づける処理は
(例えばTave = 200 K, Tobj = 100 Kなら、200K→190K→181 K→172.9K→…)

 

 \displaystyle v' = \sqrt{\frac{10 \cdot T_{\rm obj}}{9 \cdot T_{\rm ave}+T_{\rm obj}}}\cdot v 

になります。

(2010.1.27 逆数になっていたのを修正しました。thanks→平岡君)

この処理を10ステップに1回実行することにしましょう。

これまでのプログラムに、下記のような処理を加えます。

float temp = 50; //目標温度
int step = 0; //ステップ数

void draw(){
  step++;
  background(0);
  for(int i=0; i<particle_number; i++){
    particle[i].display();
    particle[i].move();
    particle[i].reflect();
    particle[i].collide();
  }
  float average;

  //ここに系の温度averageを求める処理を書く

  if(step%10==0){
    for(int i=0;i<particle_number;i++){
      if(particle[i].visible>0){
        particle[i].vx *= 1; //ここに速度を変える処理を書く
        particle[i].vy *= 1; //vx, vyを同じ倍率で変更すればよい
      }
    }
  }
}
  • 1行目 目標温度を変数tempに設定します。
  • 2行目 ステップ数を数えるための変数stepを設定します。(5行目で1回実行するたびに1づつ増やしている)
  • 17行目
    %は割り算の余りを出す演算子です。(整数型の変数にのみ有効)
    ステップ数を10で割った時の余りが0のとき、すなわち割り切れるときだけ{}内を実行します。
  • 18、19行目
    分子のうち存在しているもの(visible>0)全てについて速度変更の処理を行います。

画面に温度を表示するルーチンはそのまま残しておくとよいでしょう。
系の温度は50 Kになったでしょうか。
この温度であれば、窒素は固体になるはずです。