フェーズ5完全 DFN の実装と検証 — PyBaMM で答え合わせ

いよいよ本丸である。フェーズ1で導出した DFN(Doyle–Fuller–Newman)モデルの5式を、 一切の簡略化なしに連立して解く。フェーズ4で学んだ有限体積法(FVM)・DAE・Newton 法という道具は、 すべてこの章のためにあった。自作コードで完全 DFN を解き、 実績あるオープンソースライブラリ PyBaMM と突き合わせて「答え合わせ」を行い、 最後に SPM / SPMe / DFN のモデル階層を同じ土俵で比較して「どのモデルをいつ使うべきか」に決着をつける。

5.1 全体系の再掲と未知数の数え上げ

まず、これから解く方程式系の全体像を一枚にまとめる。フェーズ1で1本ずつ導出した5式を、 「どの変数を」「どの領域で」「どんな型の方程式で」決めるのか、という視点で整理し直す。 実装に取りかかる前にこの整理をやっておくと、コードのどの行がどの式に対応するのかが常に明確になり、 デバッグの効率がまったく違ってくる。

5.1.1 未知場の一覧

DFN モデルの未知場(unknown fields)は次の5つである。座標系はフェーズ1と同じく、 $x$ をセル厚み方向(負極集電体が $x=0$、正極集電体が $x=L=L_n+L_s+L_p$)、 $r$ を粒子中心からの半径とする。符号規約も同じで、放電電流を正($I>0$)、 酸化方向(粒子から Li が出る方向)の界面電流を正($j>0$)とする。 放電時は負極で $j>0$、正極で $j<0$ である。

「微分変数」と「代数変数」の区別はフェーズ4で学んだ DAE(微分代数方程式)の言葉である。 $c_s$ と $c_e$ は保存則から来る放物型 PDE に従い、時間微分 $\partial/\partial t$ をもつ。 一方、$\phi_s$、$\phi_e$、$j$ には時間微分が現れない。電位は「電荷の再配置が拡散よりも桁違いに速い」 という近似(電気的中性)の帰結として、各時刻で瞬時に釣り合いの取れた値を取る。 つまり電位分布と反応分布は、その瞬間の濃度分布に対して拘束条件として決まる。 これが DFN が ODE 系ではなく DAE 系になる理由である。

5.1.2 支配方程式と境界条件の一覧表

5本の支配方程式を境界条件とセットで再掲する(導出はフェーズ1)。 以下、実効物性は Bruggeman 補正 $D_e^{\mathrm{eff}} = D_e \varepsilon_e^{b}$、 $\kappa^{\mathrm{eff}} = \kappa \varepsilon_e^{b}$、$\sigma^{\mathrm{eff}} = \sigma \varepsilon_s^{b}$($b=1.5$)、 比界面積は $a = 3\varepsilon_s/R_s$ [1/m] である。

どの $x$ でも固相電流と液相電流の和は印加電流に等しい:

これは式 3 と式 4 の発散の式($\partial i_s/\partial x = -aj$ と $\partial i_e/\partial x = +aj$、 符号は電子と Li⁺ で電荷の運び手が入れ替わるため逆)を足すと $\partial(i_s+i_e)/\partial x = 0$ となり、集電体の境界条件($x=0$ で $i_s=I,\ i_e=0$)から 定数が $I$ に決まる、という論法で出てくる。実装後の検算として非常に有用なので覚えておこう (各 $x$ で $i_s + i_e - I$ を評価してゼロになるか確認できる)。

5.1.3 方程式と未知数の数が合うことの確認

連立系を組む前に、「式の数=未知数の数」を確認する。まず連続の(離散化前の)レベルで数える。 電極領域の1点 $x$ に注目すると:

セパレータ領域には固相が存在しないので、未知数は $c_e$ と $\phi_e$ の2つ、 式は $aj=0$ とした式2と式4の2本。ここでも 2 ↔ 2 で対応する。 したがって形式的には式の数と未知数の数はぴったり一致する。 ただし1つだけ落とし穴がある。それが次の「電位の基準(ゲージ)」の問題である。

5.1.4 電位の基準(ゲージ)と端子電圧

式3〜5をよく見ると、電位は $\eta = \phi_s - \phi_e - U$ と勾配 $\partial\phi/\partial x$ の形でしか現れない。 つまり $\phi_s$ と $\phi_e$ の両方に同じ定数 $C$ を足しても、方程式系はまったく変化しない ($\eta$ は差なので不変、勾配は定数の微分なので不変)。これは電磁気学でおなじみの 「電位の原点は自由に選べる」というゲージ自由度(gauge freedom)がそのまま残っていることを意味する。 数値的には、この自由度のせいで連立方程式のヤコビアンが特異(singular)になり、 Newton 法の線形ソルバが破綻する。だから基準を1つ、人為的に固定する必要がある。

本教材では負極集電体側の固相電位を基準に取る。すなわち負極集電体に接するノードの $\phi_s$ を 0 に固定するゲージ条件を1本加える(離散化での具体的な扱いは 5.3.6 節)。 ここで大事な注意がある。基準を固定しても、物理的な観測量である端子電圧はあくまで「差」で定義する:

「$\phi_s(0)=0$ と固定したのだから $V = \phi_s(L)$ でよいのでは?」と思うかもしれないが、そうしない。 理由は2つある。第一に、有限体積法ではノード(セル中心)は $x=0$ ちょうどには置かれず、 集電体表面 $x=0$ の電位はノード値から境界流束($i_s = I$)を使って外挿した値になる。 ゲージで固定するのは「最初のノードの値」であって「$x=0$ の値」ではないから、 一般に $\phi_s(0) \ne 0$ である(その差は集電体からノード中心までの半セル分のオーム降下)。 第二に、差で定義しておけばゲージの選び方に依存しない量になり、 ゲージ条件を書き換えても(たとえば別のノードを基準にしても)$V$ は変わらない。 観測量はゲージ不変量で書く — これは電磁気学と同じ作法である。

5.2 SPMe からの差分 — 何をやめて何を足すか

フェーズ3で実装した SPMe は「電極ごとに代表粒子1個+電解液の1次元拡散」というモデルだった。 DFN への拡張は、ゼロから別物を作るのではなく、SPMe で置いていた仮定を1つずつ外していく 作業として理解するのが見通しがよい。外す仮定は本質的に1つ: 「反応は電極内で一様に起こる」である。

SPMe では界面電流密度 $j$ を電極厚み方向に一様と仮定し、電流バランス $\int_0^{L_n} a\, j\, dx = I$ から $j$ を与えられた定数 $j_n = I/(a_n L_n)$(負極)、$j_p = -I/(a_p L_p)$(正極)としていた。 $j$ が既知なら、粒子拡散(式1)の境界流束も、電解液(式2)のソース項もすべて既知になり、 5式のうち式1と式2だけを解けばよかった。式3・4・5は「解く」対象ではなく、 電圧を計算するときの後処理(解析的な補正項)に格下げされていた。

DFN ではこの仮定をやめる。$j(x)$ は各 $x$ で Butler–Volmer 式(式5)を満たすように解かれる未知数になる。 すると BV 式に入る $\phi_s$ と $\phi_e$ も場として必要になるから、式3と式4も連立系に昇格する。 つまり「$j$ の一様性を外す」という1つの決断が、芋づる式に5式全部の連立を要求するのである。

この差分がもたらす計算コストの増加は劇的である。SPMe の未知数は粒子2個分+電解液で高々百個程度だったが、 DFN では「$x$ ノードの数 × 粒子シェルの数」が支配的になり、後述のグリッドでは約 400 個になる。 さらに ODE ではなく DAE なので、フェーズ4の 4.5〜4.6 節で学んだ 「質量行列をもつ残差形式+Newton 法」が必須になる。 その見返りとして、SPMe では原理的に見えなかった電極内の不均一性 ($c_e$ の枯渇、反応の局在化、表面充填率のばらつき)が見えるようになる。これは 5.6 節で目撃する。

5.3 離散化の設計 — グリッド・変数ベクトル・DAE 構造

数値解法の道具立てはすべてフェーズ4で学んだものを使う:空間は有限体積法(FVM)、 時間は後退オイラー、非線形連立は Newton 法である。この節では「どの変数をどう並べ、 どの残差をどう書くか」という設計図を固める。5.4 節のコードはこの設計図の忠実な翻訳である。

5.3.1 空間グリッド

$x$ 方向は3領域をそれぞれ一様分割する:負極 $n_n = 12$ ノード、セパレータ $n_s = 6$ ノード、 正極 $n_p = 12$ ノード、計 $N_x = 30$ ノード。セル幅は領域ごとに異なり、 $\Delta x_n = L_n/n_n$、$\Delta x_s = L_s/n_s$、$\Delta x_p = L_p/n_p$ である。 ノードは各セルの中心に置く(セル中心 FVM)。 $r$ 方向は各電極ノードの粒子を $n_r = 12$ シェルに一様分割する($\Delta r = R_s/n_r$)。 フェーズ2で作った球対称 FVM(シェル体積 $\frac{4}{3}\pi(r_{i+1/2}^3 - r_{i-1/2}^3)$ で重み付け)を そのまま各粒子に適用する。

グリッドを粗めに取っているのは教材としての実行時間(1C 放電で数十秒)を優先したためである。 収束性の確認(グリッドを倍にして解が変わらないこと)はフェーズ4の 4.7 節の作法どおり、 各自で必ずやってほしい。この解像度でも電圧誤差は PyBaMM 比で 10 mV 程度に収まる。

5.3.2 変数ベクトル $\mathbf{y}$ のレイアウト

すべての未知数を1本のベクトル $\mathbf{y}$ に並べる。並べ方は自由だが、 一度決めたら全コードで厳守する。本教材の並び順は次のとおり:

粒子ブロック $c_s^n$ の内部は「$x$ ノードごとに $r$ シェルを連続に置く」二重ループ順、 すなわち成分 $(i,\,m)$($i$:$x$ ノード番号、$m$:シェル番号)が添字 $i \cdot n_r + m$ に入る。 各ブロックの大きさと通し添字を図にすると:

ブロックサイズを合計すると $144 + 144 + 30 + 12 + 12 + 30 + 12 + 12 = 396$。 そのうち微分変数は $c_s^n,\ c_s^p,\ c_e$ の $144+144+30 = 318$ 成分、 代数変数は $\phi_s^n,\ \phi_s^p,\ \phi_e,\ j^n,\ j^p$ の $12+12+30+12+12 = 78$ 成分である。 5.1.3 節の連続レベルの数え上げが、離散化後もそのまま 「318 本の ODE 残差+78 本の代数残差 = 396 本」として保存されていることを確認してほしい (ゲージによる 1 本の入れ替えは 5.3.6 節)。

5.3.3 $x$ 方向 FVM と領域界面の調和平均

$c_e$ の FVM 離散化はフェーズ4の 4.2 節と同じく、セル $i$ について

と書ける($F_{i\pm1/2}$ はセル界面の拡散流束)。問題は界面の実効拡散係数 $D^{\mathrm{eff}}_{i+1/2}$ をどう取るかである。領域の内部では両隣のセルの物性が同じなので 単純平均でも調和平均でも大差ない。しかし負極/セパレータ、セパレータ/正極の境界をまたぐ界面では、 空隙率 $\varepsilon_e$ が不連続に変わる(0.25 → 0.47 → 0.335)ため $D_e^{\mathrm{eff}} = D_e \varepsilon_e^{1.5}$ が両側で数倍違う。 ここで単純算術平均を使うと流束を過大評価する。フェーズ4の 4.2.4 節で導いたとおり、 「2つの抵抗の直列」と見なして距離重み付き調和平均(harmonic mean)を使うのが正しい:

導出の要点だけ再掲する:界面をはさむ2つの半セル(幅 $\Delta x_i/2$ と $\Delta x_{i+1}/2$)を 流束が連続に通過すると要請し、界面値を消去すると、直列合成抵抗 $\frac{\Delta x_i/2}{D_i} + \frac{\Delta x_{i+1}/2}{D_{i+1}}$ が現れる。 これを1つの等価拡散係数に読み替えたものが上式である。 片側の $D$ が極端に小さいとき流束がそちらに律速される、という物理を正しく再現する。 $\kappa^{\mathrm{eff}}$(液相電荷保存の式4)にもまったく同じ扱いを適用する。 なお $D_e(c_e)$ と $\kappa(c_e)$ は濃度依存(params.js の De_fun, kappa_fun)なので、 各セルの $D^{\mathrm{eff}}_i$ は毎ステップそのセルの $c_{e,i}$ で評価してから調和平均を取る。

5.3.4 集電体境界とセパレータ境界の扱い

• $c_e$:両集電体でゼロ流束。FVM では境界界面の流束を単に 0 と置くだけでよい(ゴーストセル不要)。
• $\phi_e$:両集電体で $i_e = 0$。これも境界界面の電流を 0 と置く。
• $\phi_s$:集電体で $i_s = I$。負極の $x=0$ 側界面と正極の $x=L$ 側界面で、 界面電流をちょうど $I$ と置く(放電 $I>0$ のとき、電子の流れとして整合する符号はフェーズ1の図1.4 を参照)。
• $\phi_s$:セパレータ境界で $i_s = 0$。固相はセパレータに続いていないので、 負極ブロックの右端・正極ブロックの左端では固相電流がゼロ。 つまり $\phi_s^n$ と $\phi_s^p$ は互いに独立な2つの小さい楕円型問題になり、 両者をつなぐのは $\phi_e$ と $j$ を経由した間接的な結合だけである。
• 粒子表面:$r=R_s$ の界面流束を $j_i/F$ と置く(その $x$ ノードの $j_i$ を使う)。 粒子中心 $r=0$ はゼロ流束。フェーズ2と同一。

境界条件はすべて「界面の流束を指定する」Neumann 型なので、セル中心 FVM とは相性がよく、 離散式の該当する流束項を置き換えるだけで実装できる。Dirichlet 型の処理(ゴーストセルや外挿)は ゲージ条件(次項)を除いて一切登場しない。

5.3.5 DAE の質量行列構造

離散化後の系は、フェーズ4の 4.5 節で導入した標準形

に収まる。質量行列 $M$ は対角で、微分変数の行に正の係数、代数変数の行に 0 が入る:

($c_e$ の行に $\varepsilon_e$ が付くのは式2の左辺が $\varepsilon_e\, \partial c_e/\partial t$ だから。 実装では両辺を $\varepsilon_e$ で割って 1 に正規化してもよい。) 後退オイラーで時間を離散化すると、時刻 $t_{k+1}$ の解 $\mathbf{y}_{k+1}$ は非線形連立方程式 $\mathbf{F}(\mathbf{y}_{k+1}) = 0$ の根として決まる。残差 $\mathbf{F}$ の作り方が DAE のポイントで、 行ごとに型が違う:

代数行には $\Delta t$ が現れないことに注意。電位と反応分布は「前の時刻の値」を参照せず、 その時刻の濃度場と整合するように毎回釣り合い直される。 初期条件にも同じ注意が要る:$t=0$ で勝手な $\phi_s, \phi_e, j$ を与えることはできず、 初期濃度場と無矛盾な(consistent な)値を代数方程式を解いて求める必要がある (フェーズ4の 4.5.3 節「無矛盾初期化」)。実装では、初期濃度を固定したまま 代数変数だけの小さい Newton 問題を解く。よい初期推定は SPMe の一様解 ($j = \pm I/(aL)$、$\eta$ は BV の逆関数、$\phi_e = 0$ 一様)である。

5.3.6 ゲージ条件の離散的な実装

5.1.4 節で述べたゲージ自由度は、離散系では次の形で現れる: $\phi_s^n,\ \phi_s^p,\ \phi_e$ の FVM 残差(計 $12+12+30=54$ 本)は Neumann 境界条件しか持たないため、 各ブロックの残差の総和を取ると境界流束と $\sum a j \Delta x$ だけが残る。 3つのブロック和のあいだには恒等関係(負極の和 $+$ 正極の和 $=$ 液相の和、いずれも電荷保存で結ばれる) が1本成立し、残差 54 本のうち独立なのは 53 本しかない。 未知数側もゲージ自由度で1つ余っているから、対処は「従属な残差1本をゲージ条件で置き換える」ことである。 具体的には、負極集電体に接する最初のノードの $\phi_s$ 残差を

に置き換える。これで方程式 396 本・未知数 396 個、かつヤコビアンは正則になる。 端子電圧は集電体表面まで外挿した値の差

で計算する(境界セル中心から集電体表面までの半セル分のオーム降下を $i_s = I$ で外挿。 符号は $i_s = -\sigma^{\mathrm{eff}}\partial\phi_s/\partial x$ と $i_s = I > 0$ から: 負極では $x$ が増える向きに $\phi_s$ が下がるので、表面 $x=0$ はノードより高い…と早合点しそうになるが、 落ち着いて微分の向きに沿って確認してほしい。放電では電子は負極集電体へ流れ出るから $x=0$ に向かって $\phi_s$ は下がる。この符号確認も理解度チェックで問う)。

5.4 教材コードの設計と完全リスティング

設計図(5.3 節)をコードに翻訳する。使うのは Python + numpy + scipy(疎行列モジュール scipy.sparse)だけで、DAE ソルバのような既製品には頼らない。 フェーズ4で自作した部品(FVM、後退オイラー、Newton 法)を、そのままスケールアップして組み上げる。 「ライブラリのブラックボックスを一度も通らずに DFN が解けた」という経験こそ、この教材の到達点である。

5.4.1 コードの全体構成

コードは次の部品からなる。どれも 5.3 節の設計項目と1対1に対応する。

5.4.2 ヤコビアンをどう作るか — 有限差分と列のグループ化

Newton 法にはヤコビアン $J_{ik} = \partial F_i/\partial y_k$ が要る。フェーズ4では小さい系だったので 解析的に書き下したが、396×396 の DFN で全成分を手で微分するのは現実的でない (BV 式は $c_e,\ c_{s,\mathrm{surf}},\ \phi_s,\ \phi_e,\ j$ の5変数に依存し、 $D_e^{\mathrm{eff}}(c_e)$ の濃度依存まで含めると項の数が爆発する)。 教材コードでは有限差分近似を使う:第 $k$ 成分だけを $\varepsilon$ だけ動かして

とすれば第 $k$ 列が得られる。素朴にやると残差評価が $N+1 = 397$ 回必要で遅い。 そこで疎性パターンを利用する。$F_i$ が $y_k$ に依存するのは、 離散式が局所結合(隣のセル・隣のシェル・同じ $x$ ノードの変数)だからで、 ヤコビアンの非ゼロ位置はあらかじめ列挙できる。 非ゼロ行を共有しない列同士は、同時に摂動しても互いの差分を汚さない。 そこで列をそのような「同色グループ」に貪欲法で分けておき、1グループ=1回の残差評価で 複数の列を同時に取得する(グラフ彩色によるヤコビアン圧縮。CPR 法とも呼ばれる)。 DFN の結合構造ならグループ数は数十で済み、ヤコビアン1回の構築コストは 残差評価数十回分にまで落ちる。得られた $J$ は scipy.sparse の CSC 形式に詰めて spsolve で $J\,\delta = -\mathbf{F}$ を直接解く(396 元の疎行列なら1回数ミリ秒)。

5.4.3 完全コードリスティング

以下が完全 DFN の教材実装である(検証済み。5.5 節の方法で PyBaMM と突き合わせてある)。 長く見えるが、半分以上は 5.3 節の設計をそのまま式で書き写した残差関数である。 読むときは「変数レイアウト(5.3.2)→ 残差関数 → 時間ループ」の順に追うとよい。

"""
ゼロから書く DFN(Doyle–Fuller–Newman)モデル
=================================================
必要なもの: numpy, scipy のみ。パラメータは LG M50(Chen et al. 2020)。

解法(フェーズ4で学んだ道具の組み合わせ):
  - 空間: 有限体積法(x 方向 FVM + 各電極ノードに粒子の r 方向 FVM)
  - 時間: 後退オイラー(BDF1)
  - 各ステップ: 全変数まとめて Newton 法で残差 F(y)=0 を解く
  - ヤコビアン: 疎構造を明示的に列挙し、彩色(coloring)付き有限差分で構築

変数ベクトル y のレイアウト:
  [ c_s^n (n_n×n_r) | c_s^p (n_p×n_r) | c_e (N) |
    φ_s^n (n_n) | φ_s^p (n_p) | φ_e (N) | j^n (n_n) | j^p (n_p) ]
微分変数: c_s, c_e(後退オイラー残差)/ 代数変数: φ_s, φ_e, j(拘束残差)
"""
import numpy as np
from scipy import sparse
from scipy.sparse.linalg import spsolve

# ============ 1. パラメータ(Chen et al. 2020, LG M50) ============
F, Rg, T = 96485.33212, 8.314462618, 298.15
L_n, L_s, L_p = 85.2e-6, 12.0e-6, 75.6e-6
R_n, R_p = 5.86e-6, 5.22e-6            # 粒子半径 [m]
eps_e_n, eps_e_s, eps_e_p = 0.25, 0.47, 0.335   # 空隙率
eps_s_n, eps_s_p = 0.75, 0.665          # 活物質体積分率
brug = 1.5                              # Bruggeman 指数(電解液側)
Ds_n, Ds_p = 3.3e-14, 4.0e-15           # 固相拡散係数 [m^2/s]
sig_n, sig_p = 215.0, 0.18              # 固相伝導率 [S/m](Chen2020 は実効値扱い)
cs_max_n, cs_max_p = 33133.0, 63104.0
cs0_n, cs0_p = 29866.0, 17038.0         # 初期濃度(100% SOC)
ce0, t_plus = 1000.0, 0.2594
k_n, k_p = 6.716e-12, 3.545e-11         # 反応速度定数 [m^2.5 mol^-0.5 s^-1]
A_cell, cap_Ah = 0.1027, 5.0
I_1C = cap_Ah / A_cell                  # 1C の電流密度 [A/m^2]
V_CUT = 2.5                             # 放電終止電圧 [V]

def U_n(t):   # 黒鉛 OCP [V](Chen2020 フィット)
    return (1.9793*np.exp(-39.3631*t) + 0.2482
            - 0.0909*np.tanh(29.8538*(t-0.1234))
            - 0.04478*np.tanh(14.9159*(t-0.2769))
            - 0.0205*np.tanh(30.4444*(t-0.6103)))

def U_p(t):   # NMC811 OCP [V]
    return (-0.8090*t + 4.4875 - 0.0428*np.tanh(18.5138*(t-0.5542))
            - 17.7326*np.tanh(15.7890*(t-0.3117))
            + 17.5842*np.tanh(15.9308*(t-0.3120)))

def De_fun(ce):   # 電解液拡散係数 [m^2/s](Nyman 2008)
    c = ce/1000.0
    return 8.794e-11*c**2 - 3.972e-10*c + 4.862e-10

def kappa_fun(ce):  # イオン伝導率 [S/m](Nyman 2008)
    c = ce/1000.0
    return 0.1297*c**3 - 2.51*c**1.5 + 3.329*c

# ============ 2. グリッド ============
n_n, n_s, n_p, n_r = 12, 6, 12, 12
N = n_n + n_s + n_p                     # x 方向の全ノード数
dx = np.r_[np.full(n_n, L_n/n_n), np.full(n_s, L_s/n_s), np.full(n_p, L_p/n_p)]
eps_e = np.r_[np.full(n_n, eps_e_n), np.full(n_s, eps_e_s), np.full(n_p, eps_e_p)]
a_n, a_p = 3*eps_s_n/R_n, 3*eps_s_p/R_p           # 比界面積 [1/m]
a_x = np.r_[np.full(n_n, a_n), np.zeros(n_s), np.full(n_p, a_p)]
el_n = np.arange(0, n_n)                 # 負極ノードの x インデックス
el_p = np.arange(n_n+n_s, N)             # 正極ノード

def particle_grid(Rs):
    """球の FVM: 殻中心・体積・面の面積・殻幅"""
    rf = np.linspace(0, Rs, n_r+1)
    Vi = 4/3*np.pi*(rf[1:]**3 - rf[:-1]**3)
    Af = 4*np.pi*rf**2
    return Vi, Af, Rs/n_r

Vi_n, Af_n, dr_n = particle_grid(R_n)
Vi_p, Af_p, dr_p = particle_grid(R_p)

# ---- 変数ベクトルのオフセット ----
o_csn, o_csp = 0, n_n*n_r
o_ce   = o_csp + n_p*n_r
o_fsn  = o_ce + N
o_fsp  = o_fsn + n_n
o_fe   = o_fsp + n_p
o_jn   = o_fe + N
o_jp   = o_jn + n_n
NV     = o_jp + n_p                      # 未知数の総数

def unpack(y):
    return (y[o_csn:o_csp].reshape(n_n, n_r), y[o_csp:o_ce].reshape(n_p, n_r),
            y[o_ce:o_fsn], y[o_fsn:o_fsp], y[o_fsp:o_fe],
            y[o_fe:o_jn], y[o_jn:o_jp], y[o_jp:NV])

# ============ 3. 残差 F(y) = 0 ============
def cs_surf_of(cs, j, Ds, dr):
    """表面濃度: 最外殻の値を境界流束で線形外挿(-Ds dc/dr = j/F)"""
    return cs[:, -1] - (dr/2) * j / (F*Ds)

def particle_res(cs, cs_old, j, Ds, Vi, Af, dr, dt):
    """球拡散(式E1)の後退オイラー残差(電極の全ノード分をまとめて)"""
    fl = np.zeros((cs.shape[0], n_r+1))
    fl[:, 1:-1] = -Ds * (cs[:, 1:] - cs[:, :-1]) / dr   # 内部面のフラックス
    fl[:, -1] = j[:, None][:, 0] / F                    # 表面: 反応流束
    dcdt = (Af[:-1]*fl[:, :-1] - Af[1:]*fl[:, 1:]) / Vi
    return (cs - cs_old)/dt - dcdt

def residual(y, y_old, dt, I_app):
    csn, csp, ce, fsn, fsp, fe, jn, jp = unpack(y)
    csn_o, csp_o, ce_o = unpack(y_old)[:3]
    res = np.empty(NV)

    # ---- 式1: 固相拡散(微分変数) ----
    res[o_csn:o_csp] = particle_res(csn, csn_o, jn, Ds_n, Vi_n, Af_n, dr_n, dt).ravel()
    res[o_csp:o_ce]  = particle_res(csp, csp_o, jp, Ds_p, Vi_p, Af_p, dr_p, dt).ravel()

    # ---- 式2: 液相物質保存(微分変数) ----
    j_x = np.zeros(N); j_x[el_n] = jn; j_x[el_p] = jp
    De_c = De_fun(ce) * eps_e**brug                     # セルの実効拡散係数
    # 面のフラックス(2点式、界面は調和平均に相当する抵抗和で)
    Rf = dx[:-1]/(2*De_c[:-1]) + dx[1:]/(2*De_c[1:])
    Nf = np.zeros(N+1)
    Nf[1:-1] = -(ce[1:] - ce[:-1]) / Rf                 # 両端は流束ゼロ
    src = (1 - t_plus) * a_x * j_x / F
    res[o_ce:o_fsn] = eps_e*(ce - ce_o)/dt - (Nf[:-1] - Nf[1:])/dx - src

    # ---- 式3: 固相電荷保存(代数変数)  ∂i_s/∂x = -a j ----
    for (fs, jj, sig, idx, first_is, last_is) in (
            (fsn, jn, sig_n, el_n, I_app, 0.0),         # 負極: x=0 で I、セパレータ側 0
            (fsp, jp, sig_p, el_p, 0.0, I_app)):        # 正極: セパレータ側 0、x=L で I
        dxe = dx[idx]
        isf = np.zeros(len(idx)+1)
        isf[0], isf[-1] = first_is, last_is
        isf[1:-1] = -sig * (fs[1:] - fs[:-1]) / (0.5*(dxe[:-1]+dxe[1:]))
        r = (isf[1:] - isf[:-1]) + a_x[idx]*jj*dxe      # (i_s の増分) + a j dx = 0
        if idx is el_n:
            res[o_fsn:o_fsp] = r
        else:
            res[o_fsp:o_fe] = r

    # ---- 式4: 液相電荷保存(代数変数)  ∂i_e/∂x = +a j ----
    kap = kappa_fun(ce) * eps_e**brug
    Rk = dx[:-1]/(2*kap[:-1]) + dx[1:]/(2*kap[1:])
    dif = (2*Rg*T/F)*(1 - t_plus)                       # 拡散電位の係数(熱力学因子=1)
    ief = np.zeros(N+1)                                 # 両端は i_e = 0
    ief[1:-1] = (-(fe[1:] - fe[:-1]) + dif*(np.log(ce[1:]) - np.log(ce[:-1]))) / Rk
    r_fe = (ief[1:] - ief[:-1]) - a_x*j_x*dx
    r_fe[0] = fsn[0] - 0.0    # ゲージ固定: φ_s(負極集電体側ノード) = 0
                              # (φ の一様シフト自由度を消す。落とした式は他の式の和から自動で成立)
    res[o_fe:o_jn] = r_fe

    # ---- 式5: Butler–Volmer(代数変数) ----
    for (jj, cs, Ds, dr, k, cs_max, Ufun, idx, off) in (
            (jn, csn, Ds_n, dr_n, k_n, cs_max_n, U_n, el_n, o_jn),
            (jp, csp, Ds_p, dr_p, k_p, cs_max_p, U_p, el_p, o_jp)):
        css = np.clip(cs_surf_of(cs, jj, Ds, dr), 1.0, cs_max-1.0)
        ce_loc = np.maximum(ce[idx], 1.0)
        i0 = F*k*np.sqrt(ce_loc*css*(cs_max - css))
        fs = fsn if idx is el_n else fsp
        eta = fs - fe[idx] - Ufun(css/cs_max)
        res[off:off+len(idx)] = jj - 2*i0*np.sinh(0.5*F*eta/(Rg*T))
    return res

# ============ 4. ヤコビアンの疎構造(手で列挙)と彩色 ============
def build_pattern():
    rows, cols = [], []
    def add(r, c): rows.append(r); cols.append(c)
    # 式1: 粒子内三重対角 + 表面殻は j に依存
    for (off, n_e, o_j) in ((o_csn, n_n, o_jn), (o_csp, n_p, o_jp)):
        for k in range(n_e):
            base = off + k*n_r
            for m in range(n_r):
                for mm in (m-1, m, m+1):
                    if 0 <= mm < n_r: add(base+m, base+mm)
            add(base + n_r-1, o_j + k)
    # 式2: c_e 三重対角 + j(電極ノード)
    for i in range(N):
        for ii in (i-1, i, i+1):
            if 0 <= ii < N: add(o_ce+i, o_ce+ii)
        if i < n_n: add(o_ce+i, o_jn+i)
        if i >= n_n+n_s: add(o_ce+i, o_jp+(i-n_n-n_s))
    # 式3: φ_s 三重対角(電極内)+ j
    for (off, n_e, o_j) in ((o_fsn, n_n, o_jn), (o_fsp, n_p, o_jp)):
        for k in range(n_e):
            for kk in (k-1, k, k+1):
                if 0 <= kk < n_e: add(off+k, off+kk)
            add(off+k, o_j+k)
    # 式4: φ_e 三重対角 + c_e 三重対角(κ(c_e)・ln項)+ j / 行0はゲージ
    add(o_fe+0, o_fsn+0)
    for i in range(1, N):
        for ii in (i-1, i, i+1):
            if 0 <= ii < N:
                add(o_fe+i, o_fe+ii); add(o_fe+i, o_ce+ii)
        if i < n_n: add(o_fe+i, o_jn+i)
        if i >= n_n+n_s: add(o_fe+i, o_jp+(i-n_n-n_s))
    # 式5: j は 自分・φ_s・φ_e・c_e・表面殻 c_s に依存
    for (o_j, n_e, o_fs, off_cs, xoff) in ((o_jn, n_n, o_fsn, o_csn, 0),
                                           (o_jp, n_p, o_fsp, o_csp, n_n+n_s)):
        for k in range(n_e):
            add(o_j+k, o_j+k); add(o_j+k, o_fs+k)
            add(o_j+k, o_fe + xoff + k); add(o_j+k, o_ce + xoff + k)
            add(o_j+k, off_cs + k*n_r + n_r-1)
    return np.array(rows), np.array(cols)

ROWS, COLS = build_pattern()

def color_columns():
    """同じ行を共有しない列は同色にできる(グラフ彩色の貪欲法)"""
    col_rows = [set() for _ in range(NV)]
    for r, c in zip(ROWS, COLS): col_rows[c].add(r)
    colors = -np.ones(NV, dtype=int)
    used_rows = []                    # 色ごとの占有行集合
    for c in range(NV):
        for col_id, occ in enumerate(used_rows):
            if not (col_rows[c] & occ):
                colors[c] = col_id; occ |= col_rows[c]; break
        else:
            colors[c] = len(used_rows); used_rows.append(set(col_rows[c]))
    return colors

COLORS_ = color_columns()
N_COLORS = COLORS_.max() + 1
# 変数のスケール(有限差分の刻み幅に使う)
SCALE = np.empty(NV)
SCALE[o_csn:o_ce] = 3e4; SCALE[o_ce:o_fsn] = 1e3
SCALE[o_fsn:o_jn] = 1.0; SCALE[o_jn:NV] = 1.0

def jacobian(y, y_old, dt, I_app, F0):
    """彩色付き有限差分で疎ヤコビアンを構築"""
    data_rows, data_cols, data_vals = [], [], []
    col_rows = {}
    for r, c in zip(ROWS, COLS): col_rows.setdefault(c, []).append(r)
    for g in range(N_COLORS):
        cols_g = np.where(COLORS_ == g)[0]
        h = 1e-7 * SCALE[cols_g]
        yp = y.copy(); yp[cols_g] += h
        F1 = residual(yp, y_old, dt, I_app)
        dF = F1 - F0
        for c, hc in zip(cols_g, h):
            for r in col_rows[c]:
                data_rows.append(r); data_cols.append(c)
                data_vals.append(dF[r]/hc)
    return sparse.csr_matrix((data_vals, (data_rows, data_cols)), shape=(NV, NV))

# ============ 5. Newton 法と時間積分 ============
def newton(y, y_old, dt, I_app, tol=1e-8, maxit=10):
    for it in range(maxit):
        F0 = residual(y, y_old, dt, I_app)
        if not np.all(np.isfinite(F0)): return y, False
        J = jacobian(y, y_old, dt, I_app, F0)
        dy = spsolve(J, -F0)
        y = y + dy
        if np.max(np.abs(dy)/SCALE) < tol: return y, True
    return y, False

def terminal_voltage(y, I_app):
    _, _, _, fsn, fsp, _, _, _ = unpack(y)
    # セル中心 → 集電体面へ半セル分だけオームの法則で外挿
    phi0 = fsn[0] + I_app*dx[0]/(2*sig_n)
    phiL = fsp[-1] - I_app*dx[-1]/(2*sig_p)
    return phiL - phi0

def initial_state(I_app):
    y = np.empty(NV)
    y[o_csn:o_csp] = cs0_n; y[o_csp:o_ce] = cs0_p; y[o_ce:o_fsn] = ce0
    un, up = U_n(cs0_n/cs_max_n), U_p(cs0_p/cs_max_p)
    y[o_fsn:o_fsp] = 0.0; y[o_fsp:o_fe] = up - un; y[o_fe:o_jn] = -un
    y[o_jn:o_jp] = I_app/(a_n*L_n); y[o_jp:NV] = -I_app/(a_p*L_p)
    return y

def solve_dfn(crate=1.0, n_snap=8, verbose=True):
    I_app = crate * I_1C
    y = initial_state(I_app)
    t, dt = 0.0, 5.0
    t_hist, V_hist, snaps = [], [], []
    t_end = 1.1*3600/crate
    next_snap = np.linspace(0, 3600/crate, n_snap)
    k_snap = 0
    while t < t_end:
        y_new, ok = newton(y.copy(), y, dt, I_app)
        if not ok:
            dt /= 2
            if dt < 1e-3: raise RuntimeError("Newton 法が収束しません")
            continue
        y, t = y_new, t + dt
        V = terminal_voltage(y, I_app)
        t_hist.append(t); V_hist.append(V)
        if k_snap < n_snap and t >= next_snap[k_snap]:
            csn, csp, ce, fsn, fsp, fe, jn, jp = unpack(y)
            snaps.append(dict(t=t, ce=ce.copy(), phie=fe.copy(),
                              jn=jn.copy(), jp=jp.copy(),
                              thn=cs_surf_of(csn, jn, Ds_n, dr_n)/cs_max_n,
                              thp=cs_surf_of(csp, jp, Ds_p, dr_p)/cs_max_p))
            k_snap += 1
        if V <= V_CUT:
            if verbose: print(f"  終止電圧到達: t = {t:.0f} s")
            break
        dt = min(dt*1.3, 30.0)
    cap = np.array(t_hist)*I_app*A_cell/3600      # 放電容量 [Ah]
    return np.array(t_hist), np.array(V_hist), cap, snaps

if __name__ == "__main__":
    import time
    print(f"未知数: {NV}(彩色数 {N_COLORS} → ヤコビアン1回 = 残差 {N_COLORS} 回分)")
    t0 = time.time()
    t, V, cap, snaps = solve_dfn(crate=1.0)
    print(f"1C 放電: {len(t)} ステップ, {time.time()-t0:.1f} 秒")
    print(f"放電容量: {cap[-1]:.3f} Ah, 最終電圧: {V[-1]:.3f} V")
    # 結果を保存(検証スクリプトで PyBaMM と比較する)
    np.savez("dfn_scratch_result.npz", t=t, V=V, cap=cap)
    # 可視化(ローカルで):
    # import matplotlib.pyplot as plt
    # plt.plot(cap, V); plt.xlabel("容量 [Ah]"); plt.ylabel("電圧 [V]"); plt.show()

5.4.4 実行方法と実行時間の目安

このコードはブラウザ(Pyodide)では実行しない。グリッドを教材用に粗くしても 1回の放電シミュレーションに数十秒かかり、ブラウザ内実行には重すぎるためである (この章の実行ブロックに「▶ ブラウザで実行」ボタンがないのはそのため)。 手元の Python 環境で走らせてほしい:

1. Python 3.10 以降を用意し、依存パッケージを入れる: pip install numpy scipy matplotlib
2. 上のリスティングを dfn.py として保存する。
3. python dfn.py で実行。既定では 1C 放電を解き、 $V(t)$ と各スナップショットを出力する。

実行時間の目安は、ふつうのノート PC で 1C 放電(約 3600 s の物理時間)が数十秒程度。 時間の大半はヤコビアン構築(有限差分の残差評価)が占める。 レートを上げると放電時間が短くなるぶんステップ数が減り、計算も速く終わる。 もし発散する・Newton が収束しない場合は、$\Delta t$ を半分にする、 摂動幅 $\varepsilon$ を変数のスケールに合わせる(濃度は $10^{-2}$、電位は $10^{-8}$ など相対摂動にする)、 初期化(無矛盾初期化)が成功しているか確認する、の順に疑うとよい。

5.5 PyBaMM で答え合わせ

5.5.1 PyBaMM とは

PyBaMM(Python Battery Mathematical Modelling)は、 電池の連続体モデルを解くためのオープンソース Python ライブラリである。 SPM・SPMe・DFN をはじめとする主要モデルと、文献由来のパラメータセット (この教材が使ってきた Chen et al. 2020 の LG M50 セルは "Chen2020" という名前で収録されている)を備え、 空間離散化(既定は有限体積法)と DAE ソルバ(CasADi / SUNDIALS)を自動で面倒みてくれる。 研究の現場では「自作コードの検証相手」あるいは「そのまま本番採用」される標準的な道具である。

私たちにとっての PyBaMM の役割は審判である。フェーズ4の 4.7 節で立てた検証の作法 — 「自作コードは、独立に実装された信頼できる実装と、同一条件で突き合わせて初めて信用できる」 — をここで実行する。インストールは1行:

pip install pybamm

5.5.2 最小コード — Chen2020 で DFN を解く

まず PyBaMM 単体で DFN を解いてみる。モデル選択・パラメータ・実験条件(1C 放電、2.5 V 終止)を 指定して解くだけなら、実質 20 行で書ける:

# PyBaMM で DFN を解く最小例(Chen2020 = LG M50、1C 放電、2.5 V 終止)
# 事前に: pip install pybamm
import pybamm

model = pybamm.lithium_ion.DFN()                # DFN モデル
params = pybamm.ParameterValues("Chen2020")     # Chen et al. 2020 パラメータ

experiment = pybamm.Experiment(["Discharge at 1C until 2.5 V"])
sim = pybamm.Simulation(model, parameter_values=params, experiment=experiment)
sol = sim.solve()

t   = sol["Time [s]"].entries
V   = sol["Voltage [V]"].entries
cap = sol["Discharge capacity [A.h]"].entries
print(f"放電時間   : {t[-1]:.0f} s")
print(f"終止電圧   : {V[-1]:.3f} V")
print(f"放電容量   : {cap[-1]:.3f} Ah")

# 対話的な可視化(ローカルで実行すると別ウィンドウが開く)
sim.plot(["Voltage [V]", "Electrolyte concentration [mol.m-3]"])

pybamm.Experiment に英語の文字列で運転条件を書くと、C レート換算・終止条件の監視まで やってくれる。イベント検出(電圧が 2.5 V を横切った瞬間で積分を止める)は フェーズ4の 4.6.4 節で自作したものと同じ機能である。

5.5.3 C レートを変えて放電曲線を重ねる

# C レート 0.5 / 1 / 2 / 3C の放電曲線を重ねる(PyBaMM, Chen2020)
import pybamm
import matplotlib.pyplot as plt

params = pybamm.ParameterValues("Chen2020")
fig, ax = plt.subplots(figsize=(6, 4))

for crate in [0.5, 1, 2, 3]:
    model = pybamm.lithium_ion.DFN()
    exp = pybamm.Experiment([f"Discharge at {crate}C until 2.5 V"])
    sim = pybamm.Simulation(model, parameter_values=params, experiment=exp)
    sol = sim.solve()
    cap = sol["Discharge capacity [A.h]"].entries
    V = sol["Voltage [V]"].entries
    ax.plot(cap, V, label=f"{crate}C")

ax.set_xlabel("Discharge capacity [Ah]")
ax.set_ylabel("Terminal voltage [V]")
ax.legend()
ax.grid(alpha=0.3)
plt.tight_layout()
plt.show()

レートを上げるほど電圧曲線全体が下がり(オーム損・過電圧の増大)、 かつ取り出せる容量そのものが減る(電解液の輸送律速)ことが確認できるはずである。 この「高レートで容量が急減する」挙動こそ SPM では再現できない DFN の見せ場であり、 5.6〜5.7 節でデータとともに掘り下げる。

5.5.4 教材コード vs PyBaMM — RMSE で定量比較

いよいよ答え合わせである。5.4 節の教材コードの出力(ここでは np.savez で $t$ と $V$ を保存してあるとする)と、同一条件の PyBaMM 解を共通の時間グリッドに補間して重ね、 二乗平均平方根誤差 RMSE を計算する:

# 教材 DFN(5.4 節)と PyBaMM の電圧曲線を重ねて RMSE を出す
import numpy as np
import pybamm
import matplotlib.pyplot as plt

# --- (1) 教材コードの結果を読み込む ---
# 5.4 節の教材コードは実行時に dfn_scratch_result.npz(キー: t, V, cap)を保存する
data = np.load("dfn_scratch_result.npz")
t_my, V_my = data["t"], data["V"]

# --- (2) PyBaMM を同一条件で解く ---
model = pybamm.lithium_ion.DFN()
params = pybamm.ParameterValues("Chen2020")
exp = pybamm.Experiment(["Discharge at 1C until 2.5 V"])
sim = pybamm.Simulation(model, parameter_values=params, experiment=exp)
sol = sim.solve()
t_pb = sol["Time [s]"].entries
V_pb = sol["Voltage [V]"].entries

# --- (3) 共通の時間グリッドへ線形補間して RMSE ---
t_end = min(t_my[-1], t_pb[-1])          # 早く終止した方に合わせる
t_common = np.linspace(0.0, t_end, 400)
V1 = np.interp(t_common, t_my, V_my)
V2 = np.interp(t_common, t_pb, V_pb)
rmse = np.sqrt(np.mean((V1 - V2) ** 2))
print(f"RMSE = {1000 * rmse:.1f} mV")

# --- (4) 重ね描き ---
plt.plot(t_my / 60, V_my, label="educational DFN (5.4)")
plt.plot(t_pb / 60, V_pb, "--", label="PyBaMM DFN (Chen2020)")
plt.xlabel("Time [min]")
plt.ylabel("Terminal voltage [V]")
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

教材コードのグリッド($n_n=12,\ n_s=6,\ n_p=12,\ n_r=12$、$\Delta t = 10$ s)では、 1C 放電の RMSE はおおむね10 mV 以下に収まる。 完全一致しないのは、(i) 空間・時間解像度の差(PyBaMM の既定はより細かい)、 (ii) 時間積分の次数(教材は1次の後退オイラー、PyBaMM は可変次数 BDF)、 (iii) 終止イベント検出の精度、が主因である。 グリッドと $\Delta t$ を細かくしたときに RMSE が系統的に減っていくなら、 差は「離散化誤差」であって「式の間違い」ではないと判断できる — これがフェーズ4の 4.7 節で学んだ 収束による切り分けである。

5.6 DFN の結果を読む — インタラクティブダッシュボード

ここからは解を「読む」練習である。下の図5.1 は、PyBaMM(Chen2020)で事前計算した 0.5C / 1C / 2C / 3C の定電流放電を、時間スライダーで自由に行き来できるダッシュボードである。 4枚のパネルはそれぞれ (a) 端子電圧、(b) 電解液濃度 $c_e(x)$、(c) 表面充填率 $\theta_{\mathrm{surf}}(x) = c_{s,\mathrm{surf}}/c_{s,\max}$、(d) 反応電流密度 $j(x)$ を示す。 横軸はセル内の規格化位置 $x/L$ で、背景の色分けは負極(青)/セパレータ(灰)/正極(赤)である。

このダッシュボードには、SPM/SPMe では原理的に存在しなかった情報が詰まっている。 順に読み解こう。

5.6.1 (i) 高レートで正極の $c_e$ が枯渇する

パネル (b) で 3C を選ぶと、放電が進むにつれて $c_e(x)$ が大きく傾き、 正極の集電体側($x/L = 1$ 付近)で初期値 1000 mol/m³ からほぼゼロまで払底するのが見える (放電終盤のスナップショットでは、正極内に $c_e \approx 0$ の区間が広がる)。 放電中、Li⁺ は負極で生成され(酸化反応 $j>0$)正極で消費される($j<0$)。 消費に補給が追いつくかどうかは電解液内の拡散・泳動の速さ次第で、 レートが上がると輸送が追いつかず、正極の奥から塩が払底していく。 $c_e$ が下がると、その場のイオン伝導率 $\kappa^{\mathrm{eff}}(c_e)$ も交換電流密度 $i_0 \propto c_e^{1/2}$ も同時に下がるため、電圧損失が雪だるま式に増える。 3C の電圧曲線(パネル a)が終盤で急落し容量が伸びない直接の原因がこれである。 フェーズ3の SPMe も $c_e(x)$ を持っていたが、ソース項が一様と仮定していたため 枯渇が反応分布に跳ね返るフィードバックは表現できなかった。

5.6.2 (ii) 反応電流 $j(x)$ がセパレータ側に偏る — 反応の局在化

パネル (d) を見ると、低レート(0.5C)では $j(x)$ は各電極内でほぼ平坦で、 SPMe の「一様反応」仮定が良い近似であることがわかる。 ところが 3C では、放電初期から両電極ともセパレータに近い側で $|j|$ が大きい。 理由は電流の通り道を考えるとわかる。反応がセパレータから遠い側(集電体側)で起こると、 Li⁺ はセパレータまでの長い距離を液相で運ばれなければならず、液相のオーム損 ($\kappa^{\mathrm{eff}}$ は固相の $\sigma^{\mathrm{eff}}$ よりずっと小さい)が大きくつく。 系は全体の損失が小さくなるように電流経路を配分するので、 液相を通る距離が短くて済むセパレータ側に反応が集中するのである。 さらに放電後半の高レートでは、正極奥の $c_e$ 枯渇(見どころ i)が奥での反応を物理的に不可能にし、 局在化が一層強まる。時間を動かしながら (b) と (d) を見比べてほしい。

5.6.3 (iii) 表面充填率 $\theta_{\mathrm{surf}}(x)$ の不均一

パネル (c) は各 $x$ の粒子の「表面の」充填率である。反応が不均一(見どころ ii)なら、 Li の出入りも不均一なので、同じ電極の中でも粒子ごとに充電状態がずれていく。 3C の放電後半では、負極内で $\theta_{\mathrm{surf}}$ がセパレータ側と集電体側で大きく違う値を取る。 これは同じ電極の中に SOC の違う粒子が共存しているということであり、 「電極=1個の代表粒子」という SPM/SPMe の描像がもはや成立していないことの直接の証拠である。 なお放電終了後に電流を切って放置すると、この不均一は粒子間の緩和(液相を介した再配分)で ゆっくり解消していく — 実験で放電直後の開回路電圧がしばらく浮き上がる現象の一因である。

5.6.4 (iv) 電圧の内訳を読む

パネル (a) の電圧曲線を、他の3パネルと突き合わせると「電圧がどこで失われているか」を読み取れる。 端子電圧はおおまかに

と分解できる(放電ではどの損失項も電圧を押し下げる向きに働く)。 低レートの電圧曲線はほぼ OCV そのもの(熱力学項が支配)。 レートを上げたときの下げ幅のうち、放電の初期から一定量ある分は反応過電圧と固相・液相オーム損、 放電が進むにつれて増える分は $c_e$ 勾配の発達による濃度過電圧 ($\Delta\phi_e$ の中の $\ln c_e$ 項、式4の第2項)である。 3C で終盤に電圧が崖のように落ちるのは、OCV の坂($\theta$ が端に近づく)と $c_e$ 枯渇による損失増大が同時に来るためで、(b)(c) を見ればどちらがどれだけ効いているか判別できる。 この「内訳を読む」技術は、実測の放電曲線から律速過程を推定するときの基本動作になる。

5.7 モデル階層の答え合わせ — SPM / SPMe / DFN

最後に、フェーズ2〜5で積み上げてきたモデル階層を同じ土俵で比較する。 図5.2 は、同一パラメータ(Chen2020)・同一条件で SPM・SPMe・DFN を解いた放電曲線の重ね描きである。 下段には DFN を基準とした差 $\Delta V = V_{\mathrm{model}} - V_{\mathrm{DFN}}$ を示す。

読み取れることを整理する。1C では SPMe と DFN はほとんど区別がつかず(最大でも 10 mV 程度)、 SPM のずれも数十 mV に収まる。反応分布はほぼ一様(5.6.2 節)なので、 一様反応を仮定した簡約モデルが十分に機能する領域である。 3C では事情が一変し、しかも2つの簡約モデルは逆方向に破綻する。 電解液を持たない SPM は、$c_e$ 枯渇による損失と容量制限をまったく表現できないため、 電圧を最大 1 V 近くも高く見積もり、DFN が 2.3 Ah で終止するところを 4.5 Ah まで放電できると答えてしまう(過大評価)。 一方 SPMe は、$c_e$ の PDE は持っているのに反応を一様のまま固定しているため、 正極の奥まで一定ペースで塩を消費し続ける。奥の $c_e$ が払底すると濃度過電圧の補正項 ($\ln c_e$ を含む)が発散気味に効いて電圧が崩落し、わずか 0.25 Ah で 2.5 V に達して終止する(極端な過小評価)。 実際のセル(と DFN)がそこで止まらないのは、5.6.2 節で見たとおり 反応が電解液の残っているセパレータ側へ移動して放電を続けられるからで、 この「反応の再配置」という自由度こそ SPMe に欠けているものである。 つまりモデル階層は「どこまでの物理を諦めるか」の階層であり、 諦めた物理が効かない運転条件でだけ、簡約モデルは正しい。 しかも適用域の外での壊れ方は、SPM のように楽観側とは限らず、SPMe のように突然かつ悲観側のこともある。 「どの向きに、いつ壊れるか」は DFN と突き合わせて初めてわかる — 図5.2 はその実演である。

5.8 理解度チェックと次章への橋渡し