/rocket_sizing.py
Last active Jul 28, 2019
多段ロケットの最適質量配分（サイジング）問題の計算
| # -*- coding: utf-8 -*- | |
| # ====== | |
| # 多段ロケットの最適質量配分（サイジング）問題の計算 | |
| # 必要な軌道速度に空力損失、重力損失、推力損失、制御損失を追加し、 | |
| # トータルの⊿Vを事前に算出し、その軌道速度に必要なサイジングを行う。 | |
| # 初期検討段階にのみ使用可能。 | |
| # | |
| # 入力: | |
| # 各段のIsp[秒] | |
| # 各段の構造比(0.0~1.0)（各段の全備重量と推進剤以外の割合） | |
| # 推力[N](オプション) | |
| # 消費推進剤割合（0.0~1.0）（オプション）(搭載推進剤に対する消費推進剤の割合) | |
| # 投棄物重量[kg]（オプション） | |
| # エンジン基数[基]（オプション） | |
| # 海面上推力の値を出力・使用するかどうか(bool)（オプション） | |
| # エンジン1基あたりのノズル出口面積[m2]（オプション） | |
| # | |
| # if __name__ == '__main__':後のUSER INPUT部分を編集すると | |
| # パラメータの変更が可能 | |
| # | |
| # cf. 半揚 稔雄(2014) 「ミッション解析と軌道設計の基礎」 | |
| # 8.5.1 ロケットの飛翔運動 ― 最適質量配分問題 p.193 | |
| # | |
| # Copyright (c) 2016-2017 Takahiro Inagawa | |
| # This code is released under the MIT License. | |
| # http://opensource.org/licenses/mit-license.php | |
| # ====== | |
| from __future__ import unicode_literals | |
| from __future__ import print_function | |
| from __future__ import division | |
| import sys | |
| import numpy as np | |
| from scipy import optimize | |
| import pandas as pd | |
| g0 = 9.80665 | |
| class Rocket: | |
| def __init__(self, Isp, stracture_ratio, | |
| thrust=0, propellant_consumption_rate=1.0, | |
| jettison=0, number_of_engine=1, | |
| use_SeaLevel = False, nozzle_exit_area=0): | |
| """ロケットクラス | |
| Args: | |
| Isp (float) : 真空中比推力[sec] | |
| stracture_ratio (float) : 構造係数、消費推進剤以外の割合 (0.0~1.0) | |
| thrust (float, optional) : 推力 [N] | |
| propellant_consumption_rate (float, optional) : 推進剤消費割合(0.0~1.0) | |
| jettison (float, optional) : 投棄物質量 [kg] | |
| number_of_engine (int, optional) : エンジン基数 | |
| use_SeaLevel (bool, optional) : 海面上推力を考慮するかどうか (default : False) | |
| nozzle_exit_area (float, optional) : エンジン1基あたりのノズル出口面積 [m2/基] | |
| """ | |
| self.Isp = Isp | |
| self.s = stracture_ratio | |
| self.mass = 0 | |
| self.mass_ratio = 0 | |
| self.upper_mass = 0 | |
| self.thrust = thrust | |
| self.burn_time = 0 | |
| self.acc_liftoff = 0 | |
| self.acc_cutoff = 0 | |
| self.pc_rate = propellant_consumption_rate | |
| self.residual_propellant = 0 | |
| self.jettison = jettison | |
| self.num_engine = number_of_engine | |
| self.use_SeaLevel = use_SeaLevel | |
| self.nozzle_exit_area = nozzle_exit_area | |
| self.payload = 0 | |
| self.thrust_SL = 0 | |
| def mass_ratio_calc(self, rambda): | |
| temp = rambda * g0 * self.Isp | |
| self.mass_ratio = (1 + temp) / (temp * self.s) | |
| def mass_calc(self): | |
| self.mass = (self.mass_ratio - 1) / (1 - self.s * self.mass_ratio) * \ | |
| self.upper_mass | |
| self.mass_stracture = self.mass * self.s | |
| self.mass_prop = self.mass - self.mass_stracture | |
| self.mass_prop_gross = self.mass_prop / self.pc_rate # 搭載推進剤 | |
| self.mass_prop_residual = self.mass_prop_gross - self.mass_prop # 残渣推進剤 | |
| self.mass_stracture_net = self.mass_stracture - self.mass_prop_residual - self.jettison # 構造重量 | |
| def deltaV_calc(self): | |
| self.m0mf = (self.mass + self.upper_mass) / (self.mass_stracture + self.upper_mass) | |
| self.deltaV = self.Isp * g0 * np.log(self.m0mf) | |
| def burn_time_calc(self): | |
| if (self.use_SeaLevel): | |
| self.thrust_SL = self.thrust - self.nozzle_exit_area * self.num_engine * 101300 | |
| self.burn_time = self.mass_prop * self.Isp * g0 / self.thrust | |
| self.acc_liftoff = self.thrust / g0 / (self.upper_mass + self.mass) | |
| self.acc_cutoff = self.thrust / g0 / (self.upper_mass + self.mass_stracture) | |
| def sizing_Lagrange_multiplier(x, velocity, rocket_array): | |
| temp = 1 | |
| Isp_std = rocket_array[0].Isp | |
| for r in rocket_array: | |
| temp_Isp = r.Isp / Isp_std | |
| temp *= ((1.0 / (x * g0 * r.Isp * r.s)) + (1.0/r.s)) ** (temp_Isp) | |
| temp -= np.exp(velocity / g0 / Isp_std) | |
| return temp | |
| if __name__ == '__main__': | |
| # print("==== Optimal Rocket Sizing for Multi-Stage Rocket ====") | |
| print("==== 多段ロケットの最適質量配分問題 ====") | |
| # ==== USER INPUT ==== | |
| rocket_name = "test rocket" | |
| payload = 140 # [kg] | |
| velocity = 9.6 # [km/s] | |
| # stage = Rocket(Isp[s], stracture_ratio[-], (optinal:thrust[N] | |
| # propellant_consumption_rate [-], | |
| # jettison [kg], number_of_engine [-], | |
| # use_SeaLevel(bool), nozzle_exit_area [m2])): | |
| stage1 = Rocket(290.0, 0.11, 288000, 0.99, 50, 9, True, 0.05) | |
| stage2 = Rocket(310.0, 0.16, 32000, 0.98, 0, 1, False, 0.1) | |
| r_array = [stage1, stage2] | |
| # ==== USER INPUT END ==== | |
| total_mass = payload | |
| total_mass_stracture = 0 | |
| total_mass_prop = 0 | |
| velocity *= 1000 # km/s -> m/s | |
| min_Isp = np.inf | |
| for r in r_array: | |
| min_Isp = min(min_Isp, r.Isp) | |
| limit = -1.0 / g0 / min_Isp # base of exponential must not be negative | |
| try: | |
| sol = optimize.brentq(sizing_Lagrange_multiplier, | |
| -100, limit, args = (velocity, r_array)) # solver | |
| except ValueError: | |
| print(u"*** NO SOLUTION ***") | |
| print(u"Please modification Isp and structure ratio") | |
| sys.exit() | |
| temp = payload | |
| for r in r_array: | |
| r.mass_ratio_calc(sol) | |
| for r in reversed(r_array): | |
| r.upper_mass = temp | |
| r.mass_calc() | |
| temp = r.upper_mass + r.mass | |
| r.deltaV_calc() | |
| if(r.thrust != 0): | |
| r.burn_time_calc() | |
| stage = 0 | |
| data_col = [("項目", "単位"), | |
| ("質量m0", "kg"), ("質量mf", "kg"), | |
| ("Isp(vac)", "秒"), ("構造係数", "-"), | |
| ("質量比", "-"), ("各段質量m0", "kg"), | |
| ("各段質量mf", "kg"), ("構造重量", "kg"), | |
| ("残渣推進剤", "kg"), ("投棄物","kg"), | |
| ("ペイロード", "kg"), | |
| ("正味の構造効率(構造重量/各段m0)","-"), | |
| ("消費推進剤重量", "kg"), ("搭載推進剤重量", "kg"), ("推進剤消費率", "%"), | |
| ("delta V", "m/s"), | |
| ("推力(vac)", "N"), ("燃焼時間", "秒"), | |
| ("エンジン数", "基"), ("出口面積", "m2/基"), | |
| ("推力(S.L.)", "N"), | |
| ("加速度_ignition", "G"), ("加速度_cutoff", "G")] | |
| data_col_multi = pd.MultiIndex.from_tuples(data_col, names=['項目', '単位']) | |
| df = pd.DataFrame(columns=data_col_multi) | |
| for r in r_array: | |
| stage += 1 | |
| print("%d: 質量m0 = %.1f [kg]" % (stage, r.upper_mass + r.mass)) | |
| print(" 質量mf = %.1f [kg]" % (r.upper_mass + r.mass_stracture)) | |
| print(" 各段質量m0 = %.1f [kg]\t構造比 = %.3f" % (r.mass, r.s)) | |
| print(" 各段質量mf = %.1f [kg]\t質量比 = %.3f" % (r.mass_stracture, r.mass_ratio)) | |
| print(" Isp(vac) = %d [s]\t\t消費推進剤質量 = %.1f [kg]" % (r.Isp, r.mass_prop)) | |
| print(" delta_V = %d [m/s]" % (r.deltaV)) | |
| print(" 構造重量 = %.1f [kg]" % (r.mass_stracture_net)) | |
| if(r.pc_rate != 1.0): | |
| print(" 残渣推進剤 = %.1f [kg]\t推進剤消費率 = %.1f [%%]" % ( r.mass_prop_residual, r.pc_rate * 100)) | |
| if(r.jettison != 0.0): | |
| print(" 投棄物 = %.1f [kg]" % (r.jettison)) | |
| print(" 正味の構造効率(構造重量/各段m0) = %.3f" % (1 - r.mass_stracture_net / r.mass)) | |
| if(r.thrust != 0): | |
| print(" ----") | |
| print(" 推力(vac) = %d [N]\t燃焼時間 = %d [s]" % (r.thrust, r.burn_time)) | |
| print(" 推力(S.L.) = %d [N]" % (r.thrust_SL)) | |
| print(" エンジン基数 = %d [基]\t出口面積 = %.3f [m2/基]" %(r.num_engine, r.nozzle_exit_area)) | |
| print(" 加速度@点火 = %.2f [G]\t加速度@CutOff = %.2f [G]" % (r.acc_liftoff, r.acc_cutoff)) | |
| print(" ====") | |
| total_mass += r.mass | |
| total_mass_stracture += r.mass_stracture | |
| total_mass_prop += r.mass_prop | |
| # ↓ for output csv file | |
| output_list = [str(stage)+"段", | |
| round(r.upper_mass + r.mass, 1), | |
| round(r.upper_mass + r.mass_stracture, 1), | |
| r.Isp, r.s, round(r.mass_ratio, 2), round(r.mass, 1), | |
| round(r.mass_stracture, 1) , round(r.mass_stracture_net, 1), | |
| round(r.mass_prop_residual, 1), round(r.jettison, 1), | |
| round(r.payload), | |
| round(1 - r.mass_stracture_net / r.mass, 3), | |
| round(r.mass_prop, 1), round(r.mass_prop_gross, 1), | |
| round(r.pc_rate * 100, 1), | |
| round(r.deltaV , 1), | |
| round(r.thrust), round(r.burn_time, 1), | |
| round(r.num_engine), round(r.nozzle_exit_area, 3), | |
| round(r.thrust_SL), | |
| round(r.acc_liftoff, 2), round(r.acc_cutoff, 2)] | |
| df_temp = pd.DataFrame([output_list], columns=data_col) | |
| df = df.append(df_temp) | |
| print("ペイロード: %.1f [kg]" % (payload)) | |
| print("=====") | |
| print("Total: 全備質量 = %d [kg]" % (total_mass)) | |
| print(" 構造質量 = %d [kg]" % (total_mass_stracture)) | |
| print(" 消費推進剤 = %d [kg]" % (total_mass_prop)) | |
| print(" delta_V = %.2f [km/s]" % (velocity/1000)) | |
| output_list = ["合計", round(total_mass, 1), "", "", "", "", "", "", "", "", "", | |
| round(payload), "", "", "", "", round(velocity), "", "", "", "", | |
| "", "", ""] | |
| df_temp = pd.DataFrame([output_list], columns=data_col) | |
| df = df.append(df_temp) | |
| df.T.to_csv("rocket_sizing_"+ rocket_name + ".csv", | |
| float_format="%.4f", encoding="SHIFT-JIS", header = False) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment