はじめに
以前ルンゲクッタ法を用いた数値解析の手法について説明しました。
今回は、自由落下する物体の壁面での跳ね返りをシミュレーションしたいと思います。
跳ね返りのシミュレーション手法はいくつか存在しますが、今回は簡単に実装できるペナルティ法を用いてシミュレーションしたいと思います。
自由落下の数値解析
まず、跳ね返りなしの場合の自由落下の数値解析を行ってみます。
高さのみの1自由度系で実装します。
ルンゲクッタ法の実装はこちらを参考に
matplotlibでのアニメーションの作成はこちらを参考にして作成しています。
実際に実装してみます
import matplotlib.pyplot as plt from matplotlib.patches import Circle from matplotlib.animation import FuncAnimation from matplotlib.animation import PillowWriter import numpy as np l = 2 # 直径 g = 9.8 # 重力加速度 v_0 = 0 # 初期の速度 x_0 = 10 # 初期の変位 x = x_0 v = v_0 dt = 0.02 # タイムステップ # グラフを作成 fig = plt.figure() ax = fig.add_subplot(1, 1, 1) patch = Circle(xy=(0, x_0), radius=l/2) im = ax.add_patch(patch) ax.set_xlim(-10, 10) ax.set_ylim(0, 15) ax.set_aspect('equal') # 運動方程式の関数化 def func_normal(v, x): '''運動方程式''' return -g def runge_kutta_method(x, v, func): '''ルンゲクッタ法''' k1_x = v * dt k1_v = func(v, x) * dt k2_x = (v + k1_v / 2) * dt k2_v = func(v + k1_v / 2, x + k1_x / 2) * dt k3_x = (v + k2_v / 2) * dt k3_v = func(v + k2_v / 2, x + k2_x / 2) * dt k4_x = (v + k3_v) * dt k4_v = func(v + k3_v / 2, x + k3_x) * dt x_new = x + (k1_x + 2 * k2_x + 2 * k3_x + k4_x) / 6 v_new = v + (k1_v + 2 * k2_v + 2 * k3_v + k4_v) / 6 return v_new, x_new # 初期化関数 def init(): # only required for blitting to give a clean slate. global x, v x = x_0 v = v_0 im.set_center((0, x)) return im, # グラフ更新関数 def update_anim(frame): global x, v v_new, x_new = runge_kutta_method(x, v, func_normal) v, x = v_new, x_new im.set_center((0, x)) return im, ani = FuncAnimation( fig, # Figureオブジェクト update_anim, # グラフ更新関数 init_func=init, frames = np.arange(0, 10, dt), # フレームを設定 interval = int(dt*1000), # 更新間隔(ms) repeat = True, # 描画を繰り返す blit = True, save_count = 500, ) # アニメーションのgif形式での保存 ani.save("drop.gif", writer='pillow') plt.show()
これでとりあえずは自由落下する物体のシミュレーションができました。
これの高さ0の位置でバウンドさせるように改造したいと思います
跳ね返りの実装
跳ね返りはペナルティ法を用いて表現します。
高さ0の位置に壁を作ることを考えると、ある時刻において壁にボールがめり込む瞬間が発生します。
その際にめり込み量に比例した反力をボールに加えてあげることで跳ね返りを表現する手法です。
ペナルティ法について
下図にボールがめり込んだ際の反力のイメージ図を示します。
めり込み量に応じた反力はバネで表現します。
また、反発係数の代わりに減衰を加えることで跳ね返り時の速度減衰を表現します。
ここでめり込んだ際のバネ係数を、減衰係数を
とすると運動方程式は下記のようにあらわすことができます。
めり込んでいないとき
めり込んでいるとき
ここでめり込み量はボールの半径を[tex l / 2]とした際、です。
めり込みを少なくするためにタイプステップは小さめに設定し、バネ係数は大きめに設定することが一般的だそうです。
今回はバネ定数を2000, 減衰係数を3, 質量を1, タイムステップを0.02としました。
バネ定数やタイムステップが大きすぎたりすると反力が強すぎて発散してしまいボールが消えてしまいますので注意して調整してあげる必要があります。
パラメータの調整次第で挙動が全然変わるので、厳密解とは異なりますが、複数のボールの跳ね返りなどをシミュレーションしたい際はこのペナルティ法が比較的簡単に、高速に動作するのでよく使われるみたいです。
ペナルティ法の実装
さっそく実装してみます。
めり込んだ際の運動方程式を別の関数で定義しておきます。
c = 3 # 跳ね返り時の減衰係数 k = 2000 # 跳ね返りの弾性係数 m = 1 # 質量 # めり込み時の運動方程式の関数化 def func_bound(v, x): '''運動方程式''' return - c / m * v - k / m * (x - l / 2) - g
これを使って、の位置によって自由落下の運動方程式を使うかめり込みの運動方程式を使うか場合分けをして位置を更新することで跳ね返りのシミュレーションが可能です。
プログラム全文は下記のとおりです。
import matplotlib.pyplot as plt from matplotlib.patches import Circle from matplotlib.animation import FuncAnimation from matplotlib.animation import PillowWriter import numpy as np l = 2 # 直径 g = 9.8 # 重力加速度 c = 3 # 跳ね返り時の減衰係数 k = 2000 # 跳ね返りの弾性係数 m = 1 # 質量 v_0 = 0 # 初期の速度 x_0 = 10 # 初期の変位 x = x_0 v = v_0 dt = 0.02 # タイムステップ fig = plt.figure() ax = fig.add_subplot(1, 1, 1) patch = Circle(xy=(0, x_0), radius=l/2) im = ax.add_patch(patch) ax.set_xlim(-10, 10) ax.set_ylim(0, 15) ax.set_aspect('equal') # 運動方程式の関数化 def func_normal(v, x): '''運動方程式''' return -g # めり込み時の運動方程式の関数化 def func_bound(v, x): '''運動方程式''' return - c / m * v - k / m * (x - l / 2) - g def runge_kutta_method(x, v, func): '''ルンゲクッタ法''' k1_x = v * dt k1_v = func(v, x) * dt k2_x = (v + k1_v / 2) * dt k2_v = func(v + k1_v / 2, x + k1_x / 2) * dt k3_x = (v + k2_v / 2) * dt k3_v = func(v + k2_v / 2, x + k2_x / 2) * dt k4_x = (v + k3_v) * dt k4_v = func(v + k3_v / 2, x + k3_x) * dt x_new = x + (k1_x + 2 * k2_x + 2 * k3_x + k4_x) / 6 v_new = v + (k1_v + 2 * k2_v + 2 * k3_v + k4_v) / 6 return v_new, x_new # 初期化関数 def init(): # only required for blitting to give a clean slate. global x, v x = x_0 v = v_0 im.set_center((0, x)) return im, # グラフ更新関数 def update_anim(frame): global x, v if x < l/2: v_new, x_new = runge_kutta_method(x, v, func_bound) else: v_new, x_new = runge_kutta_method(x, v, func_normal) v, x = v_new, x_new im.set_center((0, x)) return im, ani = FuncAnimation( fig, # Figureオブジェクト update_anim, # グラフ更新関数 init_func=init, frames = np.arange(0, 10, dt), # フレームを設定 interval = int(dt*1000), # 更新間隔(ms) repeat = True, # 描画を繰り返す blit = True, save_count = 500, ) # アニメーションのgif形式での保存 ani.save("bound.gif", writer='pillow') plt.show()
いい感じにバウンドが表現できました
