多段ロケットの最適質量配分(サイジング)問題の計算
# -*- 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