Skip to content

Instantly share code, notes, and snippets.

@ina111

ina111/rocket_sizing.py

Last active August 6, 2022 05:32
Show Gist options
  • Save ina111/f4de3b3e27b88a4628f4574265f852cb to your computer and use it in GitHub Desktop.
Save ina111/f4de3b3e27b88a4628f4574265f852cb to your computer and use it in GitHub Desktop.
Download ZIP
多段ロケットの最適質量配分（サイジング）問題の計算
Raw
rocket_sizing.py
# -*- 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