Starship特集：打ち上げ能力計算ツールのご紹介

Contents

Rocket Performance Calculator ツールの紹介

本日はStarship特集第二回、Starshipのロケットシステムの考察と打ち上げ能力の推定に用いたツールについて使い方の説明をします。前回記事はこちらから：

この記事はこういう人向け

SpaceX特集第2回のツールの使い方を知りたい

液体ロケットのスペックが打ち上げ能力に及ぼす影響を知りたい

今回の記事のレーティング：

計算内容について少しだけ説明するわね。”pythonあるいは何かしらのプログラミング言語を少しだけいじったことがある”、”pipが使える”、これができれば大丈夫。難しい仕様はないので安心してね。

詳細な使い方はReadmeを参照してみてくださいね。

本ツールは一応MIT Licenseとさせて頂いておりますので、ご自由に利用、改変、配布頂いて構いません。(利用に対する責任や動作の保証は一切致しません、あくまでフリーユース、自己責任でお願い致します。)

動作環境・ツールの入れ方

ツールのDL

下記にてβ版を公開しています。

https://github.com/Gertrud-Violett/RocketSystem/tree/main/Performance_Calculator

Rocket System Calculation Tools
https://github.com/Gertrud-Violett/RocketSystem
1 forks.
1 stars.
0 open issues.

Recent commits:

githubのアカウント持っていない方は、”Performance_Calculator”フォルダごとダウンロードしてしまいましょう。必要なファイルはすべてここに入っているわ。

使用ライブラリ (Dependencies)

動作のためには、事前に下記の設定が必要です：

python 3.10以上 (mac, linux, windows 何でもokです)
ライブラリ：numpy, matplotlib, toml, sympy, pymap3d, plotly, scipy, pandas
rocketCEA：下記に導入方法を書いてあります。

rocketCEAはエンジンのIspや推力を計算するのに使用するプログラムです。下記の記事で導入方法を書いてあるので、試してみてください！

rocketCEAの入れ方、使い方が分からない。。。記事の通りやってもうまくいかない場合は、コメントやツイッターでご連絡ください。可能な限りサポートします！

使い方

前提

このツールはあくまで地表からロケットが離床し、軌道に到達するまでの計算を考えています。そのため、下記のような制約があります。

地球中心を原点とし、ECI座標を採用しています。
地球の公転運動は一切考慮していません。
極座標による計算のため、地球からの距離が離れると角度履歴の誤差が蓄積してしまう欠点があります。
3DoF計算を行っています。エンジンのジンバル角やロケットの迎角は一切考えていません。ロケットは推力の向きのみ自在に変えられる質点です。

打ち上げ能力推定用のプログラムであって、軌道設計用のプログラムではない、という点にご注意ください。

入力ファイル

3つのインプットファイルに諸元を記入します。

settings.toml

toml形式のテキストファイルです。機体やエンジンの各種諸元を定義しています。サンプルファイルを参照ください。

名称	内容
[ENGINE]	エンジンセクション
[ENGINE.Stage1]	1段目セクション
Throat_Dia	スロート径
Isp_efficiency	Isp効率 (全平衡流)
Combustion_Pressure	燃焼室圧
Nozzle_ExpansionRatio	ノズル膨張比
OF_mixture	O/F 酸燃比
Target_thrust	目的推力 (特になくても動きます)
Engine_qty	エンジン基数
[VEHICLE]	機体セクション
[VEHICLE.Stage1]	1段目セクション
Structural_Efficiency	構造効率
Initial_Weight	初期重量
Fuel_Remaining	着陸用の残渣推進剤の割合
Fuselage_Dia	機体直径、空力の計算に使います
[VEHICLE.Payload]	Payloadセクションは推力=0になります。
[ENVIRONMENT]	環境セクション
gravity	重力定数
barometric_pressure	大気圧 (MPa、初期は0.1013MPa)
[VECTOR]	初期ベクトルセクション
Latitude	初期緯度
Longitude	初期経度
altitude	初期高度
radius	地球の半径 (m)
[CALC_SETTINGS]	計算に用いる設定セクション
Casename	スタディの名称(出力ファイルに記載)
timestep	Δt時間刻み。推奨は0.1~0.01s
timelimit	計算をストップする時刻。指定しないと永遠に動いてしまいます。
betafile	attitude_beta.csv設定ファイルの名称
gammafile	attitude_gamma.csv設定ファイルの名称

attitude_beta.csv

東西(φ)方向の姿勢角履歴を記したcsvファイルです。時刻tにおいて地面に対して何度の推力を発生させるか入力します。注意点としては、機体から見た相対角ではなく、地表に対する絶対角で定義が必要です。

attitude_gamma.csv

南北(θ)方向の姿勢角履歴を記したcsvファイルです。時刻tにおいて地面に対して何度の推力を発生させるか入力します。

実行ファイル

FlightCalcPolar.py

メインの計算です。Argumentは設定ファイル名のみです。計算を実行すると、結果が./plotに保存されます。計算中は時刻や重量、加速度等の中間データがターミナル上にprintされます。この例のように、timelimitで設定した上限時間の諸元マトリクスが出力されたら無事計算が完了しています。

python FlightCalcPolar.py "settings.toml"
>>
Stage:3,Time[s]:9872.2,Weight[t]:100, Drag[kN]:0.00,Thrust[kN]:0,Alt[m]:43642305, Accel[m/s2]:0.26,Velocity[m/s]:4413
 Accel: -0.12029888702628849 -1.1586735571990027e-09 -4.4850806980064825e-09
 Stage:3,Time[s]:9999.9,Weight[t]:100, Drag[kN]:0.00,Thrust[kN]:0,Alt[m]:44176570, Accel[m/s2]:0.25,Velocity[m/s]:4394
 Downrange[km]:63593
 Weight[t]:100
 Drag[kN]:0.00
 Alt[km]:44176.57
 Orbit Speed [m/s]1338.59
 Vehicle ABS speed [m/s]:4394
 Vehicle Airspeed [Mach]:4.330
 1st Time[s],alt[km]: 182,98
 2nd Time[s],alt[km]: 434,261
 Payload Time[s],alt[km]: 10000,44177
 Theta[deg]:64.0
 DeltaV_th[m/s]:9584.93
 DeltaV1st[m/s]:3952.83
  Landing DeltaV[m/s]:3355.02
        Time[s]     ECI-X[km]     ECI-Y[km]     ECI-Z[km]  ...    drag[N]    Weight[kg]   Altitude[m]  Velocity[m/s]
0           0.0   -712.913495  -5682.003585   2792.173098  ...   0.000000  5.300000e+06  0.000000e+00     330.702310
1           0.0   -712.884016  -5682.007406   2792.173157  ...   0.000000  5.297984e+06  0.000000e+00     297.250720
...         ...           ...           ...           ...  ...        ...           ...           ...            ...
99999    9999.8  14520.452513  44907.931139 -18096.711690  ...   0.000000  1.000000e+05  4.417615e+07    4393.821190
100000   9999.9  14520.442117  44908.334057 -18096.886609  ...   0.000000  1.000000e+05  4.417657e+07    4393.806369
[100001 rows x 14 columns]
#ここで計算終了、csvファイルに結果が出力されます。

出力ファイル

計算を実行すると./plotに下記ファイルが出力されます。

result_Study N_XX E_YY_Casename_StgN.txt: N段目計算終了時の状態を出力したtxtファイル。XX、YYは緯度経度、StgNはN段目を意味します。段数分だけファイルが吐き出されます。
Study N_XX E_XX_Casename_plot.jpg: 各種計算結果のプロットです。速度や高度等の情報が出力されます。
Study N_XX E_YY_Casename_Trajectory2D.jpg：北極から見た軌道のプロット。緑線1段目、オレンジ線2段目、青線3段目です。
Study N_XX E_YY_Casename_Trajectory3D_ZZ.jpg：角度ZZから見た軌道の3Dプロット
out_XYZ_Study N_XX E_YY_Casename.csv：計算結果を時系列に保存したcsvファイル。
out_rocketarrStudy N_XX E_YY_Casename.csv：上記csvファイルのうちエンジン関係の結果のみを出力した簡易版

まずはダウンロードしたファイルを編集せずに実行して、下記のようなグラフが出力されるか確認してみてね。
動くことが確認できれば、あとは自由にいじってヨシ！

プロットツール

RunTrajectoryPlotter_polar.py

計算結果のcsvファイルをplotlyでhtml形式の三次元インタラクティブプロットに変換するものです。FlightCalcPolar.py実行後にお使いください。Argumentとして、読み込みしたいcsvファイルを入力します。

python RunTrajectoryPlotter_polar.py "out_XYZ_Study N_90.0 E_0.0_Demo_run.csv"
>>plotly3d_out_XYZ_Study N_90.0 E_0.0_Demo_run.csv.html

元のツールTrajectoryPlotterはbluedackさん作成のものをお借りしています：

https://github.com/bluedack-space/TrajectoryPlotter
1 forks.
1 stars.
0 open issues.

Recent commits:

Add files via uploadClass routines are added in TrajectoryPlotter, GitHub
Add files via upload, GitHub
Add files via upload, GitHub

その他のファイル

空力諸元以外は特に編集する必要はありませんが、それぞれの役割は下記の通りです。

earth.xlsx: TrajectoryPlotter用地球描画用ファイル。いじらないでください。
PandasHandler.py: TrajectoryPlotter用csv読み込みスクリプト
rocketthrustcalc.py: rocketCEAからエンジン諸元を拾うためのスクリプト。
USSTDATM: US Standard Atmosphere 大気情報と空力諸元を設定したファイル。空力諸元はM-Vをベースに迎角を平均化して3倍の抗力係数としていますが、修正したい場合はMVspecとaerodragメソッドをエディットしてください。
vehicle.py: 各時系列で力を計算して、次のタイムステップの出力をする計算のメイン部分

計算内容の説明

教科書的な四則演算で解いているから、コードを見れば分かるようになっているわ。pythonには行列を一発計算できる便利なライブラリもあるけれど、作者はプログラミングよわよわだから使っていないのよね。

また余計な一言を。。。テクラさんの言う通り、fortranをそのまま書き写しただけなようなコードになっています汗

極座標運動方程式について

vehicle.pyでは極座標運動方程式を解いています。各方向にかかる力としてはaccel_matメソッドにて、下記を外力として、加速度を計算しています。

r: 重力、推力、抗力
θ：betaの角度に応じた推力、抗力
φ：gammaの角度に応じた推力、抗力

差分式

ツールでは4次精度の線形近似で差分を計算しています。

古典ルンゲ・クッタ法ね。RGK_Matrixが差分計算のメソッドになるわ。

下記のまま、直接打ち込んでいます。。。

実装は下記のサイトで解説されています。

https://bombrary.github.io/blog/posts/eular_runge-kutta/

全体計算

メインの計算部であえるvehcalcメソッドでは、各時系列において与えられた状態量を運動方程式に代入、加速度を求め、差分法から次の時刻の状態を予測、これの逐次計算をしながら各時刻の結果をマトリクスに格納して行っています。

おおまかな流れとしてはこんな感じだわ。教科書的な時系列数値計算だと思うけれど、こういう使い方ができると面白いわよね。

パラメトリックスタディについて

Starshipの能力推定で行ったようなパラメータを振った感度解析を行う用に、FlightCalc_study.py全体をメソッド化して、別途感度解析用のファイルparastudy.pyを用意しました。それぞれご参考として、Parametricフォルダに入れてあるので、興味のある方は試してみてください。

ただしこちらはツール化できていないため、コードを直接編集する必要があります！ご自由に改変ください。

Q&A

計算トラブル

計算が動かない

必要な環境がすべて揃っているか、特にrocketCEAが正常動作しているか確認してみてください。

出力はされるが、エラーメッセージが連発する

下記のようなエラーの場合、rocketCEAが膨張比が計算できないと警告するエラーです。しかし計算には影響ないため、無視してしまって問題ありません。

=================  WARNING  ==========================
 Bad Solution in CalcPCoPE
 (desired Area Ratio= 5.65904043580018e+19   actual= 87798352321524.33 )
 input gam= 1.115665883934368  eps= 5.65904043580018e+19
=================  WARNING  ==========================

計算が発散する。収束しない

加速度に振幅が生じてしまう、質量が0になったり、諸元の設定がまずいと計算が発散することがあります。timestepを増やし、timelimitの値を小さくして、発散直前の挙動を確認してみてください。また初期のtimeとベクトルはコードの中で変更できるので、直前の値を代入して動作を確認することもできます。

計算が重い、速度が遅い。

どうしてもpythonの仕様上、timestepが0.01未満の状態で10000s以上計算すると数十分かかることもあります。rocketCEAのアクセス速度とnparrayのappend速度に起因すると思われます。0.1s程度でも比較的誤差が少ないと思うので、timestepを荒くしてみてください。また、地球を1周以上するような長時間のフライトはちゃんとした軌道計算ソフトをご利用ください。(NASAなどから無償で公開されています。)

おかしな軌道に投入されてしまう

attitude_beta.csv, attitude_gamma.csvがおかしな値になっていないか確認してみてください。また機体諸元が正常かどうか、初期位置を赤道沿い、attitude_gammaをすべてゼロにした状態で正常に動くか確認をしてみてください。赤道で軌道投入が成功しない場合は、他の緯度では困難な可能性が高いです。