【fortran+python】シミュレーション入門2~相互作用のある粒子~

プログラミング

シミュレーション入門の第二弾です。

第一弾はこちら

今回は「相互作用のある粒子」のモデルを考えてみたいと思います。

相互作用があるので粒子が衝突したときにしっかりと跳ね返るようになります。

スポンサーリンク

相互作用のある粒子の衝突

簡単のために一次元の問題を考えます。

半径\(R=1\)の二つの粒子が衝突するとき、つまり\(|x_2-x_1|<2R\)のときに

$$
F = C_{int}(2R – |x_1-x_2|)\frac{x_1-x_2}{|x_1-x_2|}
$$

という反発力が片方の粒子に働くと考えます。

これはつまり、粒子同士の接近距離\( 2R – |x_1-x_2| \)に比例した力ということです。

この力の比例定数\(C_{int}\)を指定してあげることでどのような反発なのか決まっています。

もう一方の粒子には作用反作用の法則から\(-F\)の力が働きます。

ではこのことをコードに落とし込んでみましょう。

fortranでの計算

ではfortranで数値計算をしていきましょう。

program main
  implicit none
  integer,parameter :: N=2
  real :: x(N),vx(N)
  real :: Fx(N),Fx_next(N)
  real :: dt = 0.01
	real :: C_int = 1000.0,R = 1.0,dx
  integer :: i
  integer :: i_max = 2000
  integer :: i_rec = 10
  integer :: fout = 11

  x(1)= 1.0
  vx(1)=1.0
  Fx(1) = 0.0
  x(2)= 7.0
  vx(2)=-1.0
  Fx(2) = 0.0

  open(fout,file = "two-int-force.txt")
  write(fout,*)0.0,x(1),x(2),vx(1),vx(2)

  do i=1,i_max

   x(1)=x(1) + vx(1) * dt + 0.5 * Fx(1) * dt * dt
   Fx_next(1) = 0.0
   vx(1) = vx(1) +0.5 * (Fx(1) + Fx_next(1)) *dt

   x(2)=x(2) + vx(2) * dt + 0.5 * Fx(2) * dt * dt
   Fx_next(2) = 0.0
   vx(2) = vx(2) +0.5 * (Fx(2) + Fx_next(2)) *dt

	dx = x(1) - x(2)
 	if (abs(dx) &lt; 2.0 * R) then
 	 Fx_next(1) = C_int*(2.0 * R - abs(dx))*dx/abs(dx)
 	 Fx_next(2) = -Fx_next(1)
 	else
 	 Fx_next(1) = 0.0
 	 Fx_next(2) = 0.0
 	endif
	
	Fx(1) = Fx_next(1)
	Fx(2) = Fx_next(2)

   if(mod(i,i_rec) == 0) then
    write(fout,*) i*dt,x(1),x(2),vx(1),vx(2)
   endif
  enddo
  close(fout)
end program main

相互作用なしのときのコードに

if (abs(dx) &lt; 2.0 * R) then
 	 Fx_next(1) = C_int*(2.0 * R - abs(dx))*dx/abs(dx)
 	 Fx_next(2) = -Fx_next(1)
 	else
 	 Fx_next(1) = 0.0
 	 Fx_next(2) = 0.0
 endif

というような相互作用する命令を書いただけなのでほとんど「【fortran+python】シミュレーション入門~相互作用のない2粒子のシミュレーション~」で書いたものと同じです。

上で説明した通り、\(|x_1-x_2| < 2R\)では近接距離に比例した力が働き、それ以外では相互作用が働かないような命令を書きました。

if文で条件分岐させるところがミソですね。

ではこれが本当に動くのかpythonで試してみましょう。

pythonで確認

pythonでシミュレーションしてみましょう。

こちらのコードもほぼ一緒ですね。(変えるのはファイル名のみ)

import numpy as np 
from vpython import *


data = np.loadtxt('two-int-force.txt') # データファイルの読み込み
t = data[:, 0] 
x1 = data[:, 1]
x2 = data[:, 2]

mydisplay = canvas(autoscale=False) # 描画スペースを作る。
mydisplay.x = mydisplay.y = 0 # パソコン画面上での描画スペースの位置
mydisplay.width = mydisplay.height = 640 # 描画スペースの幅と高さ
#mydisplay.background = [1, 1, 1] # 背景色。ここでは白色
mydisplay.background = color.white
mydisplay.foreground = [0, 0, 0] # デフォルトの前景色。ここでは黒色
#mydisplay.center = [5, 0, 0] # 描画の原点を移動
mydisplay.center = vector (5, 0, 0)
mydisplay.range = 20 # 描画する座標の範囲
#particle1 = visual.sphere(pos=[x1[0], 0, 0], radius=1, color=[0, 1, 0]) # 粒子1を初期位置に描く
particle1 = sphere(pos=vector (x1[0], 0, 0), radius=1, color=vector(0, 1, 0))
#particle2 = visual.sphere(pos=[x2[0], 0, 0], radius=1, color=[1, 0, 0]) # 粒子2を初期位置に描く
particle2 = sphere(pos=vector(x2[0], 0, 0), radius=1, color=vector(1, 0, 0))
myclock = label(pos=vector(0, -10, 0), text='time is %5.1f' % t[0]) # 時計を初期時刻で描く
for k in range(1, t.size): # ループを用いてパラパラ漫画を作る
 rate(25) # パラパラ漫画をめくる速さ(フレームレート)
 particle1.pos = vector(x1[k], 0, 0) # 粒子1の位置を変更して描く
 particle2.pos = vector(x2[k], 0, 0) # 粒子2の位置を変更して描く
 myclock.text = ('time is %5.1f' % t[k])

こちらを走らせてみると

互いに接近し…↓

衝突↓

そして跳ね返ります…↓

うまくシミュレーションできました!!

タイトルとURLをコピーしました