弾道ミサイルの軌道（１）

Ｊアラートからの続き

このたびの弾道ミサイルが，大圏コースで4600kmを飛行したということは，角度にして40度強だ。また，高度1000kmというのは地球半径の15%にあたる。さすがに地表面を平面として一様重力場で考えるというのではちょっとマズイ気がする。

そこで，地球の中心を原点とする極座標系での運動方程式を考える。運動する物体の座標を$( r(t), \theta(t) )$とする。運動方程式は，$ m ( \ddot{r} - r \dot{\theta} ) = F_r, \quad \dfrac{m}{r} \dfrac{d}{dt} (r^2 \theta) = F_\theta $ となる。ここで面積速度の2倍を$h(t) = r^2 \dot{\theta}$と定義すると，$ m ( \ddot{r} - \dfrac{h^2}{r^3} ) = F_r, \quad \dfrac{m}{r} \dot{h} = F_\theta $ となる。

(1) 万有引力だけが働く場合：　$F_r = \dfrac{GMm}{r^2} = mg \dfrac{R^2}{r^2}, \quad F_\theta =0 $となる。与えられた条件は，到達距離 L=4600km，到達時間 T=1320 s，最高高度 H=1000 km，平均水平速度 vt =3.55 km である。また，地表重力加速度 g=9.8 m/s^2，地球半径 R=6350 kmとする。Mathematicaのコードで調整するパラメータは，鉛直方向の初速度だけである。到達距離・時間の条件を満たすものとして vr = 4.15  km/s が得られる。このときの速度はv0 = √(vt^2+vr^2) = 5.46 km/s（マッハ16）となる。しかし，最高高度が，1300 kmとなってうまく合わない。

g = 0.0098; R = 6350; vr = 4.15; vt = 3.55; h = R vt

sol = NDSolve[{r''[t] == h^2/r[t]^3 - g R^2/r[t]^2, 

      r[0] == R, r'[0] == vr}, r, {t, 0, 1320}]

f[t_] := r[t] /. sol[[1]]; f[660] - R

Plot[{6350, f[t]}, {t, 0, 1320}]

図１：初速度のみ与えたモデル（横軸: t ，縦軸: 軌道半径 r ）

 (2) 初期加速度が一定時間働く場合：　加速度の値（簡単のため動径方向と角度方向は等しいと仮定），加速時間の2つをパラメタとする。先ほどのように到達距離・時間の条件を満たすようにパラメタを探すと，加速度 0.025 km/s^2 （2.5G）と加速時間 180 s の値が得られた。このときの vt = d[τ]/f[τ] = 4.38 km/s，vr = f[τ+1] - f[τ] = 2.96 km/s，v0 = 5.29 km/s（マッハ15.5）となる。また，最高高度は1020 kmとなり，この場合は全体として辻褄があうことになる。

g = 0.0098; R = 6350; τ= 180;

fr[t_,τ_] := 0.025 * HeavisideTheta[τ- t]

ft[t_,τ_] := 0.025 * r[t] * HeavisideTheta[τ- t]

sol = NDSolve[{r''[t] == h[t]^2/r[t]^3 - g R^2/r[t]^2 + fr[t,τ], r[0] == R, r'[0] == 0, 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}]

図２：初速度0から一定の加速をする場合の軌道半径（横軸: t ，縦軸: 軌道半径 r ）

なお，加速終了時の射影水平速度は，vt = R d[180]/f[180]^2 = 4.16 km/s である。そこで180秒のあいだに進む水平距離は，vt τ/2 = 370kmとなる。残りは，1030km/4.16km/s = 248 s なので，計 428秒で1400 km（津軽海峡上空）に達する。

図の印象で騙されていたが，到達距離4600 kmの確認が済んでいなかった。解けているのは 角速度 $\dot{\theta} = \dfrac{h(t)}{r(t)^2}$なので，これを積分した $R \theta(T)$ が必要なのだ。下図より角速度の平均値が 0.00059だと仮定すると，370 km + 0.00059*R*1140 s = 370 + 4270 = 4640 kmとなる。1%の誤差でOKだった。

図３：初速度0から一定の加速をする場合の角速度（横軸: t ，縦軸: 角速度$\dot{\theta}$  ）

最高高度は，t=679s で1050 kmとなった。加速終了時の速度は，vt = d[τ] / f[τ] = 4.69 km/s，vr = f[τ+1] - f[τ] = 3.31 km/s，v0 = √(vt^2+vr^2) = 5.74 km/s （マッハ16.9）であり，ほぼ報道結果が再現された（P. S. と思ったが・・・）。

弾道ミサイルの軌道（２）に続く。