シミュレーション入門の第二弾です。
第一弾はこちら
今回は「相互作用のある粒子」のモデルを考えてみたいと思います。
相互作用があるので粒子が衝突したときにしっかりと跳ね返るようになります。
相互作用のある粒子の衝突
簡単のために一次元の問題を考えます。
半径\(R=１\)の二つの粒子が衝突するとき、つまり\(|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) < 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) < 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])
こちらを走らせてみると
互いに接近し…↓
衝突↓
そして跳ね返ります…↓
うまくシミュレーションできました！！
