RocketCEAでロケットエンジン性能計算

概要

ロケットエンジンの化学平衡計算と呼ばれる性能計算をするために便利なRocketCEAというpythonモジュールの紹介。従来からよく使われているNASA CEAという難解なアウトプットのものに比べ、取り扱いしやすい。
詳細は公式ドキュメントやGithubのページへ

NASA CEA

化学平衡計算は初歩的なロケットの勉強から開発の最前線まで多用される。
化学平衡計算は従来からロケットエンジン設計者なら全員が使ってきたNASA開発のCEA（Chemical Equilibrium with Applications)がある。
これは論文もあり、内容は流石NASAというべきすごいもの。しかし作られたのが古く使いやすさの面ではイマイチ。
加えて、最近ではソースコードが非公開になり（探せばネット上落ちてはいるが見つけにくい）、オンラインで実行できるサイトだけが残っている。

また、NASA CEAはFORTRANで書かれており、ソースコードは難読かつ、出力は古い形式のテキストで出てくるので外部のプログラムからの取り扱いは難しい。

個人的にはこの難しいNASA CEAを幅広い人が使えるようにとコンパイルして外部のpythonスクリプトから呼び出す記事を書いたり、NASA CEAが非公開になったときにTwitterで会話していたりしていた。

実際の開発現場でもNASA CEAは社内ツールの中に多数組み込まれている。

RocketCEA

NASA CEAは中身が充実しているツールであるが、その中でもロケットの性能計算部分だけのラッパーとしてのpythonモジュールがRocketCEAである。2018年頃から開発されてきたいたらしい。私は2019年初頭にRocketCEAの存在を知った。当時は機能的にまだまだでスルーしていたが、2021年では改良されていて実用的になっている。

インストール

公式サイトを見るのが良い。Windowsではgfortranをインストールするのが多少めんどくさい。他のOSやWindowsでもWSLを使うとインストールはコマンド数回入れるだけで完了するぐらい簡単。

【追記】windowsでのインストールについて詳細記事でました。
https://makkiblog.com/rocketcea_intro/

サンプルコード

NASA CEAの結果を出力

まずは一番シンプルな例。NASA CEAと同じ結果を出力するサンプルコード。
Pcが燃焼室圧力、MRが混合比でO/Fとも言われるもの、epsがスロートとノズルの開口比。
RocketCEAの公式はヤード・ポンド法がデフォルトとして使われいるのでサンプルコードでは国際的な標準であるSI単位系で見れるようにしている。

sample01.py

from rocketcea.cea_obj import CEA_Obj 

ispObj = CEA_Obj( oxName='LOX', fuelName='CH4')
str = ispObj.get_full_cea_output( Pc=100.0, MR=6.0, eps=40.0, short_output=1, pc_units='bar', output='siunits')
# str = ispObj.get_full_cea_output( Pc=1000.0, MR=6.0, eps=40.0, frozen=1, frozenAtThroat=1, pc_units='bar')  # 凍結流の場合はこういうオプションをつける

print(str)

必要な値を出力

設計するのに必要な各種データを出力する。get_ホニャララの部分で様々な値が取得できる。また、それ以外にも公式の関数集を見れば何の値を取得できるかわかる。

sample02.py

import numpy as np
from rocketcea.cea_obj_w_units import CEA_Obj

ceaObj = CEA_Obj( oxName='LOX', fuelName='C2H5OH', 
                  pressure_units='MPa', cstar_units='m/s', temperature_units='K')

print(ceaObj.desc)

# ---- 燃焼ガス特性などを見る(平衡流) ----
def show_perf(Pc=1.1, eps=2.3, MR=1.45):
    IspVac, Cstar, Tc, MW, gamma = ceaObj.get_IvacCstrTc_ChmMwGam(Pc=Pc, MR=MR, eps=eps)
    print( '%5.1f %8.1f   %6.1f  %8.1f    %8.1f    %8.1f  %8.2f  %8.4f '%\
         (Pc, eps, MR, IspVac, Cstar, Tc, MW, gamma))

print(' Pc(MPa) AreaRatio  O/F   IspVac(sec)  Cstar(m/s)  Tc(degK)   MolWt    gamma')
for MR in np.arange(1, 2, 0.1):
    show_perf(Pc=1.1, eps=2.3, MR=MR)

# ---- オマケ：Ispや⊿Vを確認する ----
Wpayload = 30.0 # [kg]
Wstg = 350.0 # [kg]
Wpropellant = 800.0 # [kg]
effIsp = 0.90 # Isp efficiency

Winit = Wstg + Wpropellant + Wpayload
Wfinal = Winit - Wpropellant

def show_deltaV(Pc=1.1, eps=2.3, MR=1.45):
    IspVacTheoritical = ceaObj.get_Isp(Pc=Pc, MR=MR, eps=eps)
    IspVacMeasured = effIsp * IspVacTheoritical
    IspSL, mode = ceaObj.estimate_Ambient_Isp(Pc=Pc, MR=MR, eps=eps, Pamb=0.1, frozen=0, frozenAtThroat=0)
    IspSLMeasured = effIsp * IspSL
    deltaV = 9.80665 * IspVacMeasured * np.log( Winit / Wfinal ) # [m/s]
    print( '%5.1f %8.1f %7.2f  %8.1f    %8.1f     %8.1f   %8.1f    %8.1f'%(Pc, eps, MR, IspVacTheoritical, IspVacMeasured, IspSL, IspSLMeasured, deltaV))
    # print(mode)

print()
print('Pc(MPa) AreaRatio  O/F  IspVac(s)  IspVac実測(s)  IspSL(s) IspSL実測(s)  deltaV(m/s)')

for eps in np.arange(1, 4, 0.5):
    show_deltaV(Pc=1.1, eps=eps, MR=1.45)

リストにない推進剤を使う

RocketCEAでは様々なロケットの推進剤をデフォルトで使うことができるが、リストに無いものはNASA CEAと同様にカードを作ることで使えるようになる。カードには元素と生成エンタルピー、温度、密度が書いてあれば良い。

sample03.py

# cf. https://rocketcea.readthedocs.io/en/latest/new_propellants.html

from rocketcea.cea_obj import CEA_Obj, add_new_fuel, add_new_oxidizer, add_new_propellant

card_str = """
name H2O2(L) H 2 O 2  wt%=100.00
h,cal=-44880.0     t(k)=298.15  rho.g/cc=1.407
"""

add_new_propellant( 'MyProp', card_str )
ceaObj = CEA_Obj(propName="MyProp")

str = ceaObj.get_full_cea_output( Pc=250.0, eps=40.0, short_output=1, output='siunits')

print(str)

c*効率をグラフで可視化する

RocketCEAであれば、グラフ化も非常にやりやすい。
これまではグラフ化するのが本当に大変だったのでありがたい。。。

sample04.py

import numpy as np
import matplotlib.pyplot as plt
from rocketcea.cea_obj_w_units import CEA_Obj 

ceaObj = CEA_Obj(propName='', oxName='LOX', fuelName="CH4", cstar_units='m/s', pressure_units='MPa')

pcArr = [ 1, 5, 10, 30]

for Pc in pcArr:
    mrArr = np.arange(1.5, 4.5, 0.05)
    cstarArr = [ceaObj.get_Cstar(Pc=Pc, MR=MR) for MR in mrArr]
    # 上のリスト内包表記に慣れていない人は下のコメント行と同じと思ってもらうと比較的わかりやすい
    # cstarArr = []
    # for MR in mrArr:
    #     cstarArr.append( ceaObj.get_Cstar( Pc=Pc, MR=MR) )
    plt.plot(mrArr, cstarArr, label='Pc=%g MPa' % Pc)

plt.legend(loc='best')
plt.grid(True)
plt.title(ceaObj.desc)
plt.xlabel('Mixture Ratio')
plt.ylabel('Cstar (m/s)')
plt.show()

Ispをグラフで可視化する

当然、重要なIspの出力をすることも可能。真空中比推力や与えた外気雰囲気圧力下での比推力も出力可能。

sample05.py

import numpy as np
import matplotlib.pyplot as plt
from rocketcea.cea_obj_w_units import CEA_Obj 

ceaObj = CEA_Obj(propName='', oxName='LOX', fuelName="CH4", temperature_units='K', pressure_units='MPa')

Pc = 10.0
area_ratio = [50.0, 20.0, 10.0]

for eps in area_ratio:
    mrArr = np.arange(1.5, 4.5, 0.05)
    ispArr = [ceaObj(Pc, MR, eps ) for MR in mrArr]
    # ispArr = [ceaObj.estimate_Ambient_Isp(Pc, MR, eps, Pamb=0.1013)[0] for MR in mrArr]
    plt.plot(mrArr, ispArr, label='AreaRatio=%g'%eps)

plt.legend(loc='best')
plt.grid(True)
plt.title(ceaObj.desc + ",  Pc=%dMPa"%Pc + ', equilibrium')
plt.xlabel('Mixture Ratio')
plt.ylabel('Isp vac(sec)')
plt.show()

有限体積燃焼室、平衡流、凍結流の違いを見る

NASA CEAを触ったことある、もしくは論文を読んでいる人知っている通り、NASA CEAでは燃焼室体積が有限か、化学反応を途中で凍結frozenさせているかどうかで値が異なる。これを間違うと性能見積もりを誤るので重要な点。
デフォルトは無限体積燃焼室で平衡流になっているので注意。

sample06.py

import numpy as np
import matplotlib.pyplot as plt
from rocketcea.cea_obj_w_units import CEA_Obj 

Pc = 5.0
area_ratio = 10.0
chamber_contraction_ratio = 2.5
mrArr = np.arange(1.5, 4.5, 0.05)

ispObj = CEA_Obj(propName='', oxName='LOX', fuelName="CH4", fac_CR=chamber_contraction_ratio, temperature_units='K', pressure_units='MPa')

plt.figure()

ispArr = [ispObj.get_Isp(Pc, MR, area_ratio, frozen=0, frozenAtThroat=0) for MR in mrArr]
plt.plot(mrArr, ispArr, label='equilibrium')
ispArr = [ispObj.get_Isp(Pc, MR, area_ratio, frozen=1, frozenAtThroat=1) for MR in mrArr]
plt.plot(mrArr, ispArr, label='frozen in throat')
ispArr = [ispObj.get_Isp(Pc, MR, area_ratio, frozen=1, frozenAtThroat=0) for MR in mrArr]
plt.plot(mrArr, ispArr, label='frozen in chamber')

plt.legend(loc='best')
plt.grid(True)
plt.title(ispObj.desc + " ,Pc=%.1fMPa, AR=%.1f" % (Pc, area_ratio))
plt.xlabel('Mixture Ratio')
plt.ylabel('Isp vac(sec)')
plt.show()

いろんなグラフを出す

設計に必要そうなグラフの出し方の例を挙げておく。

sample07.py

import numpy as np
import matplotlib.pyplot as plt
from rocketcea.cea_obj_w_units import CEA_Obj 

ispObj = CEA_Obj(propName='', oxName='LOX', fuelName="CH4", temperature_units='K', pressure_units='MPa')
plt.ion()
# -----------------

Pc = 10.0
MR = 3.5

plt.figure()
epsArr = np.arange(5, 60, 1)
ispVacArr = [ispObj.get_Isp(Pc, MR, eps) for eps in epsArr]
ispSLArr = [ispObj.estimate_Ambient_Isp(Pc, MR, eps, Pamb=0.1013)[0] for eps in epsArr]
plt.plot(epsArr, ispVacArr, label='Isp vac')
plt.plot(epsArr, ispSLArr, label='Isp SL')
plt.legend(loc='best')
plt.grid(True)
plt.title(ispObj.desc + "  Pc=%dMPa"%Pc)
plt.xlabel('Nozzle Area Ratio(ε)')
plt.ylabel('Isp (sec)')

# -----------------
MR = 3.5
epsArr = [10, 20, 30, 40, 50]

plt.figure()
for eps in epsArr:
    pcArr = np.arange(1, 30, 1)
    ispVacArr = [ispObj.get_Isp(Pc, MR, eps) for Pc in pcArr]
    plt.plot(pcArr, ispVacArr, label='eps=%d'%eps)
plt.legend(loc='best')
plt.grid(True)
plt.title(ispObj.desc)
plt.xlabel('chamber pressure(MPa)')
plt.ylabel('vacuum Isp (sec)')

# -----------------
MR = 3.5
epsArr = [10, 20, 30, 40, 50]

plt.figure()
for eps in epsArr:
    pcArr = np.arange(1, 30, 1)
    ispSLArr = [ispObj.estimate_Ambient_Isp(Pc, MR, eps, Pamb=0.1013)[0] for Pc in pcArr]
    plt.plot(pcArr, ispSLArr, label='eps=%d'%eps)
plt.legend(loc='best')
plt.grid(True)
plt.title(ispObj.desc)
plt.xlabel('chamber pressure(MPa)')
plt.ylabel('sea level Isp (sec)')

# -----------------
MR = 3.5
epsArr = [10, 20, 30, 40, 50]

plt.figure()
for eps in epsArr:
    pcArr = np.arange(1, 30, 1)
    peArr = [Pc / ispObj.get_PcOvPe(Pc, MR, eps) for Pc in pcArr]
    plt.plot(pcArr, peArr, label='eps=%d'%eps)
plt.legend(loc='best')
plt.grid(True)
plt.title(ispObj.desc)
plt.xlabel('chamber pressure(MPa)')
plt.ylabel('Nozzle exit pressure(MPa)')
plt.ylim([0, 0.12])

plt.show()

固体ロケット（その１）

固体ロケットの場合はモノプロペラントのように各組成を記述していく。wtのところに重量割合を入れる。

sample08.py

from rocketcea.cea_obj import CEA_Obj, add_new_fuel, add_new_oxidizer, add_new_propellant

card_str = """
name NH4CLO4(I) wt=100 t,k=298
name AL(cr) wt=25 t,k=298
name HTPB wt=75 t,k=298
h,kj/mol=23.99 C 4 H 6 O 0.04 H 0.04 
"""
add_new_propellant( 'MySolid', card_str )
ceaObj = CEA_Obj(propName="MySolid")

s = ceaObj.get_full_cea_output( Pc=50.0, eps=3.0, short_output=1, pc_units='bar', output='siunits')
print(s)

固体ロケット（その２）

固体ロケットでのグラフ化の例

sample09.py

import numpy as np
import matplotlib.pyplot as plt
from rocketcea.cea_obj import add_new_fuel, add_new_oxidizer, add_new_propellant
from rocketcea.cea_obj_w_units import CEA_Obj

card_str = """
name NH4CLO4(I) wt=100 t,k=298
name AL(cr) wt=25 t,k=298
name HTPB wt=75 t,k=298
h,kj/mol=23.99 C 4 H 6 O 0.04 H 0.04 
"""
add_new_propellant( 'MySolid', card_str )
ceaObj = CEA_Obj(propName="MySolid", cstar_units='m/s', pressure_units='MPa')

pcArr = np.arange(0.5, 15, 0.05)
cstarArr = [ceaObj.get_Cstar(Pc=Pc) for Pc in pcArr]
plt.plot(pcArr, cstarArr)
plt.grid(True)
plt.xlabel('Chamber pressure (MPa)')
plt.ylabel('Cstar (m/s)')
plt.show()

固体ロケット（その３）

固体ロケットで組成を変えるときは少し工夫が必要。新しくカードを作り、その割合を変えていく。

sample10.py

import numpy as np
import matplotlib.pyplot as plt
from rocketcea.cea_obj import add_new_fuel, add_new_oxidizer, add_new_propellant
from rocketcea.cea_obj_w_units import CEA_Obj

def make_ceaObj_from_newprop(ox_ratio=0):
    wtOxid = ox_ratio
    wtFuel = 100 - ox_ratio
    card_str = """
    name NH4CLO4(I) wt=%d t,k=298
    name HTPB wt=%d t,k=298
    h,kj/mol=23.99 C 4 H 6 O 0.04 H 0.04 
    """ % (wtOxid, wtFuel)
    prop_name = 'MySolid' + str(ox_ratio) # 名前は新しくつけないといけない
    add_new_propellant(prop_name, card_str )
    ceaObj = CEA_Obj(propName=prop_name, 
                     cstar_units='m/s', pressure_units='MPa', temperature_units='K')
    return ceaObj

ox_ratioArr = np.arange(70, 100, 1)
Pc = 2.0 # (MPa)
eps = 10

ispArr = []
cstarArr = []
tcombArr = []
for i in ox_ratioArr:
    ceaObj = make_ceaObj_from_newprop(i)
    cstarArr.append(ceaObj.get_Cstar(Pc=Pc))
    ispArr.append(ceaObj.get_Isp(Pc=Pc, eps=eps))
    tcombArr.append(ceaObj.get_Tcomb(Pc=Pc))

plt.figure()
plt.plot(ox_ratioArr, cstarArr)
plt.grid(True)
plt.title("Pc=%.1f MPa" % (Pc))
plt.xlabel('Oxidizer AP ratio (%)')
plt.ylabel('Cstar (m/s)')

plt.figure()
plt.plot(ox_ratioArr, ispArr)
plt.grid(True)
plt.title("Pc=%.1f MPa, eps=%.1f" % (Pc, eps))
plt.xlabel('Oxidizer AP ratio (%)')
plt.ylabel('Isp vac (sec)')

plt.figure()
plt.plot(ox_ratioArr, tcombArr)
plt.grid(True)
plt.title("Pc=%.1f MPa" % (Pc))
plt.xlabel('Oxidizer AP ratio (%)')
plt.ylabel('Combustion chamber gas temperature (K)')

plt.show()

推進剤比較

RocketCEAの公式にも推進剤比較の図があるが、単位系がヤード・ポンド法で非常にイケてない。自分で書くためにはこうすればいいという例を挙げておく。
各推進剤での比重と比推力Ispとの関係をグラフにすることによってロケットシステムにおいての推進剤の選択肢がどういう意味を持つのかの一面が理解できる。例えば、液酸液水（LOX/LH2）は比推力は高いが、比重が小さすぎる。

sample11.py

import numpy as np
import matplotlib.pyplot as plt
from rocketcea.cea_obj_w_units import CEA_Obj
from rocketcea.biprop_utils.density_at_mr import bulkDensity
from rocketcea.biprop_utils.mr_t_limits import MR_Temperature_Limits

def plot_Isp_from_prop(oxName="LOX", fuelName="LH2", Pc=5.0, eps=15):
    ceaObj = CEA_Obj(propName='', oxName=oxName, fuelName=fuelName, cstar_units='m/s', pressure_units='MPa')
    mrObj = MR_Temperature_Limits(oxName=oxName, fuelName=fuelName, MR_MAX=20.0)
    if mrObj.min_MR < 1:
        mrObj.min_MR = 1 
    mrArr = np.arange(mrObj.min_MR, mrObj.max_MR, 0.1)
    sgArr = bulkDensity(oxName, fuelName, mrArr) * 1000 # (g/mL) -> (kg/m3)
    ispArr = [ceaObj.get_Isp(Pc=Pc, MR=MR, eps=eps) for MR in mrArr]
    plt.plot(sgArr, ispArr, label="%s / %s"%(oxName, fuelName))

plt.figure()

plot_Isp_from_prop(oxName="LOX", fuelName="LH2")
plot_Isp_from_prop(oxName="LOX", fuelName="CH4")
plot_Isp_from_prop(oxName="LOX", fuelName="RP1")
plot_Isp_from_prop(oxName="LOX", fuelName="Ethanol")
plot_Isp_from_prop(oxName="H2O2", fuelName="RP1")
plot_Isp_from_prop(oxName="N2O4", fuelName="UDMH")
plot_Isp_from_prop(oxName="N2O4", fuelName="MMH")
plot_Isp_from_prop(oxName="IRFNA", fuelName="UDMH")

plt.legend(loc='best', fontsize=9)
plt.xlabel('specific gravity of propellant bulk (kg/m3)')
plt.ylabel('CEA equilibrium Isp vac (sec)')
plt.xlim([100, 1450])
plt.ylim([240, 450])
plt.yticks(np.arange(240, 450, 20))
plt.title("Pc = 5.0 MPa, nozzle expansion ratio=15")
plt.grid(True)
plt.show()