フェーズ6(任意)拡張 — 熱連成・劣化(SEI)・パラメータ同定

フェーズ5までで、DFN モデルの導出・実装・検証という当初の目標は達成した。 この最終章では「その先」に広がる3つの世界 — 温度との連成、劣化(寿命)のモデリング、 パラメータ同定 — への入口を、それぞれ基礎から丁寧に案内する。 3つとも、DFN を実務や研究で使うなら遅かれ早かれ必ず向き合うテーマである。 完全な導出をこの章だけで済ませることはできないので、 どこまでを自分の手で導出し、どこからが文献参照になるのかを随所で明示する。

6.1 熱連成 — 温度と電池は切り離せない

6.1.1 なぜ温度が重要か

ここまでの章では温度 $T$ [K] を定数(298.15 K)として扱ってきた。 しかし実際のセルでは、放電すれば発熱して温度が上がり、温度が上がれば DFN に登場するほとんどの物性値が変わる。具体的には:

• 輸送係数:固相拡散係数 $D_s$、電解液拡散係数 $D_e$、 イオン伝導率 $\kappa$ はいずれも熱活性化過程であり、温度とともに指数関数的に増大する。 経験的には「10 K 上がると 1.5〜2 倍速くなる」程度の感度を持つ。
• 反応速度:Butler–Volmer 式の反応速度定数 $k$ (交換電流密度 $i_0$ を決める)も同様に熱活性化される。 低温では $i_0$ が小さくなり、同じ電流を流すのに大きな過電圧 $\eta$ が必要になる。
• 安全性:高温では電解液の分解などの発熱副反応が加速する。 発熱が冷却を上回ると温度上昇が止まらなくなる —— これが熱暴走(thermal runaway)である。
• 低温の劣化リスク:低温・大電流充電では負極の反応が追いつかず、 黒鉛への挿入の代わりに金属リチウムが析出(Li plating)しやすくなる。

つまり「電気化学 → 発熱 → 温度 → 物性 → 電気化学」というループが常に回っている。 このループを閉じるのが熱連成(thermal coupling)である。 まずはいちばん簡単で、実務でも最もよく使われる 集中(lumped)熱モデルをエネルギー保存則から導出する。

6.1.2 集中熱モデルの導出

新しく登場する記号を定義する。

6.1.3 発熱項の内訳 — 不可逆熱と可逆熱

式 (6.1) の心臓部は $\dot{Q}_{\mathrm{gen}}$ である。 電池の発熱は大きく不可逆熱(irreversible heat)と 可逆熱(reversible heat, エントロピー熱)に分かれる。 順に導出する。

(a) 不可逆熱:$\dot{Q}_{\mathrm{irr}} = I\,A_{\mathrm{cell}}\,(V_{\mathrm{OCV}} - V)$

まず記号をひとつ定義する。

放電中の端子電圧 $V$ は $V_{\mathrm{OCV}}$ より低い。 その差 $V_{\mathrm{OCV}} - V$ が「損失の合計」に一致することを、 電圧の分解(フェーズ2・3・5 で繰り返し使った考え方)で示す。 簡単のため、各電極を代表的な反応点1点で代表し(SPM と同じ扱い)、 固相の電子伝導によるオーム損は無視する。

(b) 可逆熱:エントロピー項

もうひとつ、損失とは無関係に出入りする熱がある。 電気化学反応そのものが持つエントロピー変化に由来する可逆熱である。 熱力学の基本関係式から1ステップずつ導く。

6.1.4 アレニウス則 — 物性の温度依存

温度が物性に効く経路は、ほとんどの場合アレニウス則(Arrhenius law)で表す。 反応や拡散は「エネルギー障壁 $E_a$ を熱ゆらぎで乗り越える」過程なので、 その速さはボルツマン因子に比例する:

$A_{\mathrm{freq}}$ は頻度因子(前指数因子)[単位は対象量と同じ] である (面積 $A$ と紛らわしいので添字を付けた)。 実用上は頻度因子を消去して、基準温度 $T_{\mathrm{ref}}$ の値からの倍率 で書く形が便利で、DFN 実装でもこの形を使う:

$T > T_{\mathrm{ref}}$ なら指数の中身は正で $\psi$ は増える(速くなる)。 活性化エネルギーの目安(文献値はばらつくので桁感のみ): $D_s$ で 30〜50 kJ/mol、$k$ で 30〜60 kJ/mol、 電解液輸送($D_e$, $\kappa$)で 10〜20 kJ/mol、 SEI 成長などの副反応で 50〜80 kJ/mol。 $E_a = 30$ kJ/mol なら 25 °C → 35 °C で約 1.5 倍、 $E_a = 60$ kJ/mol なら約 2.2 倍になる。

温度と性能のフィードバックループ

アレニウス則があるせいで、温度は一方通行の「出力」ではなく、 ループの一部になる。放電中を考えると:

• 温度が上がる → $k$, $D_s$, $\kappa$ が増える → 過電圧が減る → 不可逆熱 $I A_{\mathrm{cell}}(V_{\mathrm{OCV}}-V)$ が減る ——これは負のフィードバック(自己安定化)。
• 一方、高温では発熱性の副反応(電解液分解など)が アレニウス的に加速する → 発熱が増える ——こちらは正のフィードバック。

発熱率が温度の指数関数で増えるのに対し、冷却率 $h A_{\mathrm{surf}} (T - T_{\mathrm{amb}})$ は温度の1次関数でしか増えない。 両者の交点の有無で運命が分かれる(下の模式図、Semenov 線図と呼ばれる):

冷却線が発熱曲線と交わればそこが安定な動作温度になる(青)。 冷却を弱めていくと2線はやがて接し(緑)、それより弱い冷却では交点が消えて 温度上昇が止まらない(橙)——これが熱暴走の最も単純な描像である。 実際の熱暴走は SEI 分解(約 80〜120 °C)→ セパレータ溶融 → 正極の酸素放出、 と連鎖する多段プロセスで、その定量モデルは文献参照レベルである (総説として Feng et al. 2018 など)。

6.1.5 DFN との連成のしかた

DFN と集中熱モデルの連成は、時間ステップごとに交互に解くのが標準的である (演算子分割 / staggered scheme):

1. 現在の温度 $T_k$ で物性値 $D_s(T_k),\ k(T_k),\ D_e(T_k),\ \kappa(T_k)$ を式 (6.4) で評価する。Butler–Volmer 式の指数の $RT$ も $T_k$ を使う。
2. その物性値で DFN の1ステップ(フェーズ4・5 の方法)を解き、 $V$, $\eta$, $j$ などを得る。
3. 発熱率を計算する:$\dot{Q}_{\mathrm{gen}} = \dot{Q}_{\mathrm{irr}} + \dot{Q}_{\mathrm{rev}}$ (式 (6.2)(6.3)。より厳密には局所発熱の体積積分)。
4. 式 (6.1) を1ステップ積分して $T_{k+1}$ を得る(下のデモでは RK4 を使う)。
5. 1 に戻る。

熱時定数(〜数百秒)は電気化学の最速時間スケール(〜ミリ秒)よりずっと遅いので、 この分割で精度上の問題はまず起きない。 セル内温度分布まで欲しい場合は、集中モデルの代わりに 1D〜3D の熱伝導方程式を解き、 局所発熱密度を右辺に入れる(いわゆる P2D+熱 や 3D 熱連成)。 そのためのパラメータ一式は O'Regan et al. (2022)(LG M50 の熱物性)が整備している。

6.2 劣化 — SEI 成長を例に

6.2.1 劣化モードの分類

「電池が劣化する」という現象は、DFN の言葉では モデルのパラメータや状態が時間とともに変わることとして表現できる。 劣化機構は無数にあるが、電気的な現れ方で整理すると3つに集約される (この分類は Birkl et al. 2017 による整理が標準的):

実セルの容量低下はこれらの合成である。診断には、低レート放電カーブの微分 ($\mathrm{d}V/\mathrm{d}Q$ 解析)などで「どのモードがどれだけ効いたか」を 切り分ける手法が使われる(文献参照レベル)。 本節ではこの中で最も普遍的で、かつ数式で最後まで追える SEI 成長による LLI を掘り下げる。

6.2.2 SEI とは何か

黒鉛負極の動作電位(Li/Li⁺ 基準で 0.05〜0.2 V 程度)は、 電解液(有機カーボネート系)が電気化学的に安定でいられる下限(約 1 V)より低い。 つまり負極表面では電解液の還元分解が熱力学的に常に起こりうる。 最初の充電で分解生成物が負極表面に堆積してできる薄膜が SEI(Solid Electrolyte Interphase, 固体電解質界面被膜)である。

• 良い面:SEI は Li⁺ を通すが電子をほぼ通さないため、 いったん覆ってしまえば電解液のさらなる分解を強く抑制する(不動態化)。 SEI なしには黒鉛負極は成立しない。
• 悪い面:完全には分解を止められず、SEI はゆっくり成長し続ける。 成長のたびに Li と電解液が消費される(→ LLI)、膜が厚くなり抵抗が増える (→ 抵抗増加)。初回充電(化成工程)だけで容量の 5〜10% 程度が SEI 形成に使われる。

6.2.3 成長則の導出 — なぜ膜厚は √t で増えるのか

記号の定義:

6.2.4 サイクル数 vs 容量の典型曲線

サイクル試験では横軸がサイクル数 $N$ になるが、SEI 成長が支配的なら 劣化を進める実質的な変数は「経過時間(と温度・SOC)」なので、 1サイクルあたりの時間 $t_{\mathrm{cyc}}$ を使って $t = N\, t_{\mathrm{cyc}}$ と読み替えれば式 (6.7) がそのまま使える:容量維持率は $1 - \beta\sqrt{N}$ 型の曲線になる($\beta$ は温度・DOD に依存する定数)。 下の図6.2 でこの曲線を操作できる。 実セルでは √N 劣化が数百〜数千サイクル続いたあと、 ニーで急落するパターンが多い(ニーの機構分類は Attia et al. 2022 参照)。

6.2.5 DFN への組み込み方

SEI 副反応を DFN に組み込むには、負極の界面条件を次のように拡張する。 界面を流れる総電流密度を、挿入反応(intercalation)と副反応の和にする:

それぞれの構成則は:

ここで両方の過電圧に SEI 膜のオーム損を含める: $\eta = \phi_s - \phi_e - U(\theta) - j_{\mathrm{tot}} R_{\mathrm{film}}$、 $\eta_{\mathrm{SEI}} = \phi_s - \phi_e - U_{\mathrm{SEI}} - j_{\mathrm{tot}} R_{\mathrm{film}}$。

そのうえで、各保存則の $j$ を使い分ける:

• 固相拡散の境界条件(粒子に入る Li): $-D_s \partial c_s/\partial r|_{r=R_s} = j_{\mathrm{int}}/F$ —— 粒子に挿入されるのは挿入反応の分だけ。
• 電荷保存($\phi_s$, $\phi_e$ の式):$a j \to a j_{\mathrm{tot}}$ —— 電子とイオンは副反応にも流れる。
• SEI の状態方程式:式 (6.5) で $\delta(x,t)$ を時間発展させ、 消費された Li を LLI として集計する。

充電時は $\phi_s - \phi_e$ が下がるので $|j_{\mathrm{SEI}}|$ が増える —— 「満充電付近での保存・充電で劣化が速い」という経験則が式から出てくる。

6.3 パラメータ同定 — DFN 実用最大の壁

6.3.1 何が測れて、何が難しいか

DFN の方程式系はフェーズ5で完成した。しかし実務で最初につまずくのは数値解法ではなく 「この 30 個近いパラメータをどこから持ってくるか」である。 パラメータは測定の難易度で明確に階層化される:

「難しい」側のパラメータは、結局セルの電圧応答にモデルを合わせ込んで推定 することになる。これがパラメータ同定(parameter identification / estimation)である。

6.3.2 最小二乗フィッティングの定式化

推定したいパラメータをまとめてベクトル $\boldsymbol{p} = (p_1, \dots, p_m)$ と書く(例:$\boldsymbol{p} = (D_s^+, k^+)$)。 実験で時刻 $t_i$($i = 1,\dots,M$)の端子電圧 $V_{\mathrm{exp}}(t_i)$ を測り、同じ電流プロトコルでモデル電圧 $V_{\mathrm{model}}(t_i; \boldsymbol{p})$ を計算する。 残差 $r_i(\boldsymbol{p}) = V_{\mathrm{model}}(t_i;\boldsymbol{p}) - V_{\mathrm{exp}}(t_i)$ の二乗和を最小にする $\boldsymbol{p}$ を探すのが最小二乗法である:

実装上の定石が3つある。 (1) $D_s$ や $k$ のように桁で動く正のパラメータは $\log_{10} p$ を変数にする(スケールが揃い、正値制約も自動で満たされる)。 (2) 物理的にありうる範囲で上下限を制約する。 (3) 最適化には残差ベクトルの構造を使う Gauss–Newton 系 (Levenberg–Marquardt、scipy では least_squares)を使う。 $J$ は一般に非凸なので、初期値によっては局所解に落ちる —— これは後の実行例で実際に確認できる。

6.3.3 感度と識別可能性 — フィットすれば決まる、とは限らない

最小二乗が「うまくいく」ためには、そもそも データがそのパラメータの情報を含んでいる必要がある。 それを測る道具が感度(sensitivity)である:

感度が小さいパラメータはノイズに埋もれて決まらない。 しかしもっと厄介なのは、感度は大きいのに他のパラメータと区別が付かない場合で、 これを識別可能性(identifiability)の欠如と呼ぶ。 2つのパラメータの感度曲線 $S_i(t)$ と $S_j(t)$ が時間軸上で「同じ形」をしていると、 $p_i$ を増やした効果を $p_j$ を減らして打ち消せるため、 コスト関数 $J$ の地形に細長い「谷」ができ、谷底のどこでも $J$ がほぼ同じになる。 形の類似度はコサイン類似度(相関)で定量化できる:

具体例1(構造的な識別不能):$D_s$ と $R_s$ は積でしか効かない

SPM の固相拡散方程式を無次元化してみる。 $\tilde{r} = r/R_s$、$\tilde{t} = t/\tau_d$($\tau_d$ は後で決める)とおいて $\partial c_s/\partial t = \frac{1}{r^2}\partial_r (D_s r^2 \partial_r c_s)$ に代入すると:

境界条件の流束も $j R_s / (F D_s)$ という組合せでしか入らない(SPM では $j \propto 1/(aL) \propto R_s$ なので、これも $R_s^2/D_s$ に帰着する)。 つまり $V(t)$ には $D_s$ と $R_s$ が拡散時定数 $\tau_d = R_s^2/D_s$ という 1つの組合せでしか現れない。 $R_s$ を4倍にして $D_s$ を16倍にした電池は、電圧応答が完全に同じ —— どれだけ良いデータを何回取っても、$V(t)$ だけから両者を別々に決めることは 原理的に不可能である(構造的識別不能)。 だから $R_s$ は粒度分布で独立に測り、フィットは $D_s$ に絞る、という分業になる。

具体例2(実用上の識別困難):$k$ と接触抵抗

定電流放電では界面電流密度 $j$ がほぼ一定なので、反応過電圧 $\eta \approx (2RT/F)\,\mathrm{asinh}(j/2i_0)$ は時間的にほぼ一定の電圧オフセットになる。 ところが接触抵抗 $R_{\mathrm{contact}}$ による電圧降下 $I R_{\mathrm{contact}}$ も「一定のオフセット」。2つの感度曲線はほぼ平行($\rho \approx 1$)で、 単一の定電流カーブからは事実上分離できない。 分離するには、$i_0$ の非線形性が顔を出す電流を変えた実験 (マルチレート、パルス、EIS)が必要になる。 —— 「識別可能性は実験デザインの問題」というのが本節で一番大事な教訓である。

6.3.4 GITT と EIS(原理の紹介のみ)

フィット一辺倒を避けるため、個別パラメータを狙い撃ちする測定法が使われる。 2つだけ原理を紹介する(定式化の導出は文献参照レベル)。

• GITT(Galvanostatic Intermittent Titration Technique): 小さな定電流パルスを打っては電流を切り、緩和させる、を繰り返す。 パルス中の電圧変化のうち「√t で進む成分」の傾きが固相拡散の速さを反映し、 緩和後の電圧差が OCP の傾きを与える。両者を組み合わせると $D_s$ が $\theta$ の関数として求まる(Weppner & Huggins 1977)。 フェーズ2で作った SPM は、GITT の応答を再現・検証するのにそのまま使える。
• EIS(電気化学インピーダンス分光): 微小な交流電圧(電流)を周波数を変えながら印加し、応答の複素インピーダンス $Z(\omega)$ を測る。プロセスごとに時定数が違うことを利用して、 kHz 帯のオーム抵抗・SEI 抵抗、Hz 帯の電荷移動抵抗(→ $i_0$、つまり $k$)、 mHz 帯の拡散(Warburg 応答、45° の直線)を周波数軸上で分離できる。 DFN を線形化すればモデル側のインピーダンスも計算でき、EIS を使った同定もできる。

6.3.5 過学習と検証データの分離

パラメータを増やせば $J$ はいくらでも下がるが、それは ノイズや系統誤差まで暗記した(過学習)だけかもしれない。 機械学習と同じ規律が要る:

• 訓練と検証を分ける:例えば 0.5C と 2C の放電カーブでフィットし、 1C とパルス応答で検証する。検証誤差が訓練誤差より大幅に悪ければ過学習を疑う。
• プロトコルをまたいで検証する:定電流だけで合わせたモデルが パルスやドライブサイクルでも合うかを見る。合わなければ、 パラメータではなくモデル構造が間違っている可能性が高い。
• 物理的な妥当性チェック:推定値が文献レンジ (例:黒鉛の $D_s$ で $10^{-15}$〜$10^{-13}$ m²/s)から桁で外れたら、 たとえ $J$ が小さくても採用しない。パラメータが「他の誤差の吸収係数」に なっている兆候である。

6.3.6 やってみる:合成データからの最小二乗同定(この章の目玉)

仕上げに、パラメータ同定の全工程をブラウザ内で実行する。 「真のパラメータ」で合成データを作る → ノイズを載せる → scipy.optimize.least_squares で推定 → 真値と比較、という流れは、 実データで同定する前に必ずやるべきリハーサルでもある (合成データで真値を回収できない手順が、実データで機能するはずがない)。 モデルは図6.3 と同じ簡易 SPM の Python 版で、 負極は既知とし、正極の $D_s^+$ と $k^+$ の2つを推定する。

# ============================================================
# パラメータ同定のミニ実験:
#  (1) 「真の」D_s^+ と k^+ で簡易 SPM の 1C 放電カーブを合成し、
#      観測ノイズを加えて「実験データ」とする
#  (2) scipy.optimize.least_squares で 2 パラメータを推定する
#  (3) 真値・初期値・推定値と残差を比較する
# ============================================================
import numpy as np
from scipy.optimize import least_squares

F, Rg, T = 96485.33212, 8.314462618, 298.15
# --- LG M50(Chen 2020)のうち SPM に必要な値 ---
L_n, L_p, R_n, R_p = 85.2e-6, 75.6e-6, 5.86e-6, 5.22e-6
a_n, a_p = 3*0.75/R_n, 3*0.665/R_p           # 比界面積 [1/m]
csmax_n, csmax_p = 33133.0, 63104.0
Ds_n, k_n = 3.3e-14, 6.716e-12               # 負極は既知とする
ce0, I_1C = 1000.0, 5.0/0.1027               # 1C 電流密度 [A/m^2]
th_n100, th_p100 = 0.9014, 0.27              # 100% SOC の化学量論

def ocp_n(t):   # 黒鉛 OCP [V]
    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 ocp_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 surf_conc(D, Rs, c0, jF, dt, nstep, N=14):
    """球粒子拡散(FVM + Crank-Nicolson)。表面濃度の時系列を返す。
    jF = j/F [mol/(m^2 s)] は表面から出る向きが正(一定)。"""
    dr = Rs/N
    rf = np.arange(N+1)*dr                       # セル境界の半径
    vol = (rf[1:]**3 - rf[:-1]**3)/3.0           # セル体積(4π は省略)
    g = D*rf[1:-1]**2/dr                         # 境界面コンダクタンス
    Lm = np.zeros((N, N))                        # 拡散演算子 L
    for i in range(N-1):
        Lm[i, i]     -= g[i]/vol[i];   Lm[i, i+1] += g[i]/vol[i]
        Lm[i+1, i+1] -= g[i]/vol[i+1]; Lm[i+1, i] += g[i]/vol[i+1]
    src = np.zeros(N); src[-1] = -jF*rf[-1]**2/vol[-1]   # 表面流出
    A = np.eye(N) - dt/2*Lm
    M = np.linalg.solve(A, np.eye(N) + dt/2*Lm)  # CN の 1 歩: c ← M c + s
    s = np.linalg.solve(A, dt*src)
    c = np.full(N, c0); out = np.empty(nstep)
    for k in range(nstep):
        c = M @ c + s
        out[k] = c[-1] - jF*(dr/2)/D             # セル中心値 → 表面値に補正
    return out

dt, nstep = 30.0, 100                            # 0〜3000 s(1C 放電)
t_grid = dt*np.arange(1, nstep+1)
I = 1.0*I_1C
j_n, j_p = I/(a_n*L_n), -I/(a_p*L_p)             # 界面電流密度 [A/m^2]

# 負極は推定対象でないので 1 回だけ計算してキャッシュ(高速化のコツ)
cs_n = np.clip(surf_conc(Ds_n, R_n, th_n100*csmax_n, j_n/F, dt, nstep),
               1.0, csmax_n-1.0)
i0_n = k_n*F*np.sqrt(ce0*cs_n*(csmax_n-cs_n))
V_neg = ocp_n(cs_n/csmax_n) + 2*Rg*T/F*np.arcsinh(j_n/(2*i0_n))

def model_V(p):
    """p = [log10 D_s^+, log10 k^+] → 端子電圧 V(t)(対数で扱うと安定)"""
    Dsp, kp = 10.0**p[0], 10.0**p[1]
    cs_p = np.clip(surf_conc(Dsp, R_p, th_p100*csmax_p, j_p/F, dt, nstep),
                   1.0, 0.9999*csmax_p)
    i0_p = kp*F*np.sqrt(ce0*cs_p*(csmax_p-cs_p))
    eta_p = 2*Rg*T/F*np.arcsinh(j_p/(2*i0_p))
    return ocp_p(cs_p/csmax_p) + eta_p - V_neg

# --- (1) 合成データ:真値 + 観測ノイズ(σ = 5 mV)---
p_true = np.log10([4.0e-15, 3.545e-11])
rng = np.random.default_rng(0)
V_data = model_V(p_true) + rng.normal(0.0, 0.005, nstep)

# --- (2) 最小二乗推定(log10 空間、物理的にありうる範囲を制約)---
p0 = np.log10([1.0e-15, 1.0e-10])                # 初期値(真値からずらす)
bounds = ([-17.0, -13.0], [-12.0, -8.0])         # 10^-17〜10^-12 など
fit = least_squares(lambda p: model_V(p) - V_data, p0, bounds=bounds)

# --- (3) 結果の比較 ---
print("パラメータ     真値        初期値      推定値      誤差")
for i, nm in enumerate(["D_s^+ [m2/s]", "k^+   [SI]  "]):
    tv, ev = 10**p_true[i], 10**fit.x[i]
    print(f"{nm}  {tv:.3e}  {10**p0[i]:.3e}  {ev:.3e}  {100*(ev-tv)/tv:+.1f} %")
print(f"残差 RMS = {np.sqrt(np.mean(fit.fun**2))*1000:.2f} mV(ノイズ σ = 5 mV)")
print("残差がノイズ水準まで落ちれば、データはモデルで説明できたことになる。")
print("注意: 初期値を大きく外す(例: D_s^+ = 1e-17)と局所解に落ちることがある。")
print("      また k^+ の誤差が D_s^+ より大きいのは、感度が小さい(図6.3)ため。")

plot_spec = {
  "traces": [
    {"x": t_grid[::3].tolist(), "y": V_data[::3].tolist(),
     "name": "合成データ(ノイズ付き)", "mode": "markers",
     "marker": {"size": 5, "color": "#9ca3af"}},
    {"x": t_grid.tolist(), "y": model_V(p0).tolist(),
     "name": "初期値のモデル", "mode": "lines",
     "line": {"dash": "dot", "color": "#f59e0b"}},
    {"x": t_grid.tolist(), "y": model_V(fit.x).tolist(),
     "name": "推定値のモデル", "mode": "lines", "line": {"color": "#2563eb"}},
  ],
  "layout": {"title": {"text": "最小二乗フィット:1C 放電カーブ"},
             "xaxis": {"title": {"text": "時間 [s]"}},
             "yaxis": {"title": {"text": "端子電圧 [V]"}}},
}

実行すると、$D_s^+$ は数 % の誤差で回収できるのに対し、$k^+$ の誤差は 10% 近くなる —— 図6.3 で見た感度の大きさの差が、そのまま推定精度の差になって現れる。 さらに p0 を np.log10([1.0e-17, 1.0e-10]) に書き換えて 再実行してみてほしい。残差 400 mV 超の局所解(拡散が遅すぎて表面が飽和し、 どのパラメータ方向にも改善しない停留点)に捕まるはずだ。 実務では複数の初期値から出発して同じ解に収束するかを確認する (マルチスタート)のが定石である。

6.4 さらに学ぶには

6.4.1 PyBaMM エコシステム

本サイトでは学習のためにすべてを自作したが、研究・実務では PyBaMM (Python Battery Mathematical Modelling; Sulzer et al. 2021)を使うのが近道である。 フェーズ5の検証データも PyBaMM で生成した。この章で学んだ3テーマは すべて PyBaMM に実装済みで、数行で試せる:

• 実験プロトコル:pybamm.Experiment に 「1C で 2.5 V まで放電 → CCCV 充電 → 休止」のような文字列を渡すだけで、 サイクル試験や GITT を再現できる。
• 熱モデル:"thermal": "lumped" オプションで 本章 6.1 の集中熱モデルが、"x-full" で 1D 熱伝導が有効になる。
• 劣化モデル群:SEI(反応律速/溶媒拡散律速/電子トンネリングなど 複数の律速仮説)、Li 析出、粒子割れ、LAM が選択式で入っており、 組み合わせて長期サイクルを回せる(O'Kane et al. 2022)。
• パラメータセット:本サイトで使った Chen2020(LG M50)のほか、 熱物性込みの ORegan2022、劣化パラメータ込みの OKane2022 などが同梱されている。
• liionpack:PyBaMM のセルモデルを直並列に接続して モジュール/パックのシミュレーションを行う拡張ライブラリ。

# PyBaMM での「本章の内容全部入り」の例(ローカル環境で: pip install pybamm)
# 熱連成 + SEI 劣化 + サイクル実験を数行で。
import pybamm

model = pybamm.lithium_ion.DFN({
    "SEI": "solvent-diffusion limited",   # 6.2 の √t 型 SEI モデル
    "thermal": "lumped",                  # 6.1 の集中熱モデル
})
experiment = pybamm.Experiment([
    ("Discharge at 1C until 2.5 V",
     "Charge at 0.3C until 4.2 V",
     "Hold at 4.2 V until C/50"),
] * 100)                                  # 100 サイクル
sim = pybamm.Simulation(model, experiment=experiment,
                        parameter_values=pybamm.ParameterValues("OKane2022"))
sol = sim.solve()
pybamm.plot_summary_variables(sol)        # 容量・LLI・膜厚などの推移を一括表示

6.4.2 主要文献ガイド

次に読むべきものを、目的別に挙げる。

6.4.3 本サイトの総まとめ — フェーズ0〜6 で積み上げたもの

章末:理解度チェック

学習の完了と今後の道筋

ここまで来たあなたは、フェーズ0の「電池の中では何が起きているのか」という問いから始めて、 保存則から DFN を導出し(1)、SPM・SPMe という階層を自分の手で実装し(2・3)、 非線形連立系を解く数値技術を身につけ(4)、フルの DFN を組み上げて 公開実装と突き合わせ(5)、そして本章で温度・劣化・同定という 「現実の電池工学」への3本の橋を渡り始めた。 DFN モデルはもはやブラックボックスではない —— どの式がどの仮定から来て、どのパラメータがどう効き、 どこまで信じてよいかを、あなたは自分で判断できる。

これからの進み方の提案:

1. PyBaMM を触る:上のコード例から始めて、 本サイトの各フェーズの結果を PyBaMM で再現してみる。 自作コードとの差分を説明できたら、理解は本物である。
2. 論文を1本選んで実装する:6.4.2 の文献から (例えば SEI モデルの Safari 2009 を)選び、式を追い、 フェーズ2の SPM に組み込んでみる。本サイトでやった 「導出 → 実装 → 検証」のループはそのまま研究の作法である。
3. 実データで同定してみる:手に入る充放電データ (公開データセットでもよい)に 6.3 の手順を適用する。 合成データとの落差 —— モデル構造誤差、ドリフト、ヒステリシス —— こそが次に学ぶべきテーマを教えてくれる。
4. 興味の方向で分岐する:BMS・状態推定なら Plett、 材料・界面なら Newman & Thomas-Alyea、 安全性・熱なら Bernardi と熱暴走の文献群へ。 どの道でも、本サイトで積んだ「保存則と界面反応から出発する」思考法が土台になる。

長い道のりに付き合ってくれてありがとう。よい電池モデリングを。