円軌道はむずかしい

万里鏡１号・弾道ミサイルの軌道（２）からの続き

北朝鮮の弾道ミサイルの簡単なシミュレーションコードをMathematicaで作っていた。これを少しアレンジすれば人工衛星を軌道に投入するところまでできそうな気がする。

早速，以前のコードを修正してみた。まずは通常の加速直後に角度方向だけに加速度を加えるようにしたがうまくいなかい。打ち上げ加速は投射角の方向になっているので，動径速度成分が大きく残っているうえ，加速すれば軌道は膨らむ。このため，離心率の大きな長円軌道になって地表にぶつかってしまうか，地球の重力圏から脱出してしまうのだ。

次に，打ち上げ加速の直後に空白時間をおいて，動径速度成分が小さくなったところで角度方向に加速できるようにした。それでもうまくいかない。簡単な試行錯誤では周回軌道にのせるのが難しい。そもそも角度方向に加速するということは面積速度すなわち角運動量をふやし，動径方向の微分方程式で軌道半径を膨らませる方向に作用してしまう。

そこで，後期加速では衛星をその速度ベクトルの方向に加速することにした。$t=0$で速度ベクトルをとりだすところに発散があったので，これを回避するため，地球の自転による2倍面積速度$h(t)$の初期値として，$h(0) = R^2\omega=R^2 \frac{2\pi}{24*3600}=2930$ km$^2/$sを与えた。初期加速度はこれまでの$\ a=0.0446$として$30$秒加速する。その後，800秒程度休止した後に，後期加速度$\ b=0.1445$（ここを微調整した）で$250$秒加速すると，なんとか軌道に投入することができた。投射角は$s=45$度，初期加速における燃料比は$p=0.85$であり，加速方向の角度には$0.3$をかけて動径成分を抑えた。

なかなか難易度の高いゲームである。衛星の軌道高度が1200km程度の準円軌道となっている。これを500kmにしなさいといわれても，こんな単純な2段階制御ではちょっと難しい。なお，プログラムの検証のため，$r=6850$kmの宇宙空間で第一宇宙速度に相当する$v=\sqrt{gr}=8.2$ km/sを角度方向の初速度として与えると，正確に円軌道を描くことが確かめられた。

g = 0.0098; R = 6350; τ = 30; τs = τ*27; τt = 250; p = 0.85;

 a = 0.0446; b = 0.1445 a; s = 45 Degree; T = 15400; 

fs[t_] := 0.3*ArcTan[r[t]*r'[t]/ h[t]]

fr[t_, τ_] :=  a*Sin[s]*HeavisideTheta[τ - t] + 

   b*Sin[fs[t]]*HeavisideTheta[t-τs-τ]

   *HeavisideTheta[τ+τs+τt-t]

ft[t_, τ_] :=  a*Cos[s]*r[t]*HeavisideTheta[τ - t] + 

   b*Cos[fs[t]]*r[t]*HeavisideTheta[t-τs-τ]

   * HeavisideTheta[τ+τs+τt-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] == 2930 + 0*Sqrt[g] R^(3/2)}, {r, h}, {t, 0, T}]

f[t_] := r[t] /. sol[[1, 1]]

d[t_] := h[t] /. sol[[1, 2]]

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

Plot[{f[t + 1] - f[t], d[t]*R/f[t]^2, d[t]/f[t]},

 {t, 0, T}, PlotRange -> {-5, 15}]

 

図：苦労すると有難みがわかる衛星の準円軌道のグラフ

P. S. もう少しがんばると，軌道高度650km（r=6980km）の準円軌道まで達成できた。

g = 0.0098; R = 6350; τ = 25; τs = τ*15.3; τt = 350;  p = 0.85;

a = 0.0446; b = 0.1275 a; s = 45 Degree; T = 15400;