投稿日: 投稿者: simu-kun

【第５話】シリンダー内の粒子の運動計算アルゴリズム【ルンゲクッタで行こう！⑤】

############################################################################
#　シリンダーの中の１粒子運動シミュレーション
############################################################################
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

#図全体
fig = plt.figure(figsize=(10, 10))
#全体設定
plt.rcParams['font.family'] = 'Times New Roman' #フォント
plt.rcParams['font.size'] = 16 #フォントサイズ
plt.rcParams["mathtext.fontset"] = 'cm' #数式用フォント
#カラーリストの取得
colors = plt.rcParams['axes.prop_cycle'].by_key()['color']

#################################################################
## 物理系の設定
#################################################################
#質量
m = 1.0
M = 10.0
radius = 0.05
#箱の長さ
L = 2.0
#重力加速度
g = np.array([0.0, 0.0, -10.0])

def F ( t, m, r_t, v_t, fc ):
	return m * g + fc

#初期位置と速度
r0 = np.array([0.0 , 0.0, 0.0])
v0 = np.array([0.5, 0.0, 2.5])
R0 = np.array([0.0 , 0.0, L/2])
V0 = np.array([0.0, 0.0, 0.0])
#時間区間
t_min = 0
t_max = 10.0
#時間区間数
NT =  2000
skip =  1
#時間間隔
dt = (t_max - t_min) / NT
#座標点配列の生成
ts = np.linspace(t_min, t_max, NT + 1)

#衝突判定（戻値：法線ベクトル）
def CollisionDetection( r, R ):
	if(r[0] + radius > L/2):
		return 1, np.array([-1.0, 0.0, 0.0])
	if(r[0] - radius < -L/2):
		return 1, np.array([1.0, 0.0, 0.0])
	if(r[1] + radius > L/2):
		return 1, np.array([0.0, -1.0, 0.0])
	if(r[1] - radius < -L/2):
		return 1, np.array([0.0, 1.0, 0.0])
	if(r[2] + radius> R[2]):
		return 2, np.array([0.0, 0.0, -1.0])
	if(r[2] - radius < -L/2):
		return 1, np.array([0.0, 0.0, 1.0])
	#衝突なし
	return False, np.array([0.0, 0.0, 0.0])
	
######################################
#　４次のルンゲ・クッタ　クラス
######################################
class RK4:
	#コンストラクタ
	def __init__(self, m, dt = 0.01, r0 = np.array([0.0, 0.0, 0.0]), v0 = np.array([0.0, 0.0, 0.0]) ):
		self.m = m
		self.r = r0.copy()
		self.v = v0.copy()
		self.dr = np.array([0.0, 0.0, 0.0])
		self.dv = np.array([0.0, 0.0, 0.0])
		self.dt = dt
		self.fc = np.array([0.0, 0.0, 0.0])
		self.__v1 = np.array([0.0, 0.0, 0.0])
		self.__v2 = np.array([0.0, 0.0, 0.0])
		self.__v3 = np.array([0.0, 0.0, 0.0])
		self.__v4 = np.array([0.0, 0.0, 0.0])
		self.__a1 = np.array([0.0, 0.0, 0.0])
		self.__a2 = np.array([0.0, 0.0, 0.0])
		self.__a3 = np.array([0.0, 0.0, 0.0])
		self.__a4 = np.array([0.0, 0.0, 0.0])

	#速度を与えるメソッド
	def V(self, t, r, v):
		return v.copy()
	##########################################################
	#加速度を与えるメソッド
	def A(self, t, r_t, v_t):
		return F( t, self.m, r_t, v_t, self.fc ) / self.m
	##########################################################

	#時間発展を計算するメソッド
	def timeEvolution(self, t):
		#１段目
		self.__v1 = self.V(
			t, 
			self.r, 
			self.v
		)
		self.__a1 = self.A(
			t, 
			self.r, 
			self.v
		)
		#２段目
		self.__v2 = self.V(
			t + self.dt / 2.0, 
			self.r + self.__v1 * self.dt / 2.0, 
			self.v + self.__a1 * self.dt / 2.0
		)
		self.__a2 = self.A(
			t + self.dt / 2.0, 
			self.r + self.__v1 * self.dt / 2.0, 
			self.v + self.__a1 * self.dt / 2.0			
		)
		#３段目
		self.__v3 = self.V(
			t + self.dt / 2.0, 
			self.r + self.__v2 * self.dt / 2.0, 
			self.v + self.__a2 * self.dt / 2.0					
		)
		self.__a3 = self.A(
			t + self.dt / 2.0, 
			self.r + self.__v2 * self.dt / 2.0, 
			self.v + self.__a2 * self.dt / 2.0			
		)
		#４段目
		self.__v4 = self.V(
			t + self.dt, 
			self.r + self.__v3 * self.dt, 
			self.v + self.__a3 * self.dt					
		)
		self.__a4 = self.A(
			t + self.dt, 
			self.r + self.__v3 * self.dt, 
			self.v + self.__a3 * self.dt			
		)
		#差分の計算
		self.dr = self.dt * ( self.__v1 + 2.0 * self.__v2 + 2.0 * self.__v3 + self.__v4 ) / 6.0 
		self.dv = self.dt * ( self.__a1 + 2.0 * self.__a2 + 2.0 * self.__a3 + self.__a4 ) / 6.0 


def T(v, V):
	return 1.0/2.0 * m * np.sum( v**2 ) + 1.0/2.0 * M * np.sum( V**2 )
def U(r, R):
	return  - m * np.dot(g, r) - M * np.dot(g, R)
#################################################################
## 計算開始
#インスタンス
rk4_m = RK4(m, dt, r0, v0)
rk4_M = RK4(M, dt, R0, V0)

#アニメーション作成用
ims=[]
def plot( t, r, v, R, V ):
	E = T(v, V) + U(r, R)
	print( "{:.2f}".format(t), E )
	#各コマを描画
	img  = plt.plot([r[0]], [r[2]], colors[0], marker = 'o', markersize = 20, linestyle='solid', linewidth = 0 )
	img += plt.plot([-L/2, L/2], [R[2],R[2]], colors[1], linestyle='solid', linewidth = 10 )
	time = plt.text(0.8, 1.25, "t = " +  str("{:.2f}".format(t)) + r"$[{\rm s}]$" , ha='left', va='center', fontsize=18)
	#テキストをグラフに追加
	img.append( time )
	ims.append( img )

#各時刻における運動方程式の計算
for tn in range(len(ts)):
	t = ts[tn]
	if( tn%skip == 0 ): 
		plot(t, rk4_m.r, rk4_m.v, rk4_M.r, rk4_M.v )

	#ルンゲ・クッタ法による時間発展
	rk4_m.fc = np.array([0.0, 0.0, 0.0])
	rk4_M.fc = np.array([0.0, 0.0, 0.0])
	rk4_m.timeEvolution( t )
	rk4_M.timeEvolution( t )
	#位置と速度を更新
	rk4_m.r += rk4_m.dr
	rk4_m.v += rk4_m.dv
	rk4_M.r += rk4_M.dr
	rk4_M.v += rk4_M.dv

	CollisionFlag , n = CollisionDetection( rk4_m.r, rk4_M.r )
	if( CollisionFlag ):
		#いったん巻き戻して
		rk4_m.r -= rk4_m.dr
		rk4_m.v -= rk4_m.dv
		rk4_M.r -= rk4_M.dr
		rk4_M.v -= rk4_M.dv
		#衝突力
		if( CollisionFlag == 1 ):
			barV = - rk4_m.v
			rk4_m.fc = (2.0 * m * np.dot(barV, n) * n ) / dt - m * np.dot(g, n) * n
			rk4_M.fc = np.array([0.0, 0.0, 0.0])
		elif( CollisionFlag == 2 ):
			barV = rk4_M.v - rk4_m.v
			rk4_m.fc = 2.0 * m * M / ( m + M ) * np.dot(barV, n) * n / dt
			rk4_M.fc = - rk4_m.fc

		rk4_m.timeEvolution( t )
		rk4_M.timeEvolution( t )
		#位置と速度を更新
		rk4_m.r += rk4_m.dr
		rk4_m.v += rk4_m.dv
		rk4_M.r += rk4_M.dr
		rk4_M.v += rk4_M.dv

plt.title( u"シリンダーの中の粒子の運動", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$x\,[{\rm m}]$" , fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.ylabel(r"$z\,[{\rm m}]$" , fontsize=20, fontname="Yu Gothic", fontweight=1000)
#罫線の描画
plt.grid(which = "major", axis = "x", alpha = 0.8, linestyle = "-", linewidth = 1)
plt.grid(which = "major", axis = "y", alpha = 0.8, linestyle = "-", linewidth = 1)

#描画範囲を設定
plt.xlim([-1.2, 1.2])
plt.ylim([-1.2, 1.2])

#箱の概形
plt.vlines([-L/2], -L/2, L/2, "black", linestyles='solid', linewidth = 10 )
plt.vlines([L/2], -L/2, L/2, "black", linestyles='solid', linewidth = 10 )
plt.hlines([-L/2], -L/2, L/2, "black", linestyles='solid', linewidth = 10 )

#アニメーションの生成
ani = animation.ArtistAnimation(fig, ims, interval=10)
#ani.save("output.html", writer=animation.HTMLWriter())
#ani.save("output.gif", writer="imagemagick")
ani.save("output.mp4", writer="ffmpeg", dpi=300)

#グラフの表示
plt.show()

コメントを残す

メールアドレスが公開されることはありません。 ※ が付いている欄は必須項目です