Pythonで1Dシミュレーション入門｜揚水発電所とロケットエンジンを自作ソルバーで解いてみる

Contents

時間刻みと計算の発散
燃焼振動：数値の振動と本物の振動
ツールを使ってみよう
- 最後に
参考図書・引用一覧
- 揚水発電所

時間刻みと計算の発散

ロケットエンジンのモデルを、時間刻み2 msで、3つの積分法で計算してみます。揚水発電所は60秒刻みで問題なかったので、2 msなら十分細かそうに見えます。

同じモデルと同じ時間刻み2 msでの、オイラー法・2段ルンゲ＝クッタ法・4段ルンゲ＝クッタ法の比較 — 同じモデル・同じ時間刻み2 msでの、積分法による違い

ところが、オイラー法と2段ルンゲ＝クッタ法は、燃焼圧が上がっていく途中の0.3〜0.4秒で、正しい解から外れ始めます。オイラー法は1ステップごとに値が上下に大きく振れ続けます（下側が大気圧で止まっているのは、モデルに入れた圧力の下限という歯止めのためです）。2段ルンゲ＝クッタ法は上に外れたあと一気に跳ね上がり、0.39秒で、ポンプの動作点が求められなくなって計算が止まりました。計算の誤差がステップごとに増幅されて、正しい解に戻れなくなることを、計算の発散と呼びます。正しく解けたのは4段ルンゲ＝クッタ法だけです。

同じモデルなのに、計算方法が違うだけで発散するの？

原因は、燃焼圧が変化する速さ（時定数）です。主燃焼室の容積は39 Lしかないのに、毎秒約284 kgのガスが通り抜けます。2ページ目の主燃焼室のモデルに出てきた、中のガスが入れ替わる時間 \( \tau \) に数字を入れると、次のとおりです。

\[ \begin{aligned} \tau&=\frac{m}{\dot m}=\frac{P_cV_c/(R\,T_c)}{\dot m}\\ &=\frac{0.39\ \mathrm{kg}}{284\ \mathrm{kg/s}}\approx1.4\ \mathrm{ms} \end{aligned} \]

さらに、圧力が上がると噴射器の差圧（前後の圧力差）が小さくなって流入が減る、というブレーキも働くので、圧力の変化はもっと速くなります。そこで、タンクの残量を除く6つの状態量について、状態方程式を定常点（つり合った状態）のまわりで線形化し、固有値を求めました。固有値は、状態量の組み合わせ（モード）ごとに、定常値に戻る速さを表す値で、逆数が時定数です。いちばん速いものは \( \lambda \approx -1{,}170\ \mathrm{s^{-1}} \)、時定数で約0.85 msでした。

今の傾きから先を予測する方法（陽解法）が発散しない条件は、この固有値と時間刻みの積で決まります。

安定に解ける条件（固有値が負の実数の場合）

オイラー法と2段ルンゲ＝クッタ法：

\[ \begin{gathered} \Delta t\,|\lambda|\le 2\\ \Rightarrow\quad \Delta t\le\frac{2}{1{,}170}\approx1.7\ \mathrm{ms} \end{gathered} \]

4段ルンゲ＝クッタ法：

\[ \begin{gathered} \Delta t\,|\lambda|\le 2.785\\ \Rightarrow\quad \Delta t\le\frac{2.785}{1{,}170}\approx2.4\ \mathrm{ms} \end{gathered} \]

実際に時間刻みを変えて計算すると、オイラー法と2段ルンゲ＝クッタ法は1.7 msまで安定で2.0 msで発散、4段ルンゲ＝クッタ法は2.3 msまで安定で2.5 msで発散しました。見積もりと合っています。最初の図の発散は物理現象ではなく、数値計算が作ったものです。

2段ルンゲ＝クッタ法の誤差は、時間刻み0.2 msではオイラー法の約1/26、1 msでも約1/6です。それでも発散する限界はオイラー法と同じです。誤差の小ささ（精度）と発散しないこと（安定性）は別の性質なので、分けて確認する必要があります。

精度が高い方法なら、発散もしにくいんだと思ってたよ。

また、図の0.2〜1 msの範囲では、オイラー法の誤差は時間刻みにほぼ比例し（1次の精度）、2段と4段の誤差は時間刻みの2乗にほぼ比例しています（2次の精度。時間刻みを半分にすると誤差が1/4。斜めの破線が、それぞれの傾きです）。4段が4次に届かないのは、主弁が開いて燃焼圧が立ち上がり始める点など、式が滑らかでない箇所がモデルに残っているためです。誤差が最大になるのも、主弁を開き始める0.2秒の直後です。1ページ目で誤差の減り方を「式が滑らかな場合の理論値」と書いたのは、このためです。

それより大きい時間刻み（1.5 msや1.7 ms）では、4段の誤差が不規則に上下しています。主弁を開き始める0.2秒が、ステップの境目ではなく途中に来るので、その位置しだいで誤差が変わるためです。オイラー法と2段は、安定の限界にも近づいています。

入力がある時刻でいきなり切り替わる場合は、影響がもっとはっきり出ます。試しに、起動トルクがまだ約140 N·m効いている0.3秒で、いきなり0に切ってみると、2段も4段も1次の精度に落ちました。0.3秒ちょうどで終わるステップで、最後の傾き（2段の \( \boldsymbol{k}_2 \)、4段の \( \boldsymbol{k}_4 \)）だけが、切ったあとのトルクで計算されるためです。商用ソフトの多くは、こうした切り替えの時刻でステップを区切り直して（イベント処理）、精度が落ちないようにしています。

速い現象（燃焼室、約1 ms）と遅い現象（タンクが空になるまで、約150秒）が同居する問題を、硬い（stiff）問題と呼びます。1Dシミュレーションではごく普通に出てきます。対策をまとめます。

迷ったら4段ルンゲ＝クッタ法を使う（onedsimの既定）。1ステップの手間はオイラー法の4倍だが、同じ時間刻みなら安定な範囲が広く、精度もずっと高い
時間刻みは、安定の限界（今回は4段で約2.4 ms）より十分小さくとる。過渡を正確に追いたいときは、いちばん速い時定数の1/5程度が目安。時定数のおおよその値は「容器の中の質量 ÷ 流量」で見積もれる（流入にブレーキが働くと、実際はこれより短くなる）
今回は、起動と停止の過渡を0.1〜0.2 ms、定常値と150秒の燃焼を0.5 msで計算した
結果が振動したり急に崩れたりしたら、まず時間刻みを半分にして計算し直す。結果が変われば数値の問題、変わらなければモデルの式が表している挙動（それが物理なのか、モデル化の都合なのかは別に考える）
商用ソフトの多くは陰解法（少し先の状態を未知数にして方程式を解く方法。BDF法と、それを実装したソルバーのDASSLが代表例）を備えている。大きな時間刻みでも発散しにくく、速い変化が落ち着いたあとは刻みを自動で広げられる

前編で紹介したEcosimPro（ESPSS）も、DASSLなどの陰解法ソルバーを積んでいます。陰解法が必要な理由は、一度自分でオイラー法を発散させてみるとよく分かります（汗）

燃焼振動：数値の振動と本物の振動

今回出た振動は計算のせいだったけど、本物のエンジンでも燃焼室が振動することはあるの？　それも1Dで計算できる？

燃焼室には本物の振動もあります。燃焼振動（燃焼不安定）です。燃焼の発熱と圧力の揺れのタイミングが合うと振幅が育ち、短時間で噴射器や燃焼室の壁を壊すことがあります。

燃焼振動は、どの固有振動数と結びつくかで分類されます（周波数の区分はStanford大学の講義資料による）。燃焼室そのものだけでなく、配管を含めたシステム全体の固有振動数も関わります。

種類	周波数の目安	何の固有振動数か	1Dで扱えるか
高周波（スクリーミング）	1,000 Hz以上	燃焼室そのものの音響（縦方向・接線方向・半径方向）	×　音響解析・3D CFD・燃焼試験
中周波（バズ）	約400〜1,000 Hz	噴射器・供給配管の音響と燃焼の結合	△　1Dと音響解析の組み合わせ
低周波（チャギング）	400 Hz程度まで	燃焼室の容量（上の時定数）、噴射器の差圧、燃焼の遅れ、供給配管やポンプまで含めたシステム全体	○　1Dの得意分野
（参考）POGO	数Hz〜数十Hz	燃焼振動ではなく、推進系の配管と機体構造の連成振動	○　構造モデルと連成

高周波の振動数は、燃焼ガスの音速と燃焼室の寸法で決まります。

燃焼室の音響の固有振動数（目安）

\[ f_{1L}=\frac{a}{2L},\quad f_{1T}=\frac{1.841\,a}{\pi D} \]

\( a \) は燃焼ガスの音速（今回の条件で約1,210 m/s。化学平衡計算の値1,245 m/sに、燃焼効率97.5%を掛けた値）、\( L \) は燃焼室の長さ、\( D \) は直径です。1Lは縦方向1次、1Tは接線方向1次のモードで、1.841はベッセル関数から来る係数です。

今回のモデルの燃焼室について、断面積をスロートの3倍、特性長さを1.0 m、絞り部の半角を30°と仮定すると、直径は約39 cm、噴射面からスロートまでは約38 cmになります。これを当てはめると、1Lモードは約1.6 kHz、1Tモードは約1.8 kHzになります（実寸は非公表なので目安です）。

燃焼室の形そのものが効く現象なので、1Dではなく音響解析や3Dの出番です。実際のエンジン開発では、噴射器の設計や、バッフル（噴射面の仕切り板）・音響キャビティで安定性を確保したうえで、燃焼中の燃焼室に小さな爆薬で圧力の衝撃を与え、揺れが収まるかを確かめます（ボム試験）。

低周波のチャギングはシステム全体の問題です。先ほどのブレーキ（燃焼室の圧力が上がると流入が減る）は、推進剤が燃えるまでの時間のぶんだけ遅れて効きます。この時間遅れを含んだフィードバックが、条件によって発振します。

チャギング（低周波の燃焼振動）が起きるフィードバックの模式図 — チャギングのしくみ。燃焼の遅れをはさんだフィードバックが発振する

そこに、上流の配管を流れる液体の慣性や、ポンプの特性が加わります。配管の慣性は次の式で表せます。

配管の慣性（長さ l、断面積 A）

\[ \begin{aligned} \frac{d\dot m}{dt}&=\frac{A}{l}\bigl(P_{\mathrm{up}}-P_{\mathrm{down}}\\ &\qquad-\Delta P_{\mathrm{loss}}\bigr) \end{aligned} \]

対策の基本は、噴射器の差圧を大きくとることです。燃焼室の圧力の揺れに対して流量が鈍感になり、燃焼室と供給系の結びつきが弱まります。

チャギングが起きるかどうかは、配管を慣性と容量を持つ部品に置き換え、燃焼の遅れを入れたうえで、システム全体を線形化して固有値を調べれば分かります。実部がプラスの固有値が出れば不安定です。燃焼の遅れは「一定の時間だけ遅れて効く」要素なので、そのままでは状態方程式に乗りません。1次遅れ（時定数を持つ追従）をいくつか直列につないで近似します。今回のモデルには配管の慣性も燃焼の遅れも入れていないので、固有値はすべて負の実数で、チャギングは出ません。

チャギングは、公開したコードに配管の慣性と燃焼の遅れを足せば試せます。このあとの「試してほしい改造」に入れておきました。

ツールを使ってみよう

コード一式はGitHubで公開しています（MITライセンス）。

GitHubで onedsim を見る（1D-simulation-demo）

ダウンロードするときは、リポジトリ（SharedTools）のトップページ右上の「Code」→「Download ZIP」で全体を取得し、中の 1D-simulation-demo フォルダを使います。そのフォルダで pip install -r requirements.txt を実行すると、python pumped_storage.py で揚水発電所のモデルが動きます。

uvを使っている場合は、uv run --with-requirements requirements.txt pumped_storage.py の1行で、ライブラリの用意から実行までできます。

ファイル	内容
onedsim.py	ソルバー本体（Table1D / bisect / Model / simulate / sweep）
pumped_storage.py	揚水発電所のモデル。python pumped_storage.py で7日間の要約を表示
gg_rocket.py	ガスジェネレータサイクルエンジンのモデル。python gg_rocket.py で設計計算（教科書の値に置き換えた場合も）と、起動から停止までの結果を表示
cea_tables.py	燃焼ガスの表をRocketCEAで作り直すスクリプト（任意。RocketCEAが必要）
run_studies.py	この記事のグラフと主な数値をまとめて再計算（10〜15分）
test_onedsim.py test_pumped_storage.py test_gg_rocket.py	ソルバー・揚水発電所・ロケットエンジンのモデルの確認（積分法の精度と安定の限界、水やガスの質量保存、エネルギー保存など）。pytestで実行
requirements.txt	必要なライブラリ（numpy / matplotlib）
README.md	使い方と、自分のモデルを書くときの雛形

RocketCEAの使い方は、こちらの記事で紹介しています。

試してほしい改造を挙げておきます。

揚水：n_units=6 にして、神流川の完成形（282万kW）でどこまで需要に応えられるか見る
揚水：固定速のポンプを可変速にして、揚水の入力を余った電力に合わせる制御を入れる
ロケット：gg_valve_t、main_valve_t、spin_torque を変えて、回転数の行き過ぎを出さずに、推力の立ち上がりをどこまで速くできるか探す
ロケット：主弁を開く条件を、時刻ではなく「ポンプの吐出圧が決めた値を超えたら」に書き換えて、起動トルクがばらついたときの行き過ぎがどう変わるか見る
ロケット：cutoff_t、gg_cut_t、main_cut_t、gg_lox_cut_lag を変えて、停止のあとに残る力積と、吐出圧のピークを両方小さくする手順を探す
ロケット：design_point() のタービン効率や入口温度、ポンプの吐出圧、燃焼圧を変えて、比推力とターボポンプの動力がどう動くか確かめる
ロケット：タービンの出口に面積一定の排気ノズルを付けて、出口の圧力が流量で変わるようにし、推力を絞ったときの落ち方がどう変わるか見る
ロケット：gg_n2_fill（起動前のGGに入っている窒素の量）を変えて、弁を開く順番がずれたときの温度のピークがどう変わるか見る
ロケット：供給配管の慣性と燃焼の遅れ（1次遅れを数段つないで近似）を追加して、噴射器の差圧を小さくするとチャギングが出るか試す
自分の分野のシステムをひとつ書いてみる（冷却水回路でも油圧回路でも、部品と保存則は同じ）

最後に

前後編を通して伝えたかったのは、1Dシミュレーションは3Dの簡易版ではなく、部品の特性をつないでシステム全体の挙動と限界を見る、別の道具だということです。たった100行のソルバーでも、考え方自体や基本骨格は商用ソフトと同じです。揚水発電所の供給の限界も、ロケットエンジンの酸燃比や、起動と停止の手順の難しさも、簡易モデルでここまでグラフにできました。

次にP&ID（配管計装図）や系統図を見たとき、「これは derivatives() に書けるな」と思えたら、もう1Dシミュレーションを使う側の人です。是非ご自分のシステムで回してみてください！

前編はこちらです。