Pythonで学ぶ常微分方程式の数値解（応用編：ロケットの運動）

物理学の具体的な計算にPythonを使って、物理学もPythonも同時に学んでしまいましょう。今回はPythonを使って、ロケットの運動を記述する常微分方程式の解を数値的に求めてみたいと思います。本記事で記載している図の描画の確認は、JupyterNotebook上で行っています。インターフェイスなどの動作環境の違いによって適宜変更点があるかもしれません。

今回の記事は以下の記事の発展的な内容になっているので、あわせて読んでいただければさいわいです：

pianofisica.hatenablog.com

運動方程式
モデルロケット
数値計算による運動の解析（空気抵抗を無視した場合）
数値計算による運動の解析（空気抵抗を加えた場合）

運動方程式

時刻 $t=0$ で発射する質量 $m(t=0)=M$ のロケットの運動

$\quad \displaystyle{m(t)\frac{d^2x}{dt^2}=f(t)-m(t)g}$

を考えましょう。ここで、燃料の噴射によってロケットは推進力 $f(t)$ を受け、それと同時にロケットは質量を失っていき、その質量は時刻 $t$ の関数 $m(t)$ になります。簡単のために推進力は時間的に一定であるとし、燃料の噴射量も時間的に一定であると仮定します。数式で表せば

$\quad \displaystyle{\begin{aligned}m(t)&=\begin{cases}\ M- m_{\rm fuel}\,\frac{t}{t_0} &\qquad(0\leq t\leq t_0)\\ \ M-m_{\rm fuel} &\qquad (t>t_0) \end{cases}\\ f(t)&=\begin{cases}\ f \qquad (0\leq t\leq t_0) \\ \ 0 \qquad (t>t_0)\end{cases}\end{aligned}}$

です。ここで $t_0$ は燃料の噴射時間です。燃料の噴射時間が短いほど、単位時間あたりの噴射量は大きくなってより大きな推進力が得られ、逆に噴射時間が短いと推進力は小さくなります。つまり、燃料から得られるトータルの推力は燃料の種類・量から決まるものですが、それを推進力として利用するときには、それをどれくらいの時間・ペース配分で使うかという要素がエンジンの出力特性を決める問題に関わってきます。ですがここでは簡単のために、推進力は噴射時間に反比例し、時間的に一定のペースで推進力を出し続けると仮定します：

$\quad \displaystyle{f=\frac{\varphi}{t_0}}$

ここで $\varphi [{\rm N}\cdot{\rm s}]$ はトータルの推力です。

モデルロケット

いま、より現実の問題と近づけるために『モデルロケット』の具体的なデータを参考にして運動を解析してみましょう。

www.ja-r.net

全質量 $M$ は 50グラム、燃料 $m_{\rm fuel}$ は 8グラム、合計推進力 $\varphi$ は 3 $[\rm{N}\cdot\rm{s}]$ 、噴射時間 $t_0$ は 3.5秒、重力加速度 $g$ は 9.8 $[ \rm{m}/\rm{s}^2]$ としましょう。

数値計算による運動の解析（空気抵抗を無視した場合）

過去記事

pianofisica.hatenablog.com

で解説した常微分方程式の数値解法を今回の問題設定にあてはめてみます：

h = 0.01   #  時間間隔（s）
N = 1100   #  時間間隔 h を単位にとったときの運動を追跡する時間

M = 50*10**(-3)   # ロケットの質量（kg）
mf = 8*10**(-3)   # 燃料の質量（kg）
phi = 3    # 積算推進力（Ns）
t0 = 350   # 時間間隔 h を単位にとったときの燃焼時間
g = 9.8   # 重力加速度(m/s^2)
fc = phi/(t0*h)   # 推進力

T = [0]  # 時刻
Mt = []   # 各時刻での質量
ft = []   # 各時刻での推進力
def t(k): return h+T[k-1]
def m(k): 
     if k <= t0: return M - mf*k/t0
     else: return M - mf
def f(k): 
     if k <= t0: return fc
     else: return 0
for k in range(N):
      T.append(t(k+1))
      Mt.append(m(k))
      ft.append(f(k))

x_pos = [0]  # 各時刻でのロケットの位置
p_vel = [0]   # 各時刻でのロケットの速度
def x(k): return x_pos[k-1]+h*(p_vel[k-1])+h**2/2*(ft[k-1]/Mt[k-1]-g)
def p(k): return p_vel[k-1]+h*(ft[k-1]/Mt[k-1]-g)
for k in range(N):
      x_pos.append(x(k+1))
      p_vel.append(p(k+1))
      
import matplotlib.pyplot as plt
plt.plot(T,x_pos)
plt.show()
print('最高到達点(m)=', max(x_pos), '最高速度(km/h)=', max(p_vel)*60**2/1000)

によって得られる時間（横軸）と高度（縦軸）の関係を図示したのが次のグラフです：

f:id:pianofisica:20200928122144p:plain

滞空時間 11秒程度、最高高度 101メートル、最高速度は時速 112キロと計算されます。

数値計算による運動の解析（空気抵抗を加えた場合）

いま、運動方程式の右辺に、推進力、重力のほかに空気抵抗

$\quad \displaystyle{D=\frac{1}{2} \rho v^{2} S C_{\mathrm{D}}}$

を加えましょう。ここで $\rho$ は流体密度で 1.2 $[\rm{kg}/\rm{m}^3]$ 、 $v$ は物体の速度、 $S$ は物体の断面積、 $C_{\mathrm{D}}$ は効力係数で円柱の場合 1.2とされます。いまモデルロケットを半径 2 $[\rm{cm}]$ の円柱だとみなすことにしましょう。以下のような数値計算

N = 1100
h = 0.01

M = 50*10**(-3)
mf = 8*10**(-3)
phi = 3
t0 = 350
g = 9.8
fc = phi/(t0*h)
rho = 1.2
Cd = 1.2
d = 1/2*rho*Cd*(2*10**(-2))**2*3.14

T = [0]
Mt = []
ft = []
def t(k): return h+T[k-1]
def m(k): 
     if k <= t0: return M - mf*k/t0
     else: return M - mf
def f(k): 
     if k <= t0: return fc
     else: return 0
for k in range(N):
      T.append(t(k+1))
      Mt.append(m(k))
      ft.append(f(k))

x_pos = [0]
p_vel = [0]
def x(k): return x_pos[k-1]+h*(p_vel[k-1])+h**2/2*(ft[k-1]/Mt[k-1]-g-d*p_vel[k-1]**2)
def p(k): return p_vel[k-1]+h*(ft[k-1]/Mt[k-1]-g-d*p_vel[k-1]**2)+h**2*d*p_vel[k-1]
for k in range(N):
      x_pos.append(x(k+1))
      p_vel.append(p(k+1))
      
import matplotlib.pyplot as plt
plt.plot(T,x_pos)
plt.show()
print('最高到達点(m)=', max(x_pos), '最高速度(km/h)=', max(p_vel)*60**2/1000)