はじめに
- 1Dシミュレーションを、自分の手で動かしてみたい
- 商用ソフトを買う前に、1Dツールの仕組みを知りたい
- 揚水発電所やロケットエンジンの過渡挙動を、Pythonで解いてみたい
- 今回の記事のレーティング：
前編で「Pythonで100行もあれば書ける」って言ってたけど、本当？ 商用の1Dツールって、個人では手が出ない値段だよね。。。
みなさんこんにちは！本日は1Dシミュレーションの後編です。かずおさん、本当ですよ。骨格だけなら約100行で書けます。商用ソフトが高価なのは、部品ライブラリや高機能なソルバー、画面まわりまで作り込んであるからなんです。
（前編に引き続き、MATLABを日頃からいじっている人は生暖かい目で見守ってください笑）
前編では、1D CAEとは何か、そして3Dだけでは扱いきれないシステム全体の挙動を、なぜ1Dなら扱えるのかを紹介しました。今回の後編では、前編で予告した揚水発電所とロケットエンジンを、Pythonで1Dシミュレーションのモデルにして動かします。
今回作るもの
作るのは次の3つです。コードはすべてGitHubで公開していて、この記事のグラフもすべてそのコードの計算結果です（リンクは記事の最後にあります）。
- onedsim：numpyだけで動く1Dシミュレーションの骨格。
- 日本最大級の揚水発電所のモデル：東京電力リニューアブルパワーの神流川発電所の公表値を参考に、1週間の運転を解いてみる
- Merlin1Dロケットエンジンのモデル：SpaceXのMerlin 1Dを題材に、ガスジェネレータサイクル（推進剤の一部を燃やしたガスでポンプを回す方式）のエンジンの起動・停止と、燃焼終了までの通しの運転を解いてみる
ロケットエンジンの燃焼圧や効率などの前提は、SpaceX特集 第3回の2ページ目「Raptor 3の強みの検証」とそろえました。
前提知識
- 質量保存：タンクの中身が変わる速さ ＝ 入る流量 − 出る流量
- エネルギー保存：軸の回転エネルギーが変わる速さ ＝ タービンの出力 − ポンプの動力
- 時間積分：「今の状態」と「変化率」から「少し先の状態」を求める
- Python 3 と numpy、matplotlib が動く環境
保存則って難しそうに聞こえるけど、「入る量 − 出る量」なら家計簿と同じだね。それなら追えそう。
この記事で使う用語をまとめておきます。
| 用語 | 意味 |
|---|---|
| 状態量 | 時間とともに変化し、微分方程式で追いかける量（タンクの水量、軸の回転数、燃焼室の圧力など） |
| 性能曲線・マップ | 部品の性能を表すデータ。ポンプなら、流量と揚程（水を押し上げられる高さ）の関係を表すQ-H曲線や、効率曲線。本来は試験で求める。1Dモデルはこれを表として引く |
| 動作点 | 部品の性能曲線と、システムの要求曲線が交わる点。システム全体を連立して初めて決まる |
| 過渡応答 | 起動・停止・外乱のあと、システムが新しい定常状態に落ち着くまでの時間変化 |
| パラメータスタディ | 条件をひとつずつ変えて何度も計算し、感度や限界を調べること |
Pythonで書く1Dソルバー「onedsim」の中身
1Dソルバー（方程式を、時間を追って解くプログラム）の中身は、突き詰めると3つの要素でできています。
① Table1D：試験データの表を引く
前編で、ポンプの性能は最終的に実測で決まると書きました。その実測値を1Dモデルに取り込む要素がTable1Dで、やっていることは線形補間です。
試験データの点 \( (x_i,\,y_i) \) の間を直線で結びます。
\[ \begin{gathered} y(x)=y_i+\frac{y_{i+1}-y_i}{x_{i+1}-x_i}\,(x-x_i)\\ (x_i\le x\le x_{i+1}) \end{gathered} \]範囲の外は外挿せず、端の値で止めます。試験していない範囲の値を作らないためです。
\[ \begin{gathered} y(x\lt x_0)=y_0\\ y(x\gt x_N)=y_N \end{gathered} \]Table1D(x, y)：試験データ（カタログのQ-H曲線など）を登録するtable(xq)：任意の点 xq での値を返す
試験データの表さえあれば、部品の中の流れが分からなくてもモデルに組み込めます。前編で書いた、速い現象を実測から得た特性として畳み込むという話の正体が、この表です。
② Model：状態量とその変化率
システム全体を「状態量」と「その変化率を返す関数」で表します。部品をつなぐとは、この関数の中に質量保存とエネルギー保存を書くことです。
状態量をまとめたベクトルを \( \boldsymbol{s} \) とすると、1Dモデルはこの形の連立常微分方程式になります。
\[ \frac{d\boldsymbol{s}}{dt}=\boldsymbol{f}(t,\,\boldsymbol{s}),\quad \boldsymbol{s}(0)=\boldsymbol{s}_0 \]発電出力や流量のように、状態量から計算できる量は、記録用の量として別に持ちます。
\[ \boldsymbol{y}=\boldsymbol{g}(t,\,\boldsymbol{s}) \]initial_state()：初期値 \( \boldsymbol{s}_0 \) を返すderivatives(t, s)：変化率 \( \boldsymbol{f}(t,\boldsymbol{s}) \) を返す。部品の特性と保存則はすべてここに書くoutputs(t, s)：記録したい量 \( \boldsymbol{g}(t,\boldsymbol{s}) \) を返すclamp(s)：物理的にありえない値を防ぐ（水量はマイナスにならない、など）
「部品をつなぐ」って、画面の上でアイコンを線で結ぶイメージだったよ。
商用ソフトのあの画面も、裏側ではこの関数を自動で組み立てています。今回はそれを手で書く、というだけの違いです。
③ simulate：時間積分
simulate は、前提知識で挙げた時間積分を、時間刻みΔt秒ごとに繰り返す部分です。方法は3つ用意しました。3つある理由は、4ページ目の「時間刻みと計算の発散」で説明します。
今の傾きがΔt秒間続くと仮定して進めます。
\[ \boldsymbol{s}_{n+1}=\boldsymbol{s}_n+\Delta t\,\boldsymbol{f}(t_n,\,\boldsymbol{s}_n) \]このあとも式が続くの。。。？ 3つの違いだけ先に教えて！
違いは「1ステップで傾き（変化率）を何回計算するか」だけです。下の2つの式は読み飛ばしても、この先の話は追えますよ。式のあとに、図でも比べます。
オイラー法で一度仮に進め、その先の傾きも計算して、始点と終点の傾きを平均します。
\[ \begin{aligned} \boldsymbol{k}_1&=\boldsymbol{f}(t_n,\,\boldsymbol{s}_n)\\ \boldsymbol{k}_2&=\boldsymbol{f}(t_n+\Delta t,\,\boldsymbol{s}_n+\Delta t\,\boldsymbol{k}_1)\\ \boldsymbol{s}_{n+1}&=\boldsymbol{s}_n+\frac{\Delta t}{2}\left(\boldsymbol{k}_1+\boldsymbol{k}_2\right) \end{aligned} \]区間の始点で1回、中点で2回、終点で1回、合わせて4回傾きを計算し、重みを付けて平均します。
\[ \begin{aligned} \boldsymbol{k}_1&=\boldsymbol{f}(t_n,\,\boldsymbol{s}_n)\\ \boldsymbol{k}_2&=\boldsymbol{f}\!\left(t_n+\tfrac{\Delta t}{2},\,\boldsymbol{s}_n+\tfrac{\Delta t}{2}\boldsymbol{k}_1\right)\\ \boldsymbol{k}_3&=\boldsymbol{f}\!\left(t_n+\tfrac{\Delta t}{2},\,\boldsymbol{s}_n+\tfrac{\Delta t}{2}\boldsymbol{k}_2\right)\\ \boldsymbol{k}_4&=\boldsymbol{f}(t_n+\Delta t,\,\boldsymbol{s}_n+\Delta t\,\boldsymbol{k}_3)\\ \boldsymbol{s}_{n+1}&=\boldsymbol{s}_n+\frac{\Delta t}{6}\bigl(\boldsymbol{k}_1+2\boldsymbol{k}_2\\ &\qquad+2\boldsymbol{k}_3+\boldsymbol{k}_4\bigr) \end{aligned} \]1ステップの計算量は \( \boldsymbol{f} \) を呼ぶ回数で決まり、オイラー法が1回、2段が2回、4段が4回です。手間が増える代わりに精度は上がり、Δtを半分にしたときの誤差は、それぞれ約1/2、1/4、1/16に減ります（式が滑らかな場合の理論値）。
simulate(model, t_end, dt, method)：時間積分を実行する。method は euler（オイラー法）、rk2（2段）、rk4（4段。既定）record_every：Nステップごとに記録する（長い計算でもメモリを食わない）stop_when：条件を満たしたら打ち切る（例：タンクが空になった）sweep(...)：条件を変えて simulate を繰り返すbisect(f, a, b)：二分法（答えをはさんだ範囲を、半分ずつ狭めていく方法）で、部品どうしのつり合い点（ポンプの動作点など）を探す。範囲が答えをはさんでいなければエラーにする
3つの方法が1ステップをどう進めるかを、図で比べます。例にしたのは、\( ds/dt=s \)（最初の値は1）を時間刻みΔt＝1で1ステップだけ進める計算です。図の濃い色の線が正しい解です。
1ステップ後の正しい値は2.718です。オイラー法は2.000、2段は2.500、4段は2.708で、誤差は0.718、0.218、0.010と小さくなります。
骨格はこれだけです。あとは部品の特性と保存則を derivatives() に書けば、揚水発電所もロケットエンジンも同じ simulate() で回ります。
実践①：揚水発電所
モデルの構成
揚水発電所は、電力が余る夜に下池の水を上池へくみ上げておき（揚水）、電力が足りない昼にその水を落として発電します。
題材は東京電力リニューアブルパワーの神流川発電所です。上池は南相木ダム、下池は上野ダムで、有効貯水量はどちらも1,267万m³です。水位の標高は、上池が満水位1,527 m・最低水位1,500 m、下池が満水位843.4 m・最低水位814.0 mです（電力土木技術協会の水力発電所データベース）。有効落差（上池と下池の水位差から、水路での損失を引いた値）は、基準が653 m、最高が675 m、最低が617 mです（東芝レビュー、2006年）。以下、単に「落差」と書くときは、損失を引く前の上池と下池の水位差を指します。
現在は2台・94万kW（940 MW）が運転中で、最終的には6台・282万kWと世界最大級の揚水発電所になる計画です。今回は2台の構成で、この数字を参考にモデルを組みます。
写真は別の揚水発電所の例です。床から出ているのが発電電動機で、その真下にポンプ水車がつながっています。
モデルに入れたのは、次の6項目です。
- 上池・下池：最低水位より上の水量を状態量とするタンク。面積は一定とし、有効貯水量を、満水位と最低水位の差（上池27.0 m、下池29.4 m）で割った値にした（上池は約47万m²、下池は約43万m²）。水位の標高 ＝ 最低水位 ＋ 水量 ÷ 面積で、落差は約657〜713 mの範囲で変わる（上池も下池も水量が半分のとき684.8 m。このあとの1週間の計算では662〜695 m）
- ポンプ水車：Q-H曲線と効率曲線を表で持つ。実機の曲線は非公表なので、形は仮定した。水車は、公表されている仕様の、有効落差653 m・流量85 m³/sで出力482 MWという点に合わせてある。ポンプの公表仕様は揚水量60 m³/s・最大軸入力464 MWで、モデルのポンプは、両池の水量が半分のとき約57 m³/s
- 発電電動機：効率97.5%。水車出力482 MWのときの発電出力は470 MWで、2台なら公表値の940 MWになる
- 水路：損失 ＝ k×Q²（Qは1台あたりの流量）。損失は公表されていないので、85 m³/sを流したときの有効落差の範囲が、公表されている最高675 m・最低617 mに合うように決めた（損失は約39 m。モデルの有効落差の範囲は約618〜674 m）。有効落差が基準の653 mになるのは、落差が約692 mのとき
- 運用のルール：深夜22〜6時は、電力系統で余った電力（平日900 MW、週末630 MWと想定）でまかなえる台数だけ揚水する（1台あたり約450 MW）。昼13〜16時と夕方18〜21時は900 MWの発電を要求する（週末は5〜6割に減らす）。上池の水量が有効貯水量の5%（水位1.35 m）以下になったら発電を止め、下池が5%以下になるか、上池が満水になったら揚水を止める
- 外乱：降雨（上池に雨が集まる範囲＝集水域は6.2 km²、降った雨のうち池に流れ込む割合＝流出率は0.8）、蒸発4 mm/日、猛暑日の需要増。下池に流れ込む川の水は、同じ量を下流へ流すものとして、モデルには入れていない
\[ \begin{aligned} \frac{dV_u}{dt}&=Q_{\mathrm{pump}}-Q_{\mathrm{turbine}}\\ &\quad+Q_{\mathrm{rain}}-Q_{\mathrm{evap}}\\ &\quad-Q_{\mathrm{spill},u}\\[6pt] \frac{dV_l}{dt}&=Q_{\mathrm{turbine}}-Q_{\mathrm{pump}}\\ &\quad-Q_{\mathrm{evap}}-Q_{\mathrm{spill},l} \end{aligned} \]
\( V_u \)、\( V_l \) は上池と下池の水量、\( Q \) は流量です。\( Q_{\mathrm{spill}} \) は、池が満水のときにあふれる分です（満水でなければ0）。コードでは、発電量などの積算値も同じ形の式で積分しています。
式が2本だけ？ 揚水発電所って、もっと複雑なんじゃないの？
複雑さは「今この瞬間、ポンプ水車がどれだけ流すか」を決める部分に入っています。次の動作点の話です。
ポンプの動作点
ポンプの流量は、ポンプだけでは決まりません。ポンプのQ-H曲線と、システムの要求曲線（落差＋水路の損失）の交点で決まります。落差は上池と下池の水位で変わるので、動作点も時々刻々と動きます。
コードでは、この交点を二分法（bisect）で毎ステップ求めています。商用ソフトの多くは、こうした交点をシステム全体の連立方程式としてまとめて解きます。解き方は違いますが、つり合う点を探しているという中身は同じです。
def pump_operating_point(self, H_static):
"""ポンプ曲線 H_pump(Q) と、システムの要求 H_static + k Q^2 の交点（1台あたり）。"""
f = lambda Q: self.pump_H(Q) - (H_static + self.k_loss * Q ** 2)
if f(0.0) <= 0:
return 0.0
return bisect(f, 0.0, 80.0)カタログに書いてある「揚水量」は、ある落差のときの交点なんだね。
発電のときは逆に、要求された出力になるように流量を決めます（実機ではガイドベーンの開度で調整します）。動かすのは要求をまかなえる最小の台数で、流量には1台あたり85 m³/sの上限があります（落差によらず一定と仮定）。
1週間の運転
月曜から日曜までの7日間を、時間刻み60秒で計算しました。水曜は猛暑日（昼間ずっと発電し、ピークには定格いっぱいの940 MWを要求）、土曜は12時間の豪雨（30 mm/h）という設定です。計算は1秒ほどで終わります。
平日は夜に揚水、昼と夕方に発電を繰り返します。平日の揚水の入力は約860〜940 MWの間で動き、余った電力の900 MWにぴったり合わせることはできません。固定速のポンプ水車は、落差が決まれば流量も入力も決まってしまうためです。
この弱点を補うのが、回転数を変えて揚水の入力を調整できる「可変速機」です。発電電動機の回転子に、周波数変換器（インバータやサイクロコンバータ）で作った低い周波数の交流を流して、回転数を数%の幅で変えます。可変速機は1990年に実用化されています（東芝レビュー、2025年）。「モデルの構成」で写真を載せた奥清津第二発電所でも、2号機が可変速機です。神流川発電所の1・2号機は回転数が一定の同期機なので、今回のモデルも固定速にしています。
グラフの水位は、最低水位（標高1,500 m）から測った高さです。猛暑日は11時間発電し続け、上池の水位は17.4 mから3.3 mまで下がりました。運用の下限とした1.35 m（有効貯水量の5%）までは下がりませんが、落差は664 mまで小さくなりました。
このモデルで2台が940 MWを出せるのは、落差が約692 m以上のときです。猛暑日に940 MWを要求した時間帯は、落差がそれより小さく、出力は最大でも約930 MW、最後には約900 MWまで下がりました。940 MWの要求に対して、約150 MWhの供給不足が出ています。
水が足りなくなるのは想像できるけど、「落差が下がると出力が出ない」は計算しないと気づかないね。
1週間の往復効率（同じ量の水を上げ下げしたときの、発電量 ÷ 揚水に使った電力量）は約70%でした。実際の揚水発電のエネルギーロスは約3割（資源エネルギー庁）なので、ほぼ同じ水準です。ただし、今回のモデルには、変圧器や補機の損失は入っていません。なお、これらの値は、水路の損失に大きく左右されます。損失を10 m小さくすると、猛暑日の不足は約60 MWh、往復効率は約72%になります。10 m大きくすると、不足は約360 MWh、往復効率は約69%です。
条件を変えて限界を探す
1条件を変えて、どこまで耐えられるかを調べます。まずは豪雨です。猛暑日のない通常の週で、土曜の朝から24時間降り続くとして、強さを0〜100 mm/hで変えました。夜の揚水を予定どおり続ける場合と、雨のあいだは揚水しない場合の2通りを計算しています。
上池に収まる限界は、予定どおり揚水を続けると約57 mm/h、雨のあいだ揚水を止めると約70 mm/hで、これより強い雨では上池からあふれます（越流）。
実際のダムには洪水吐（こうずいばき）があるので、あふれても水は安全に流せます。豪雨が貯水池に貯まった分は一切の揚水仕事をしなくて済むので、実質的には無料で電気が手に入ることになります。
天気予報で運転パターンを決めるのも電力会社の仕事の一つと言うことですね！
次に、最初の1週間（水曜が猛暑日、土曜が豪雨）に戻して、ピーク時の発電要求を600〜1,200 MWで変えました。猛暑日は11時間ずっと、ほかの平日は昼と夕方に3時間ずつ、この値を要求します（週末はその5〜6割）。
不足が出始めるのは、要求が約900 MWを超えたあたりからです。それより上では、上池の水位が下がって落差が小さくなった時間帯から、出力が足りなくなります。2台の定格の940 MWを要求すると、不足は約0.75 GWh（約750 MWh）です。最初の1週間（約150 MWh）より多いのは、猛暑日だけでなく平日のピークすべてで940 MWを要求するためです。両池の水量が半分のとき（落差684.8 m）に2台で出せるのは約930 MWで、940 MWには届きません。定格を超えた分はそのまま不足に上乗せされ、1,000 MWでは約2.8 GWhです。
一方、1週間のうちで上池の水位がいちばん下がったときの値は、要求が940 MW以上では1.8 mのまま変わりません。すでに設備の上限で運転していて、要求を上げても使う水の量が変わらないためです。
「1週間の上池の水位」のグラフをもう一度見てください。猛暑日に使った水は、一晩8時間の揚水では戻しきれません。木曜と金曜は、水位が低いまま運転することになります。1週間を通して解いてみると、こういうところが見えてきます。
2022年3月の電力需給ひっ迫では、東京電力管内の揚水発電の「発電可能量」（くみ上げた水の残量）が、需要のピークを越えるまでもつかどうかが焦点のひとつになりました（資源エネルギー庁）。
発電所の細かい特性や毎日の運用までは、外からは分からないよね。それでも、こういう話を自分で確かめられるの？
確かめられます。公表されている落差・流量・貯水量と、気象庁の降水量のデータがあれば、このくらいの粗さの計算はできます。ニュースで「ダムの貯水量が減っている」「真夏の電力需要がピーク」といった話題が出たときも、あと何時間発電できるのか、雨でどれだけ水が戻るのかを、自分の手で見積もれますよ！
ちなみに2.8GWhと言うと、蓄電池容量に換算すると数百億円規模、コストの高い日本の場合1000億円以上の規模の蓄電設備が必要になるのおで、揚水発電がいかに多くの電力を貯蔵しているのかが良く分かりますね。
※このページの値は、神流川発電所の公表値を参考にした簡易モデルの計算結果です。曲線の形や運用のルールは仮定で、実際の発電所の性能や運用を表すものではありません。
次ページへ
次のページでは、同じソルバーでロケットエンジンのモデルを解いてみます。時間のスケールは、日からミリ秒へ変わります。
