　　　　On a Thread of the Web: 弾道ミサイルの軌道（２）

弾道ミサイルの軌道（１）からの続き

今回のモデルによるシミュレーションで，軌道を再現する加速時間は180秒=3分だったが，これを少し変えるだけで，軌道は大きく変わってしまう。つまり，ミサイル発射直後の軌道推定は困難であり，少なくとも初期加速が終る時点までの数分間は待つ必要がある。と思ったのだが，コロラド先生は燃焼時間が1分だといっていた。

そこで，各方向の加速度を0.069 km/s^2 (7G)，加速時間を60 sにしたところ，vt = 4.07 km/s，vr = 3.58 km/s，v0 = 5.42 km/s（マッハ15.9）が得られた。このときの，最高高度が1160 km，到達距離が4270 km である。これでよいかと思ったが，調べてみると，そもそも推進剤の質量が全質量に対して非常に大きな割合を占めることがわかった。やり直し。

質量が変化する場合の極座標の運動方程式では，$m \ddot{\bm{r}}\ $の項に$\ \dot{m} \dot{\bm{r}}\ $が加わることになる。これを極座標にすれば，$\dot{m} (\dot{r} \bm{e}_r + r \dot{\theta} \bm{e}_\theta)\ $である。そこで，運動方程式の各成分は次の通り。

$\ddot{r} - \dfrac{h^2}{r^3} + \dfrac{\dot{m}}{m} \dot{r} = \frac{1}{m} F_r = -g \dfrac{R^2}{r^2} + \alpha H(\tau- t)$

$\dfrac{1}{r }(\dot{h} + \dfrac{\dot{m}}{m} h) = \frac{1}{m} F_\theta = 0 + \beta H(\tau-t) $

なお，$H(x) = 1 (x>0) ; =0 (x<0)$はヘヴィサイドの階段関数である。$\alpha, \beta$は，加速開始時刻 $t=0$から加速終了時刻 $t=\tau$ まで動径方向と角度方向に加わる加速度を表わす。

出発前の弾道ミサイルの全質量を$m_0$，全質量に対する推進剤の割合を $p$とすると，加速中（$0 \le t \le \tau$）の単位時間当たりの噴射質量は，$\dfrac{p m_0}{\tau}$となる。そこで，時刻 $t$における弾道ミサイルの質量は，$m(t) = m_0 - \dfrac{p m_0}{\tau} t = m_0 (1 - p\dfrac{t}{\tau})$，$\dot{m}(t) = - m_0 \dfrac{p}{\tau}$。

したがって，$\dfrac{\dot{m}}{m} = - \dfrac{p}{\tau - p t} \quad (0 \le t \le \tau)$ であり，$t>\tau$ では $\ \dfrac{\dot{m}}{m} = 0\ $となる。

なお，最高高度は，$\dot{r}(t)=0$となる時刻 $t_p$に対応する$r(t_p)$で与えられ，到達距離は，$\int_0^T \dfrac{R h(t)}{r(t)^2}dt$で得られる。ただし，$T$は到達時間である。

準備ができたので再計算してみる。全重量に対する推進剤の比率は，p=0.75（3/4）と仮定した。加速時間 60s，加速度 0.0446 km/s^2（4.6 G）で到達時間1320 sが再現される。最高高度は，t=682sのとき 1060 km。加速終了時の速度は，vt = 4.69 km/s，vr = 3.32 km/s，v0 = 5.75 km/s （マッハ16.9）である。また，到達距離は4930 kmであり，6-7%の誤差で報告値に一致した。

g = 0.0098; R = 6350; τ = 60; p = 0.75; a= 0.0446; s = Pi/4
fr[t_,τ_] := a * Sin[s] * HeavisideTheta[τ - t]
ft[t_,τ_] := a * Cos[s] * r[t] * HeavisideTheta[τ- t]
fm[t_,τ_] := -p / (τ- p * t) * HeavisideTheta[τ - t]
sol =
NDSolve[{r''[t] == -fm[t,τ] * r'[t] + h[t]^2/r[t]^3 -
g R^2/r[t]^2 + fr[t,τ], r[0] == R, r'[0] == 0,
h'[t] == -fm[t,τ] * h[t] + ft[t,τ], h[0] == 0},
{r, h}, {t, 0, 1320}]
f[t_] := r[t] /. sol[[1, 1]]
d[t_] := h[t] /. sol[[1, 2]]
Plot[{6350, f[t]}, {t, 0, 1320}]
Plot[{f[t + 1] - f[t], d[t]*R/f[t]^2, d[t]/f[t]}, {t, 0, 1319}, PlotRange -> {-4, 5}]
(f[t] - R) /. FindRoot[D[f[t], t] == 0, {t, 660}]
{d[τ]/f[τ], f[τ+1] - f[τ]}
{Sqrt[(d[τ]/f[τ])^2 + (f[τ+1] - f[τ])^2],
Sqrt[(d[τ]/f[τ])^2 + (f[τ+1] - f[τ])^2]/.34}
{NIntegrate[R d[t]/f[t]^2 , {t, 0, τ}],
NIntegrate[R d[t]/f[t]^2 , {t, 0, 1320}]}