実践②：ガスジェネレータサイクルのロケットエンジン
題材と前提
題材はSpaceXのMerlin 1Dです。Falcon 9の1段目に9基使われている、液体酸素（LOX）とケロシンのガスジェネレータ（ガス発生器。以下GG）サイクルのエンジンです。
ポンプが送り出した推進剤（燃料と液体酸素）の一部をGGで燃やし、そのガスでタービンを回してポンプを駆動します。タービンとポンプを1本の軸でつないだものがターボポンプです。タービンを回したガスは捨てます。SpaceX特集の第3回では、この捨てる分が燃焼圧（主燃焼室の圧力）とともに増えることを、時間変化のない計算で確かめました。今回は燃焼圧や効率などの前提を第3回とそろえたうえで、起動や停止のときに何が起きるかを、燃焼終了まで時間を追って解いてみます。
| 項目 | 値 | 根拠 |
|---|---|---|
| 海面推力・ノズル開口比 | 845 kN（86 tf）・16 | Wikipedia（SpaceX Merlin） |
| 燃焼圧 | 12.9 MPa（約130 bar） | 第3回と同じ推定値（現行型は非公表。後述） |
| 主燃焼室の酸燃比（O/F。混合比とも呼ぶ） | 2.36 | 第3回と同じ仮定 |
| ポンプ効率 | 75% | 第3回と同じ仮定 |
| GGタービンの効率・入口温度 | 45%・900 K | 第3回と同じ仮定 |
| ポンプ吐出圧・GGの圧力 | 燃焼圧の1.4倍（18.1 MPa）・0.8倍（10.3 MPa） | 第3回と同じ仮定 |
| タービン出口圧・タンク圧（ポンプの入口圧） | 0.2 MPa（流量によらず一定）・0.3 MPa | 第3回と同じ仮定 |
| 主燃焼室の容積 | 39 L（特性長さ ＝ 容積 ÷ スロート面積 ＝ 1.0 m） | 今回の仮定 |
| ターボポンプの回転数 | 36,000 rpm | Wikipedia（SpaceX Merlin。どの型の値かは不明。後述） |
| 燃焼効率（c*効率）・ノズル効率 | 97.5%・98% | 今回の仮定 |
| GGの容積・軸の慣性モーメント | 4 L・0.06 kg·m² | 今回の仮定 |
| 起動トルク（スピンスタートを模擬） | 止まっているとき300 N·m。回転数とともに下がり、設計回転数で0。0.5秒で止める | 今回の仮定 |
| 起動前のGGの中のガス | 窒素（1気圧・300 K） | 今回の仮定 |
| ノズルの流れが壁からはがれる目安 | 壁の圧力が大気圧の0.4倍 | Wikipedia（Rocket engine nozzle）にある出口の圧力の目安を、壁の圧力に当てはめた |
スロートはノズルのいちばん細い部分、ノズル開口比は出口とスロートの面積比です。あとで出てくる比推力は、推進剤を毎秒1 kg使うときに何kgfの推力が出るかを表す値で、燃費に当たります（単位は秒）。
非公表の値がけっこう多いね。
エンジンの中身はほとんど公表されていないので、仮定はこうして表にしておきます。大事なのは絶対値より、何を変えるとどれだけ動くかです。
部品のモデル
部品ごとのモデルと、実際の開発ではその中身を何で決めるかを表にまとめます。
| 部品 | モデル | 実際の開発でのデータ源 |
|---|---|---|
| ポンプ | 揚程マップと軸動力マップ（相似則を使い、回転数が変わっても1組の曲線で済むように整理したもの） | ポンプ単体試験 |
| 弁（バルブ）・噴射器・オリフィス | 流量 \( \dot m = C_dA \)\( \sqrt{2\rho\,\Delta P} \) | 水流し試験 |
| 主燃焼室 | 容器の中のガスの質量保存（圧力）と、流れ込む推進剤の燃焼温度への追従（温度） | 容積は図面 |
| GG | 容器の中のガスの質量保存。温度は、中のガスのO/Fの燃焼温度 | 容積は図面 |
| 燃焼ガスの性質 | 酸燃比O/Fの表（燃焼温度、特性排気速度c*、推力係数） | 化学平衡計算、燃焼試験 |
| タービン | 速度比（翼の周速 ÷ ガスの噴出速度）で決まる効率マップ | タービン単体試験 |
| ターボポンプの軸 | \( J\,d\omega/dt = \)（タービンの出力 − ポンプの動力）\( /\omega \) ＋ 起動トルク | 慣性は図面 |
| ノズル | 推力 \( F = \eta\,C_F P_c A_t \)\( {}- P_a A_e \)（\( C_F \) は真空中の推力係数）。燃焼圧が下がって流れが壁からはがれるときは、はがれる位置までで計算する。タービン排気の推力も足す | ノズルの設計 |
式の記号は、\( C_d \) が流量係数、\( A \) が流路の面積、\( \rho \) が密度、\( \Delta P \) が前後の圧力差、\( J \) が軸の慣性モーメント、\( \omega \) が回転の角速度、\( P_c \) が燃焼圧、\( A_t \) と \( A_e \) がスロートとノズル出口の面積、\( P_a \) が大気圧、\( \eta \) がノズル効率です。
噴射器は推進剤を燃焼室に吹き込む部品、オリフィスは流量を決める絞り穴です。c*（特性排気速度）は同じ流量でどれだけ高い燃焼圧が立つか、推力係数はノズルが燃焼圧をどれだけ推力に変えるかを表します。
右の列には、ポンプやタービンの単体試験、噴射器の水流し試験といった、部品ごとの試験が並びます。前編で書いた、実測を部品として取り込むという考え方が、ロケットエンジンにもそのまま当てはまります。今回は試験データがないので、マップの形は典型的なものを仮定しました。
開発が進んで単体試験が1つ終わるたびに、仮定のマップを実測のデータに差し替えていきます。モデルが少しずつ実機に近づいていくのが、この仕事の楽しいところです！
こうした単体試験（要素試験）を何度も繰り返すのは、部品そのものを開発するためだけではありません。システム全体のチューニングに欠かせないパラメータを得るためでもあります。最後まで全部が仮定のままのロケットを飛ばすのは無理ですからね。
遠心ポンプが吸収する動力は、流量がゼロ（締切り）に近いところでは、設計点の半分ほどしかありません（今回は締切りで45%と仮定しました）。この性質は、次のページで主弁を開くのを遅らせたときの、回転数の行き過ぎに効いてきます。
モデルは2段階で作ります。まず時間変化のないつり合いの計算で、推力845 kN・燃焼圧12.9 MPaになるようにスロートとオリフィスの大きさを決めます（以下、設計計算）。次に、その大きさを固定して時間積分します。
酸燃比と燃焼温度
ロケットエンジンの性能と温度を大きく左右するのが酸燃比（O/F）です。液体酸素と燃料を、質量で何対何の割合で燃やすかを表します。同じ推進剤でも、この割合で燃焼温度は大きく変わります。
主燃焼室の燃焼温度は、O/F 2.36で約3,690 Kです（モデルでは、燃焼効率97.5%の2乗にあたる約0.95を掛けた、約3,500 Kを使います）。GGはO/F 0.32という極端な燃料過剰で燃やして900 Kに抑え、そのガスでタービンを回します。
O/F 0.8以上の値は、第3回と同じく、化学平衡計算（燃焼ガスが反応しきった状態を仮定して、温度や組成を出す計算）で求めました。計算には、NASAの計算コードCEAをPythonから呼び出せるRocketCEAを使っています。一方、GGのような極端な燃料過剰の領域では、平衡計算の結果が教科書の値と大きく食い違います。O/F 0.32では、平衡計算だと約1,200 Kになります。これに対して、Huzel & Huangの教科書『Design of Liquid Propellant Rocket Engines』（NASA SP-125）に載っている燃料過剰ガスの代表値は約900 Kです。そこでGGの領域には、この教科書の値を使っています。
\[ \begin{aligned} \frac{dT}{dt}&=\frac{T_{\mathrm{ad}}(O/F)-T}{\tau}\\[4pt] \tau&=\frac{m}{\dot m_{\mathrm{in}}},\quad m=\frac{PV}{RT} \end{aligned} \] \[ \begin{aligned} \frac{dP}{dt}&=\frac{RT}{V}\left(\dot m_{\mathrm{in}}-\dot m_{\mathrm{out}}\right)\\ &\quad+\frac{P}{T}\,\frac{dT}{dt}\\[4pt] \dot m_{\mathrm{out}}&=\frac{P\,A_t}{c^*} \end{aligned} \]
上の式は、容器の中のガスの温度 \( T \) が、流れ込む推進剤のO/Fの燃焼温度 \( T_{\mathrm{ad}} \) に、ガスが入れ替わる時間 \( \tau \) で追いつくことを表します。燃えてできたガスが中のガスと混ざっていく様子を簡単に表した式で、熱の出入りを厳密に追うエネルギーの式ではありません。
下の式は、理想気体の状態方程式 \( P=mRT/V \) を時間で微分したもので、ガスの質量保存を表します（ガス定数 \( R \) は一定としました）。
主燃焼室は、中のガスが約1.4 msで入れ替わるので、中身のO/Fは流れ込む推進剤のO/Fとほとんど変わりません。そのため、流れ込む推進剤のO/Fで燃焼温度を決めています。
これに対してGGは、弁を開け閉めする順番がずれると、中のガスのO/Fが、流れ込む推進剤のO/Fと大きく食い違います。そこでGGでは、中のガスの質量を3つに分けて追いかけ、温度は中身のO/Fから決めます。
\[ \begin{aligned} \frac{dm_{gg}}{dt}&=\dot m_f+\dot m_o-\dot m_t\\[4pt] \frac{dm_{o,gg}}{dt}&=\dot m_o-\frac{m_{o,gg}}{m_{gg}}\,\dot m_t\\[4pt] \frac{dm_{i,gg}}{dt}&=-\frac{m_{i,gg}}{m_{gg}}\,\dot m_t \end{aligned} \] \[ \begin{gathered} T=T_{\mathrm{ad}}\!\left(\frac{m_{o,gg}}{m_{f,gg}}\right)\\ P=\frac{m_{gg}RT}{V} \end{gathered} \]
\( m_{gg} \) はGGの中のガス全体の質量、\( m_{o,gg} \) はそのうち酸化剤からできた分、\( m_{i,gg} \) は起動前から入っていた窒素です（構成図の3つの状態量）。全体からこの2つを引いた残りの \( m_{f,gg} \) が、燃料からできた分です。\( \dot m_f \)、\( \dot m_o \) はGGに流れ込む燃料とLOX、\( \dot m_t \) はタービンへ出ていくガスの流量で、出ていくガスの組成は中身と同じとします。
温度は、推進剤からできたガスのO/F（酸化剤からできた分 ÷ 燃料からできた分）の燃焼温度です。あとから流れ込んだ酸素は、中に残っている燃料過剰のガスとも燃える、と考えています。窒素が残っているあいだは、その分だけ温度が下がります（式では省略しました）。
どちらの温度の式にも、圧力が急に変わるときの圧縮や膨張による温度変化と、壁への熱の逃げは入っていません。
コードでは、これが次の形になります。
def derivatives(self, t, s):
f = self.flows(t, s) # ポンプの動作点・各流量・動力
# 軸の運動方程式
torque = (f["W_t"] - f["W_f"] - f["W_o"]) / s["omega"] + self.spin(t, s["omega"])
# 主燃焼室（圧力と温度）
T_c = self.d["eta_cstar"] ** 2 * T_AD(f["of_c"])
dPc, dTc = self.chamber(s["Pc"], s["Tc"], self.V_c, self.R_c,
f["mf_c"] + f["mo_c"], f["m_c_out"], T_c)
# GG（ガスの質量保存）: タービンへ出ていくガスは、中身と同じ組成
out = f["m_t"] / s["m_gg"]
return {"m_f": -(f["mf_c"] + f["mf_g"]), # タンクの質量保存
"m_o": -(f["mo_c"] + f["mo_g"]),
"omega": torque / self.J,
"Pc": dPc, "Tc": dTc,
"m_gg": f["mf_g"] + f["mo_g"] - f["m_t"], "mo_gg": f["mo_g"] - out * s["mo_gg"], "mi_gg": -out * s["mi_gg"]}タンクの式は、揚水発電所の池とまったく同じ形だ。部品が増えても、やることは変わらないんだね。
起動過渡
時間刻み0.2 msで、起動からの3秒間を計算しました。計算時間は数秒です。実機の起動手順は公表されていないので、手順は次のように仮定しました。
- 0秒：スピンスタート。外部からのガスでタービンを回す方式を、回転数とともに下がる起動トルクで模擬する（止まっているとき300 N·m、設計回転数で0。0.5秒で止める）
- 0.15秒：GGの燃料弁とLOX弁を開く（0.1秒で全開）
- 0.2〜0.5秒：主燃料弁と主LOX弁（以下、主弁）を開いていく
この手順は、回転数や圧力が設計値を行き過ぎないように調整したものです。主弁を早めに開いて、回転が上がりきる前にポンプに負荷をかけます。
回転数は0.52秒で設計回転数の99%に達し、行き過ぎることなく36,000 rpmに落ち着きました。
ポンプの吐出圧、GGの圧力、主燃焼室の圧力は、そろって上がっていきます。吐出圧も、設計値の18.1 MPaを超えません。
推力は0.46秒で定常値の90%に達し、0.6秒までには845 kNに落ち着きます。燃料弁とLOX弁を同時に開いているので、起動中もO/Fはほぼ設計値のままです。
思ったよりあっさり立ち上がるんだね。
実機でも、弁を開く順番と速さを決めて、各段階を監視しながら起動します。その順番や速さが少しずれるとどうなるかは、次のページでわざとずらして確かめます。
答え合わせ：Merlin 1Dの公表値と比べる
前編で、Merlin 1Dの公表値が1Dモデルの答え合わせになると書きました。設計計算の結果を比べます。
| 項目 | モデル | 公表値 |
|---|---|---|
| 海面推力 | 845 kN（合わせた値） | 845 kN |
| 燃焼圧 | 12.9 MPa（推定値を入れた） | 非公表（2011年の設計目標は9.7 MPa） |
| 海面比推力 | 282 s（278〜285 s） | 282 s |
| 真空比推力 | 306 s（302〜309 s） | 311 s |
| 真空推力 | 916 kN | 914 kN |
| 推進剤の流量 | 306 kg/s | 約306 kg/s（推力と比推力から逆算） |
| ターボポンプの動力 | 7.2 MW | 約7.5 MW（どの型の値かは不明。後述） |
| ノズル出口の直径 | 0.89 m | 非公表（エンジンの直径は0.92 m） |
| GGに回す推進剤 | 全体の7.1%（5.8〜8.8%） | 非公表 |
| スロート直径 | 22.3 cm | 非公表 |
公表値はWikipediaに載っている値です。（ ）は、第3回と同じ幅（タービン効率40〜50%、入口温度880〜920 K、ポンプ効率72〜78%）で前提を変えたときの範囲です。GGに回す推進剤について、第3回で出した「100バールで約3.9%」は、メタンと液体酸素で計算した値です。今回はケロシンと液体酸素なので、同じ値にはなりません。
推力は、スロートの大きさとGGの流量を決めるときに合わせた値です。現行型の燃焼圧は公表されていないので、第3回と同じ推定値を入れました。2011年に発表された設計目標（燃焼圧 約1,410 psi＝9.72 MPa、真空推力 約690 kN。Aviation Week）から、燃焼室とノズルの大きさは変わらず、燃焼圧が現行型の真空推力914 kNに比例して上がったと考えた値です（9.72×914÷690＝約12.9 MPa）。
答え合わせになるのは、合わせていない値です。海面比推力は282 sで公表値と同じ、真空比推力は約2%低く出ました。ターボポンプの動力は7.2 MW（公表値 約7.5 MW）で、ノズル出口の直径0.89 mも、エンジンの直径0.92 mに収まります。真空推力は916 kN（公表値914 kN）です。ただし、燃焼圧の推定に真空推力914 kNを使っているので、真空推力が近い値になるのは、独立した答え合わせとは言えません。
海面はぴったりなのに、真空だけ2%ほど低いんだね。どこがずれてるんだろう？
仮定を見直す：タービンとポンプ
比推力を大きく左右する要因のひとつが、GGに回して捨てる推進剤の量です。まずタービンの効率の仮定を変えて、設計計算をやり直しました。
タービンの効率を上げると、GGに回す量が減って比推力が上がります。海面比推力は、基準の45%のときに公表値の282 sとほぼ同じで、効率をそれより高く見ると公表値を超えます。
次に教科書の値を見てみます。Huzel & Huangの教科書には、同じLOX/ケロシンのGGサイクルの例題エンジンが載っていて、ターボポンプの値は次のとおりです。
| 仮定 | 今回の基準 | 教科書の例題エンジン |
|---|---|---|
| タービンの効率 | 45% | 58% |
| タービン入口温度 | 900 K | 1,033 K |
| ポンプ吐出圧 ÷ 燃焼圧（LOX側／燃料側） | 1.4／1.4 | 1.5／1.7 |
| ポンプ効率（LOX側／燃料側） | 75%／75% | 71%／66% |
例題のタービンは、ノズルで加速したガスの勢いを、2列の動翼で2回に分けて受け止める形式（速度複式衝動タービン）です。NASAの設計基準書（SP-8110）によると、GGサイクルなどで使われる2列のタービンの効率は35〜65%です。基準の45%と例題の58%は、この範囲の低めと高めに当たります。Merlin 1Dの値は公表されていません。
教科書の例題エンジンで燃料側の吐出圧が高い（燃焼圧の1.7倍）のは、燃料が燃焼室の冷却通路を通ってから噴射器に届くためです。この4項目を教科書の値に置き換えると、比推力は海面287 s・真空311 sになりました。真空は公表値と同じになりますが、今度は海面が約2%高くなります。GGに回す推進剤は5.0%に減ります。
基準だと海面が合って、教科書の値だと真空が合うんだ。両方いっぺんには合わないの？
合わない理由のひとつは、公表値どうしの食い違いです。推進剤の流量が同じなら、比推力の比（311÷282＝1.10）と推力の比（914÷845＝1.08）は一致するはずですが、2%ほどずれています。比推力の値には、Wikipediaでも「要更新」の注記が付いています。Falcon 9 v1.1（2013年）のページにも同じ282 s・311 sが載っていて、推力を上げたあとの値かどうかは分かりません。
燃焼効率やノズル効率の仮定でも比推力は1%程度動くので、これでエンジンの中身を特定できたわけではありません。分かったのは、公表値と大きくは矛盾しない仮定の組み合わせがある、というところまでです。
燃焼圧の前提を変えた場合
燃焼圧の推定値も確かめておきます。現行型の値は公表されていないので、次の2通りを比べます。
- 12.9 MPa（約130 bar）：今回の基準（第3回と同じ推定値）
- 9.7 MPa：2011年に発表された設計目標の値のまま、と考えた場合（当時の目標は真空推力 約690 kN）
燃焼圧が9.7 MPaのままだと、海面845 kNを出すノズルの出口の直径は1.04 mになり、Wikipediaに載っているエンジンの直径0.92 mに収まりません。このモデルで出口が0.92 mに収まるのは、燃焼圧が約12.2 MPa以上のときです。
推力の公表値だけから逆算することもできます。2016年に発表された増強後の推力は、海面845 kN・真空914 kNです。この差69 kNは「大気圧 × 出口の面積」（主ノズルと、タービン排気の出口の合計）に当たるので、出口の面積が決まり、開口比16からスロートの大きさも決まります。このスロートで海面845 kNを出すのに必要な燃焼圧は、このモデルでは約13.3 MPaです。推定値の12.9 MPaとの差は3%ほどです。
表と図は、12.9 MPaと9.7 MPaの2通りの燃焼圧と、ターボポンプの2通りの仮定（今回の基準と、教科書の例題エンジン。表では「基準」「例題」と書きます）を組み合わせた、4通りの設計計算の結果です。
| 燃焼圧 [MPa] | ターボポンプ の仮定 | 比推力 [s] 海面／ 真空 | GGに回す 推進剤 [%] | ターボポンプ の動力 [MW] |
|---|---|---|---|---|
| 12.9 | 基準 | 282／306 | 7.1 | 7.2 |
| 12.9 | 例題 | 287／311 | 5.0 | 8.8 |
| 9.7 | 基準 | 278／309 | 5.7 | 5.4 |
| 9.7 | 例題 | 283／313 | 4.0 | 6.6 |
| 公表値 | — | 282／311 | 非公表 | 約7.5 |
ターボポンプの動力の公表値は、Wikipediaに載っている36,000 rpm・10,000馬力（約7.5 MW）です（出典はElon Musk氏の2018年のインタビュー）。Wikipediaでは初期型（Merlin 1A〜1C）のターボポンプの説明に続けて書かれていて、どの型の値かははっきりしません。
4通りとも、比推力と公表値の差は約2%以内です。ターボポンプの動力は5.4〜8.8 MWで、公表値の約7.5 MWをはさみます。燃焼圧を上げると、ターボポンプの動力は燃焼圧にほぼ比例して増えます。海面の比推力は大気圧の影響が小さくなるぶん上がり、真空の比推力はGGに回す推進剤が増えるぶん下がります。
比推力と動力だけでは、4通りのどれが実機に近いかまでは絞り込めません。違いがはっきり出るのはノズルの大きさで、9.7 MPaでは出口がエンジンの直径に収まりません。なお、Wikipediaのロケットエンジン比較表には、個人のサイトを出典として108 bar（10.8 MPa）とありますが、そのサイトに根拠の記載はありません（このモデルでは、10.8 MPaでも出口は0.92 mに収まりません）。
このあとのパラメータスタディは、12.9 MPaを基準にした今回のモデルで続けます。
公表値を集めるときは、それがいつの時点の値なのかも確かめる必要がありますね（汗）
※このページの値は、仮定を置いた簡易モデルの計算結果です。前提の表のうち非公表の値（燃焼圧は第3回と同じ推定値）と、起動の手順は仮定で、実機のMerlin 1Dの値や手順を表すものではありません。
次ページへ
次のページでは、このモデルの条件を変えて、エンジンが何に敏感なのかを調べます。
