フェーズ2単粒子モデル(SPM)の導出と実装

フェーズ1で組み上げた DFN モデルは、5 本の偏微分方程式が絡み合う大きな系でした。 この章では、DFN に 3 つの近似を入れて得られる最小の電池モデル 「単粒子モデル(Single Particle Model, SPM)」を導出します。 SPM は「各電極を代表粒子 1 個で表す」という大胆な簡約ですが、 低レートでは驚くほど良く実測を再現し、しかも球の拡散方程式 2 本と代数式だけで書けるため、 手を動かして電池モデルの数値計算を学ぶのに最適な題材です。 章の後半では、ブラウザ上で動く SPM シミュレータと、 そのまま Jupyter でも動く Python 実装を作ります。

2.1 SPM の仮定 — 何を無視するのか

DFN モデル(フェーズ1)は、セル厚み方向の位置 $x$ [m] ごとに 粒子を置き、電解液の濃度 $c_e(x,t)$ [mol/m³] と電位 $\phi_e(x,t)$ [V]、 固相電位 $\phi_s(x,t)$ [V]、界面電流密度 $j(x,t)$ [A/m²] をすべて解きました。 SPM はこのうち「電解液側の物理」と「電極内の位置依存性」をまとめて捨てます。 まず、何を仮定するのかを明示的に並べます。

それぞれの仮定が「いつ妥当か」を、LG M50 セル(Chen et al. 2020、 js/params.js の値)で定量的に見積もっておきます。 ここで使う実効物性は Bruggeman 補正 $\kappa^{\mathrm{eff}} = \kappa\,\varepsilon_e^{1.5}$, $D_e^{\mathrm{eff}} = D_e\,\varepsilon_e^{1.5}$ です($\kappa$ [S/m]:液相イオン伝導率、$D_e$ [m²/s]:電解液拡散係数、 $\varepsilon_e$ [−]:空隙率)。

仮定 (i) の目安 — 低 C レート・薄い電極

電解液に濃度勾配ができる速さと大きさを見積もります。 セル全厚 $L = L_n + L_s + L_p = 172.8$ μm、 $D_e(c_{e0}) = 1.77\times10^{-10}$ m²/s、セパレータの $D_e^{\mathrm{eff}} \approx 5.7\times10^{-11}$ m²/s を使うと、 電解液が定常濃度分布に達する時定数は

つまり放電開始から 10 分ほどで電解液には定常的な濃度勾配が立ちます。 その大きさは、供給すべき Li⁺ 流束と拡散のつり合いから (輸率 $t_+^0 = 0.259$ [−]、$F = 96485$ C/mol、印加電流密度 $I$ [A/m²]、放電を正)

と見積もれます($I_{1C} = 48.7$ A/m² は 1C 相当の電流密度)。 1C では $\Delta c_e$ が初期濃度 1000 mol/m³ の半分を超えるオーダーであり、 仮定 (i) が安心して使えるのは目安として 0.5C 以下、厳密には 0.2C 程度までです。 また $\Delta c_e \propto L$ なので、電極が薄いセル(高出力設計)ほど SPM は良く当たります。 液相電位の方は、1C・$\kappa^{\mathrm{eff}} \approx 0.12$–$0.31$ S/m で 電位降下が電極あたり 10–18 mV 程度(後述)。 端子電圧 3–4 V に対しては小さいものの、無視できるのはやはり低レートに限られます。

仮定 (ii) の目安 — 電位降下が反応の「感度」より小さいこと

$j(x)$ が一様になるのは、電極内で $\phi_s - \phi_e$ がほぼ一様なときです (界面反応は $\eta = \phi_s - \phi_e - U$ で駆動されるため)。 一様 $j$ のときの液相電位降下は電極厚さ $L^\pm$ にわたって

固相側は負極では $\sigma^{\mathrm{eff}} = 140$ S/m と大きく降下は 0.01 mV 程度で完全に無視できますが、 正極は $\sigma^{\mathrm{eff}} = 0.098$ S/m しかなく約 19 mV(1C)になります。 これらを反応の感度スケール $RT/F = 25.7$ mV(室温)と比べると、 1C ですでに同程度です。 つまり厳密には 1C で $j$ の非一様性が現れ始め、 それ以下のレート(目安 0.5C 以下)なら一様とみなせます。 伝導度が高い(または電極が薄い)ほどこの条件は緩みます。

仮定 (iii) は (i)(ii) の帰結

(i)(ii) が成り立てば、電極内の全粒子は「同じ表面流束・同じ初期濃度」なので 区別がつきません。粒径分布も無視して代表半径 $R_s^\pm$ [m] の球 1 個に集約します。 以上の適用範囲をまとめると:

2.2 DFN からの簡約 — 5 式のうち何が消えるか

フェーズ1で導出した DFN の支配方程式 5 本を再掲します(導出はフェーズ1参照)。 記号はすべてフェーズ1と共通です: $c_s(x,r,t)$ [mol/m³] は固相 Li 濃度、$r$ [m] は粒子中心からの半径、 $D_s$ [m²/s] は固相拡散係数、$a = 3\varepsilon_s/R_s$ [1/m] は比界面積 ($\varepsilon_s$ [−]:活物質体積分率)、 $i_s, i_e$ [A/m²] は固相・液相の見かけ電流密度、 $\sigma^{\mathrm{eff}}$ [S/m] は実効電子伝導率、 $i_0$ [A/m²] は交換電流密度、$k$ [m²·⁵/(mol⁰·⁵·s)] は反応速度定数、 $\alpha_a, \alpha_c$ [−] は移動係数(本教材では 0.5, 0.5)、 $R = 8.314$ J/(mol·K)、$T = 298.15$ K です。 座標は負極集電体を $x = 0$、放電電流を正 $I > 0$、 酸化(粒子から Li が出る)方向の界面電流を正 $j > 0$ とします。

① 固相拡散(各 $x$ の粒子内):

② 液相物質保存:

③ 固相電荷保存:

④ 液相電荷保存(修正オームの法則。熱力学因子は 1 とする):

⑤ 界面反応(Butler–Volmer):

どの $x$ でも $i_s + i_e = I$(電荷保存則)が成り立つことも思い出しておきます。 ここに §2.1 の仮定を順番に入れていきます。

ステップ 1: 仮定 (i) で液相の 2 式が消える

$c_e(x,t) = c_{e0}$(定数)と置くと、(DFN-2) の左辺と拡散項はゼロになります。 厳密にはこのとき生成項 $(1-t_+^0)aj/F$ が残って矛盾しますが、 仮定 (i) の主張は「この生成項が作る濃度変化 $\Delta c_e$ が $c_{e0}$ に比べて無視できるほど小さい」(§2.1 の見積もり)ということです。 つまり (DFN-2) は「解かない」と決めるのが SPM です。

(DFN-4) では $\partial \ln c_e/\partial x = 0$ となるので濃度項が消え、 $i_e = -\kappa^{\mathrm{eff}}\,\partial\phi_e/\partial x$ という普通のオームの法則に戻ります。 さらに $\partial\phi_e/\partial x = 0$ とみなす(伝導率が実質無限大という理想化)ので、 $\phi_e(x) =$ 定数となり (DFN-4) も消えます。 ただし液相電流 $i_e$ 自体は流れていることに注意してください — 「勾配ゼロで電流が流れる」のは $\kappa^{\mathrm{eff}} \to \infty$ の極限の理想化です。 このあと §2.3 で使うのは (DFN-4) の後半、 電荷保存の微分形 $\partial i_e/\partial x = a j$ だけです。 これは電位とは無関係に成り立つ「電流の帳簿」なので、SPM でも生き残ります。

ステップ 2: 固相電荷保存が自明になる

(DFN-3) は $\phi_s(x)$ の分布を決める式でした。 負極では $\sigma^{\mathrm{eff}} = 140$ S/m と大きく、 1C での電位降下は $\sim 0.01$ mV(§2.1)。 よって $\phi_s$ は各電極内で一様と置けます。 $\phi_s$ が一様なら (DFN-3) の左辺の意味は失われ、 残る内容は「両辺を電極厚さで積分した電流の帳簿」だけになります。 実際、$i_s + i_e = I$ と $\partial i_e/\partial x = aj$ から $\partial i_s/\partial x = -aj$ であり、これを負極全体 $[0, L_n]$ で積分すると

(負極集電体 $x=0$ では電流はすべて固相なので $i_s(0) = I$、 セパレータ境界 $x=L_n$ では固相電流は流れられないので $i_s(L_n) = 0$)。 つまり「電極内の反応電流を全部足すと印加電流になる」という当然の帳簿だけが残ります。 これが「固相電荷保存が自明になる」の意味です。

ステップ 3: 対応表

2.3 一様反応流束の導出 — $j^\pm = \pm I/(a^\pm L^\pm)$

SPM の入力は印加電流密度 $I$ [A/m²](放電を正)だけです。 代表粒子の境界条件になる $j^\pm$ を、電荷保存と境界条件から 1 ステップずつ求めます。 使うのは (DFN-4) の後半 $\partial i_e/\partial x = aj$ と、 「集電体では電流はすべて固相、セパレータでは電流はすべて液相」という境界条件だけです。

単位確認:$[I] = \mathrm{A/m^2}$、$[a] = 1/\mathrm{m}$、$[L] = \mathrm{m}$ なので $a L$ は無次元、したがって $[j] = \mathrm{A/m^2}$ ✓。 $j$ は「活物質表面 1 m² あたりの電流」、$I$ は「セル断面 1 m² あたりの電流」で、 同じ単位ですが面積の基準が違うことに注意してください。 その換算係数が無次元量 $aL$ です。

2.4 SPM の方程式系 — 球拡散と平均濃度の帳簿

材料が出そろいました。SPM の全方程式系を明示します。 上付き $\pm$ は正極・負極を表し、$-$ が負極です。 未知関数は $c_s^-(r,t)$ と $c_s^+(r,t)$ の 2 つだけ、 残りはすべて代数式(その場で計算できる式)です。

平均濃度の時間発展 — 積分するだけで厳密に出る保存則

(2.3) を数値で解く前に、手で厳密に分かることを 1 つ引き出しておきます。 粒子の体積平均濃度

がどう時間変化するかは、拡散方程式を体積積分するだけで求まります。 途中の積分を省略せずに追います。

2.5 Butler–Volmer の逆変換 — $\eta$ を $j$ から求める

DFN では $j$ と $\eta$ が互いに絡む非線形方程式でしたが、 SPM では $j^\pm$ が式 (2.1)(2.2) で先に決まってしまうので、 Butler–Volmer 式は「$\eta$ を求めるための式」に役割が変わります。 幸い $\alpha_a = \alpha_c = 0.5$ のとき、この逆変換は閉じた形で書けます。

交換電流密度 $i_0$ [A/m²] は表面濃度に依存します (params.js の i0_fun と同じ式。$c_e = c_{e0}$ 固定):

$k^-$ = 6.716×10⁻¹² m²·⁵/(mol⁰·⁵·s)、$k^+$ = 3.545×10⁻¹¹ m²·⁵/(mol⁰·⁵·s) です。 形に注目してください:$i_0$ は $c_{s,\mathrm{surf}} \to 0$(表面に Li がない=酸化できない)でも $c_{s,\mathrm{surf}} \to c_{s,\max}$(表面が満杯=還元できない)でもゼロに落ちます。 $i_0 \to 0$ のとき (2.8) の $\mathrm{asinh}$ の引数が発散し $|\eta| \to \infty$、つまり表面の枯渇・飽和は過電圧の暴走として電圧に現れます。 これが図2.2 で見る「高レートで電圧が急落する」メカニズムの正体です。

満充電(100% SOC)の LG M50 に 1C を流した瞬間の値を入れてみます:

• 負極:$\theta_n = 0.9014$ → $i_0^- = 0.202$ A/m²。 $\dfrac{j^-}{2i_0^-} = \dfrac{1.49}{0.404} = 3.68$ → $\eta^- = 51.4\,\mathrm{mV} \times \mathrm{asinh}(3.68) = +103$ mV
• 正極:$\theta_p = 0.27$ → $i_0^+ = 3.03$ A/m²。 $\dfrac{j^+}{2i_0^+} = -0.278$ → $\eta^+ = -14$ mV

同じ 1C でも、反応が遅い負極(黒鉛)の過電圧は正極の 7 倍あります。 このセルの反応抵抗のボトルネックは負極、 拡散のボトルネックは正極(§2.4)という非対称な分担になっています。

2.6 端子電圧の組み立て — OCV・動作窓・容量の換算

最後のパーツは端子電圧 $V$ [V] です。定義は $V = \phi_s(L) - \phi_s(0)$(正極集電体の固相電位 − 負極集電体の固相電位)でした。 過電圧の定義 $\eta = \phi_s - \phi_e - U(\theta)$ を各電極で $\phi_s = \phi_e + U + \eta$ と書き直します。 §2.2 より $\phi_s$ は各電極内で一様、$\phi_e$ は全域で一様(定数)なので:

定数 $\phi_e$ が差し引きで消えて($\phi_e$ を測る基準がどこでも良い理由がこれです)、 SPM の電圧式が得られます:

$U_n, U_p$ [V] は各電極の開回路電位(OCP)、 $\theta^{\mathrm{surf}} = c_{s,\mathrm{surf}}/c_{s,\max}$ [−] は表面の化学量論です。 平均でなく表面の値を使う点が重要です(§2.4)。

OCV 曲線と動作窓

電流を流さなければ($I = 0$ → $j^\pm = 0$ → $\eta^\pm = 0$、 十分待てば表面=平均)、(2.10) は開回路電圧(OCV)になります: $V_{\mathrm{OCV}} = U_p(\theta_p) - U_n(\theta_n)$。 ここで $U_n(\theta)$、$U_p(\theta)$ は params.js の ocpN/ocpP(Chen et al. 2020 の実測フィット)です。 黒鉛の $U_n$ は 0.1–0.2 V 付近に平坦部(ステージ構造と呼ばれる Li の秩序配列に対応)を持ち、 $\theta \to 0$ で急上昇します。NMC811 の $U_p$ は 3.5–4.3 V をなだらかに変化します。

重要なのは、各電極は $\theta$ の全域 $[0, 1]$ を使わないことです。 セルの充放電は端子電圧の上下限(このセルでは 4.2 V / 2.5 V)で区切られ、 それに対応する $\theta$ の範囲(動作窓)は:

SOC(State of Charge、充電率)[−] と $\theta$ の対応は、 この端点間の線形補間で定義します:

SOC が上がると負極には Li が入り($\theta_n$ 増)、正極からは抜けます($\theta_p$ 減) — 2 つの $\theta$ は常に逆方向に動きます。図2.1 で実際に動かして確かめてください。

セル容量と θ の換算 — 化学量論の帳簿付け

「θ が 0.9 から 0.03 まで動く」ことと「セルから 5 Ah 取り出せる」ことを結びつけます。 電極(面積 $A$ [m²]、厚さ $L^\pm$ [m]、活物質体積分率 $\varepsilon_s^\pm$ [−])に 入っている Li の物質量は $n_{\mathrm{Li}}^\pm = \varepsilon_s^\pm L^\pm A\, c_{s,\max}^\pm\, \theta^\pm$ [mol]。 $\theta$ が動作窓を端から端まで動くときに出入りする電荷が、その電極の容量です:

LG M50 の値($A = 0.1027$ m²)を入れると:

• 負極:$Q^- = \dfrac{96485 \times 0.75 \times 85.2\,\mathrm{μm} \times 0.1027 \times 33133 \times 0.8735}{3600} = 5.09$ Ah
• 正極:$Q^+ = \dfrac{96485 \times 0.665 \times 75.6\,\mathrm{μm} \times 0.1027 \times 63104 \times 0.6384}{3600} = 5.58$ Ah

公称容量 5.0 Ah にいちばん近いのは負極の窓で、 このセルの放電は負極側の窓で頭打ちになる(負極律速)設計です。 また、1C 電流密度が $I_{1C} = \dfrac{5.0\ \mathrm{Ah}}{0.1027\ \mathrm{m^2}} = 48.7$ A/m² という params.js の定義もここから来ています。 平均化学量論の時間変化は、(2.5) を $c_{s,\max}$ で割れば $d\langle\theta\rangle/dt = -3j/(F R_s c_{s,\max})$ となり、 定電流なら $\langle\theta\rangle$ は時間の 1 次関数です — 容量 [Ah] と $\langle\theta\rangle$ は完全に比例します。

2.7 数値解法(最小限)— シェルに分割する

SPM で数値計算が必要なのは球拡散 (2.3) だけです。 $r$ 方向に連続な関数 $c_s(r,t)$ をそのまま計算機に載せることはできないので、 $r$ 方向を離散化します。 ここでは実装の直感だけをつかみ、理論(離散化の精度・安定性・剛性)はフェーズ4に譲ります。

有限体積法(FVM)の直感 — タマネギの皮の保存則

粒子を $N_r$ 枚の同心シェル(タマネギの皮)に分割します。 シェル $i$($i = 0, \dots, N_r - 1$)は半径 $r_i = i\,\Delta r$ から $r_{i+1} = (i+1)\Delta r$ の球殻で、$\Delta r = R_s/N_r$。 各シェルについて「Li の増減 = 内側の面から入る流れ − 外側の面から出る流れ」 という保存則をそのまま書きます:

$c_i$ [mol/m³] はシェル $i$ の平均濃度、面を通るフラックスは フィックの法則 $-D_s\,\partial c/\partial r$ を隣接シェルの濃度差で近似したものです。 境界条件も同じ言葉で書けます: 最内シェルは内側の面積が $A_0 = 0$ なので流入項が自然に消え(対称条件)、 最外シェルの外側の面には拡散フラックスの代わりに 反応流束 $-\dfrac{j}{F} A_{N_r}$ [mol/s] を与えます(境界条件 (2.4))。

method of lines — 空間だけ離散化して ODE に落とす

(2.13) は $N_r$ 本の常微分方程式(ODE)の系 $\dfrac{d\mathbf{c}}{dt} = \mathbf{A}\mathbf{c} + \mathbf{b}$ です ($\mathbf{A}$:三重対角行列、$\mathbf{b}$:表面流束の項)。 このように「空間を先に離散化して、時間方向は ODE ソルバに任せる」戦略を method of lines(線の方法)と呼びます。 時間積分のやり方は 2 通り用意しました:

• Python 実装(§2.8):scipy.integrate.solve_ivp (陰的 BDF 法)に丸投げします。拡散方程式の ODE 系は「硬い(stiff)」 — 速い成分と遅い成分が混在する — ため、陽解法だと時間刻みを極端に 小さくしないと発散します。陰的ソルバを選ぶのが定石です(理由はフェーズ4で)。
• ブラウザ実装(図2.2):Crank–Nicolson 法(半陰的、2 次精度) + Thomas 法(三重対角ソルバ)で自前で解いています。 1 ステップが三重対角連立 1 回で済むため、スライダーを動かすたびに 全放電を再計算しても 50 ms 以内に収まります。

2.8 実装 — インタラクティブ SPM と Python コード

方程式系(§2.4)を実際に動かします。まずはシミュレータで SPM の挙動を体感してください。 設定は定電流放電:初期状態 100% SOC($\theta_n = 0.9014$, $\theta_p = 0.27$)、 下限電圧 2.5 V で打ち切りです。

いくつか試してほしい実験と、結果の読み方:

• 0.1C:電圧曲線は OCV 点線とほぼ重なり、容量は窓いっぱいの約 5.1 Ah。 過電圧も表面・平均の差もほぼゼロ。「電池の熱力学だけ」が見えている状態です。
• 1C(既定):曲線全体が OCV から 0.1–0.15 V 下に平行移動 (主に負極の反応過電圧 $\eta^- \approx +0.1$ V)し、容量は約 4.96 Ah。 正極プロファイルには目に見える勾配が立ちます。
• 5C:初期電圧から大きく下がり、約 3.5 Ah で急落。 終了時の読み出しで $\theta_p$ の表面が 1 近く、平均が 0.7 程度 — 「中身は空いているのに表面が満杯」の典型例です。
• $D_s$ 0.1 倍(1C):5C と同様の早期打ち切りが起きます。 レートを上げることと拡散を遅くすることは、粒子から見ると同じことです。

Python 実装(実行可能)

同じ SPM を、§2.7 の方針(FVM で $N_r = 20$ シェル + solve_ivp)で ゼロから実装したコードです。パラメータは js/params.js と同一の値を ハードコードしてあり、「▶ ブラウザで実行」でこのページ内で走ります (初回は Pyodide のダウンロードで数十秒かかります)。 そのままコピーすればローカルの Jupyter でも動きます (末尾のコメントに matplotlib 版の描画コードがあります)。

# =====================================================================
# 単粒子モデル(SPM)の完全実装 — LG M50 セル(Chen et al. 2020)1C 放電
# 必要ライブラリ: numpy, scipy(ローカル Jupyter でもそのまま動きます)
# =====================================================================
import numpy as np
from scipy.integrate import solve_ivp

# ---- 物理定数・セルパラメータ(js/params.js と同じ値) ----
F, R, T = 96485.33212, 8.314462618, 298.15   # C/mol, J/(mol K), K
A_cell  = 0.1027                              # 電極面積 [m^2]
I_1C    = 5.0 / A_cell                        # 1C 電流密度 [A/m^2] ≈ 48.69
ce0     = 1000.0                              # 電解液塩濃度 [mol/m^3](SPM では定数)

# 電極パラメータ: [負極 n, 正極 p]
L      = np.array([85.2e-6, 75.6e-6])         # 厚さ [m]
Rs     = np.array([5.86e-6, 5.22e-6])         # 粒子半径 [m]
eps_s  = np.array([0.75, 0.665])              # 活物質体積分率 [-]
Ds     = np.array([3.3e-14, 4.0e-15])         # 固相拡散係数 [m^2/s]
cs_max = np.array([33133.0, 63104.0])         # 最大 Li 濃度 [mol/m^3]
k_rate = np.array([6.716e-12, 3.545e-11])     # 反応速度定数 [m^2.5/(mol^0.5 s)]
th_100 = np.array([0.9014, 0.2700])           # θ @100% SOC(放電開始点)

a_s = 3 * eps_s / Rs                          # 比界面積 a = 3ε_s/R_s [1/m]

def ocpN(t):   # 黒鉛負極の OCP [V](Chen 2020 フィット)
    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 ocpP(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 i0_fun(k, ce, cs, csm):   # 交換電流密度 [A/m^2]
    return k * F * np.sqrt(ce) * np.sqrt(cs) * np.sqrt(csm - cs)

# ---- 放電条件と一様反応電流密度(2.3 節の式) ----
Crate = 1.0
I = Crate * I_1C                              # 放電を正 [A/m^2]
j = np.array([+I / (a_s[0] * L[0]),           # 負極: j > 0(Li が出る)
              -I / (a_s[1] * L[1])])          # 正極: j < 0(Li が入る)

# ---- r 方向の有限体積(FVM)離散化: 各電極 Nr シェル ----
Nr = 20
dr = Rs / Nr                                  # シェル幅 [m]
iarr = np.arange(Nr)
# シェル体積 V_i ∝ (r_{i+1}^3 - r_i^3)/3、面係数 G_i ∝ D_s r_{i+1}^2 / dr
# (共通因子 4π は分子分母で打ち消えるので省略)
def fvm_matrix(e):   # 電極 e(0=n, 1=p)の拡散行列 A(dc/dt = A c + b)
    Vsh = ((iarr + 1) ** 3 - iarr ** 3) * dr[e] ** 3 / 3.0
    G = Ds[e] * ((iarr[:-1] + 1) * dr[e]) ** 2 / dr[e]
    A = np.zeros((Nr, Nr))
    for i in range(Nr):
        if i > 0:
            A[i, i - 1] += G[i - 1] / Vsh[i]
            A[i, i] -= G[i - 1] / Vsh[i]
        if i + 1 < Nr:
            A[i, i + 1] += G[i] / Vsh[i]
            A[i, i] -= G[i] / Vsh[i]
    b = np.zeros(Nr)
    b[-1] = -(j[e] / F) * Rs[e] ** 2 / Vsh[-1]   # 表面フラックス項(境界条件)
    return A, b

A_n, b_n = fvm_matrix(0)
A_p, b_p = fvm_matrix(1)
J_full = np.zeros((2 * Nr, 2 * Nr))              # 定数ヤコビアン(BDF 用)
J_full[:Nr, :Nr], J_full[Nr:, Nr:] = A_n, A_p
b_full = np.concatenate([b_n, b_p])

def rhs(t, y):                                    # dy/dt = A y + b(線形!)
    return J_full @ y + b_full

def voltage(y):                                   # 状態 → 端子電圧 [V]
    cs_surf = np.array([y[Nr - 1], y[2 * Nr - 1]]) - j * dr / (2 * F * Ds)
    cs_surf = np.clip(cs_surf, 1e-3, cs_max - 1e-3)   # 数値ガード
    th = cs_surf / cs_max
    i0 = i0_fun(k_rate, ce0, cs_surf, cs_max)
    eta = (2 * R * T / F) * np.arcsinh(j / (2 * i0))  # BV の逆変換(2.5 節)
    return (ocpP(th[1]) + eta[1]) - (ocpN(th[0]) + eta[0])

def cutoff(t, y):                                 # V = 2.5 V で積分打ち切り
    return voltage(y) - 2.5
cutoff.terminal = True
cutoff.direction = -1

y0 = np.concatenate([np.full(Nr, th_100[0] * cs_max[0]),
                     np.full(Nr, th_100[1] * cs_max[1])])
t_max = 1.2 * 3600 / Crate
sol = solve_ivp(rhs, (0, t_max), y0, method="BDF", jac=lambda t, y: J_full,
                events=cutoff, t_eval=np.linspace(0, t_max, 400), rtol=1e-6)

# 打ち切り時刻の点を末尾に追加
t_arr = np.append(sol.t, sol.t_events[0])
y_arr = np.hstack([sol.y, sol.y_events[0].T])
V_arr = np.array([voltage(y_arr[:, i]) for i in range(y_arr.shape[1])])
cap = I * A_cell * t_arr / 3600.0                 # 放電容量 [Ah]

print(f"1C 放電: 打ち切り時刻 t = {t_arr[-1]:.1f} s ({t_arr[-1]/60:.1f} 分)")
print(f"放電容量 = {cap[-1]:.3f} Ah(公称 5.0 Ah)")
print(f"電圧範囲: {V_arr.min():.3f} V 〜 {V_arr.max():.3f} V")
th_sN = (y_arr[Nr-1, -1] - j[0]*dr[0]/(2*F*Ds[0])) / cs_max[0]
th_sP = (y_arr[2*Nr-1, -1] - j[1]*dr[1]/(2*F*Ds[1])) / cs_max[1]
print(f"終了時の表面化学量論: θ_n,surf = {th_sN:.3f}, θ_p,surf = {th_sP:.3f}")

# Plotly 描画(このページ用)。Jupyter では下の matplotlib 版を使ってください
plot_spec = {
    "traces": [{"x": cap.tolist(), "y": V_arr.tolist(),
                "name": "SPM 1C", "mode": "lines"}],
    "layout": {"title": {"text": "SPM 1C 放電曲線(LG M50)"},
               "xaxis": {"title": {"text": "放電容量 [Ah]"}},
               "yaxis": {"title": {"text": "端子電圧 V [V]"}}},
}
# --- matplotlib 版(Jupyter 用) ---
# import matplotlib.pyplot as plt
# plt.plot(cap, V_arr)
# plt.xlabel("放電容量 [Ah]"); plt.ylabel("端子電圧 V [V]")
# plt.title("SPM 1C 放電曲線"); plt.grid(True); plt.show()

コードの読みどころと期待される出力:

• 構造が方程式系と 1 対 1: $j^\pm$ の計算(§2.3)→ fvm_matrix が式 (2.13)(§2.7)→ voltage が BV 逆変換 (2.8) と電圧式 (2.10)(§2.5–2.6)。 教科書の式とコードの行が対応していることを確認してください。
• 系が線形:rhs は J_full @ y + b_full だけ。 定数ヤコビアンを solve_ivp に渡しているので BDF 法が高速に走ります。
• イベント検出:cutoff 関数の符号が変わった時刻を ソルバが自動で見つけ、そこで積分を打ち切ります(terminal = True)。
• 出力の目安:放電容量 ≈ 4.955 Ah、電圧範囲 2.500–4.038 V、 終了時 $\theta_n^{\mathrm{surf}} \approx 0.035$(負極表面がほぼ枯れて終了)。 図2.2(Crank–Nicolson、$N_r = 30$)の 1C の結果 4.96 Ah と 0.2% 以内で一致します — 別々に実装した 2 つの解法が一致するのは良い検算です。

2.9 SPM の限界 — どこから使えなくなるか

SPM の誤差は、§2.1 で捨てた「電解液の物理」がそのまま返ってきます。 完全 DFN(フェーズ5で計算します)と比べると、LG M50 では概ね次のようになります:

なぜずれるのかを §2.1 の見積もりで説明できます。 SPM が無視した電圧損失は主に 2 つ、いずれも電解液由来です(総称して電解液分極):

1. 液相オーム損: $\Delta\phi_e \sim \dfrac{I L}{2\kappa^{\mathrm{eff}}}$。 1C で電極・セパレータ合計 約 30 mV。レートに比例して増えます。
2. 濃度分極:$\Delta c_e \sim 570$ mol/m³ × (I/I₁C) の濃度勾配が (DFN-4) の第 2 項を通じて作る電位差 $\dfrac{2RT}{F}(1 - t_+^0)\ln\dfrac{c_e(\text{正極側})}{c_e(\text{負極側})}$。 1C で約 20 mV ですが、$\Delta c_e$ が $c_{e0}$ に近づく 2C 以上では 対数が効いて急速に成長します。 さらに $c_e$ の低下は $\kappa$ と $i_0$ 自体も下げるため、損失は雪だるま式に増えます。

合計すると 1C で 50 mV 程度 — これが上の表の誤差の正体です。 加えて高レートでは $j(x)$ の非一様性(仮定 (ii) の破れ)により、 セパレータ付近の粒子が先に働いて先に疲れる「反応の偏り」が生じますが、 SPM はこれも表現できません。

理解度チェック

次章へ:SPM は「固体の物理だけの電池」でした。 フェーズ3では、捨てた電解液の物質保存 (DFN-2) を復活させて濃度分布 $c_e(x,t)$ を解き、 その分布が作る電位損失を電圧式に足し込む SPMe を導出します。 この章で作った球拡散ソルバはそのまま再利用します。