pythonで物理シミュレーションしてみる・その1【卓球における空気抵抗の効果】

こんにちは、ヤスクマです！

今回はプログラミング言語のひとつであるpythonを使って物理現象をシミュレーションしていきたいと思います(自分の練習も兼ねて)。

物理を勉強すると最初に運動方程式を習うと思います。運動方程式を習うと鉛直投げ上げだの放物運動だのやりますが、実際にどんな運動をするのかというのをイメージするのは難しいのではないでしょうか？

そこで実際の物体の運動を計算しシミュレーションすることで理解していきたいと思います。

環境
運動方程式
微分方程式を解く
空気抵抗下の卓球ボール
おわりに

環境

anaconda環境でpython3でやっています。インストール方法が分からない人はgoogle colaboratoryでやるのがオススメです。

運動方程式

運動方程式は高校では ma=fで習うと思います。mは質量、aは加速度、fは力です。

この加速度aは速度vを微分したものです。さらに速度vは変位rの微分です。

実際には加速度などはベクトルは下のように書けます。

f:id:ys9m:20210214134929j:plain

この微分方程式を解くことによって物体の運動を予測することが出来ます。

微分方程式を解く

微分方程式を解けば運動を予測することが出来るのですが、実際に手の計算で解が求まることは、現実ではほとんどありません。

そこでプログラミングを用いて数値的に解いていきます。

今回はscipyのodeintを使って解いていきます。まず練習として下の簡単な微分方程式を解きます。

f:id:ys9m:20210214140213j:plain

まず必要なライブラリをインポートします。

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

次に微分方程式と初期値を定義します。

def func(y, t):
    dydt = -a*y
    return dydt
a = 1
y0 = 1 #初期値

次にodeintを使って微分方程式を解きます

t = np.arange(0, 6, 0.01) #tの刻み幅を0.01にして0から6までのリストを作成
y = odeint(func, y0, t) #関数,初期値,変数
plt.plot(t, y, label= 'exp')
plt.legend()
plt.show()

実行すると下のようなグラフが出てきます。このグラフは指数関数y = e^-tになっています。

f:id:ys9m:20210214140841j:plain

無事、微分方程式を解いてグラフ表示が出来ました。

空気抵抗下の卓球ボール

以上のプログラムを応用して、卓球における空気抵抗の効果を計算していきます。ボールには空気抵抗は速度の2乗に比例した空気抵抗がかかることが知られているそうです。 yの負の方向に重力がかかった状況で、x方向にボールを打つという状況を考えます。空気抵抗を考慮にいれた運動方程式は

f:id:ys9m:20210214141732j:plain

で書かれることになります。この空気抵抗ありの場合となしの場合でボールの軌道を計算していきます。

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from IPython.display import HTML

def func_1(state, t):   #空気抵抗あり
    dydt = np.zeros_like(state)
    dydt = [ state[1], -1 * k / m * state[1]2 , state[3], -G - k / m * state[3]2  ]　#[1]2, [3]2となっているところ本当は[1]2と[3]2
    return dydt

def func_2(state, t):  #空気抵抗なし
    dydt = np.zeros_like(state)
    dydt = [ state[1], 0 , state[3], -G   ]
    return dydt

#各種定数
G  = 9.8         # 重力加速度(m/s²)
m = 2.7 * 1e-3 # ボールの質量(kg)
k = m * G / 7.8**2    #空気抵抗(kg/m) 速度の2乗に比例する抵抗

#初期値
x0 = 0.0        # xの初期値
y0 = 0.2      # yの初期値
v0_km = 18.0      # 初速度 (km/h)
theta0 = 10.0      # 打ち上げ角度(deg)

#諸々の変換
v0 = v0_km * 1000 / 3600    # m/s に変換
v_x0 = v0 * np.cos( np.radians(theta0) ) # x初速度
v_y0 = v0 * np.sin( np.radians(theta0) ) # y初速度

# 初期状態
state =[ x0, v_x0 , y0, v_y0 ]

#微分方程式を解く
dt = 0.05
t = np.arange(0.0, 1, dt)

sol_1 = odeint(func_1, state, t)
sol_2 = odeint(func_2, state, t)

x_1 = sol_1[:,0]
y_1 = sol_1[:,2]
x_2 = sol_2[:,0]
y_2 = sol_2[:,2]

fig = plt.figure()
ax = fig.add_subplot(111, autoscale_on=False, xlim=(0, 1.525), ylim=(0, 0.4))
ax.set_aspect('equal')
ax.grid()

plt.plot(x_1, y_1, color = 'b', label= 'o resistance')
plt.plot(x_2, y_2, color = 'r', label= 'x resistance')
plt.legend()
plt.show()

空気抵抗の係数は https://pp-physics.com/technic/terminal-velocity/を参考にしました。

詳しい中身は少し難しいのですが、func_1で空気抵抗ありを、func_2で空気抵抗なしを計算しています。微分方程式を定義するところでリストを渡している以外は前章のものとほぼ同じです。卓球台の端の20cm上から時速18km/hで10度打ち上げてボールを打った時の軌道を計算しています。

実行して得られた軌道は↓のようになります。青が空気抵抗ありを、赤が空気抵抗なしを表しています。また横軸の幅は卓球台の幅になっています。

f:id:ys9m:20210214143049j:plain:w450

赤の空気抵抗なしではギリギリはいるボールが、空気抵抗を考慮に入れると余裕で入るようになっています。卓球において空気抵抗がいかに重要かというのが分かります。これに回転の効果が加わり実際の軌道が生まれることになります。今後回転の効果を入れて計算してみたいと思います。

ここでは一例を示しましたが、是非初速度や打ち上げ角度を変えて遊んでみてください。 Google Colaboratory ならインストールなしで実行することができてオススメです。 https://colab.research.google.com/notebooks/intro.ipynb?hl=ja