目次
はじめに
MathWorksのSimscape Batteryという製品にBattery Single Particleというブロックがある。 これはバッテリーの電気化学モデルのブロックで、内部では固相のSingle Particle modelと液相の方程式を解いている。
ヘルプページに固相の拡散方程式を差分方程式に変換した結果がある。
ここでは、この差分方程式を導出をしてみる。
なお、下記では電流の符号は放電を正として、$c_s$が正極のリチウムイオン濃度とすることを想定している。 負極の場合は電流の符号は逆にする必要がある。
電流の符号はそれぞれ流儀があるが、バッテリーを電圧源と考える場合は放電を正にしておくのが回路要素の機能としては一貫性が保たれる。
解きたい拡散方程式
Single Particle modelで解く拡散方程式はこれ:
\[ \cfrac{\partial c_s(t, r)}{\partial t} = D_s \cfrac{1}{r^2} \cfrac{\partial}{\partial r} r^2 \cfrac{\partial c_s(t, r)}{\partial r} \]差分方程式にするには、微分の中の積を展開する必要がある。
$r$ に関する$r^2 \partial c_s/\partial r$の微分は
\[ 2r c_s(t, r) + r^2\cfrac{\partial^2 c_s(t, r)}{\partial r^2} \] だから、上の拡散方程式は最終的に \[ \cfrac{\partial c_s(t, r)}{\partial t} = D_s \left[ \cfrac{2}{r} \cfrac{\partial c_s(t, r)}{\partial r} + \cfrac{\partial^2 c_s(t, r)}{\partial r^2} \right] \]となる。($c_s$は活物質内のリチウムイオン濃度、$D_s$は拡散定数、$r$は活物質を球体とした場合の中心からの位置座標である。)
これを差分方程式にする。実際には下で実施する。
境界条件
境界条件は二つある。
$r=0$ での境界条件
\[ \left. D_s \cfrac{\partial c_s}{\partial r} \right|_{r=0} = 0. \]実のところ、差分方程式では$r=0$での境界条件は考慮されないことになる。これについては後で説明する。
$r=R_s$ での境界条件
\[ \left. D_s \cfrac{\partial c_s}{\partial r} \right|_{r=R_s} = \cfrac{I}{S F}. \]ヘルプ・ページの表記は少し混乱しているように思われる。$I$は電流(電流密度ではない)、$S$は電気化学的に活性な表面積で
\[ S = \cfrac{3\varepsilon V}{R_s} \]である。$V$は電極の体積、すなわち電極の面積と電極の厚さの積である。また$\varepsilon$はvolume fractionである (volume fractionの説明は省く。)
なお、右辺をファラデー定数$F$で割っているのは、電流密度からモル密度に変換するためである。
差分方程式
まず$r$の微小量を定義しておく:
\[ \delta r \triangleq \cfrac{R_s}{N} \]MathWorksのヘルプ・ページでは小文字の$n$を使っているが、$n$はindex的な使われ方が多いので、ここでは大文字の$N$で全分割数を表すことにした。
こうすると
\[ r_i = i \delta r \]となる。
ポイント
ここでのポイントは1階の微分を$i$について対称な形にすることだ。
1階の微分を無邪気に差分化すると
\[ \cfrac{\partial c_s}{\partial r} = \cfrac{c_{s, i+1} – c_{s, i} }{\delta r} \]あるいは
\[ \cfrac{\partial c_s}{\partial r} = \cfrac{c_{s, i} – c_{s, i-1} }{\delta r} \]になる。これらは$i$について対称ではない。
これを対称にするには
\[ \cfrac{\partial c_s}{\partial r} = \cfrac{c_{s, i+1} – c_{s, i-1} }{2\delta r} \]とすれば良い。
2階の微分の差分方程式
2階の微分を差分方程式にするのは、よくやるやり方で問題ない:
\[ \cfrac{\partial^2 c_{s, i}}{\partial r^2} = \cfrac{ c_{s, i-1} -2 c_{s, i} +c_{s, i+1}}{\delta r^2}. \]差分方程式
これで準備はできた。
$$ \begin{aligned} \cfrac{\partial c_{s, i}}{\partial t} = & D_s \left[ \cfrac{2}{i \delta r} \cfrac{c_{s, i+1} – c_{s, i-1} }{2 \delta r} + \cfrac{ c_{s, i-1} -2 c_{s, i} +c_{s, i+1}}{\delta r^2} \right] \\ = & \cfrac{D_s}{\delta r ^2} \left[ \cfrac{c_{s, i+1} – c_{s, i-1} }{i } + \cfrac{i\left( c_{s, i-1} -2 c_{s, i} +c_{s, i+1} \right) }{i } \right] \end{aligned} $$最終的に極座標系の拡散方程式の差分方程式は以下のようになる:
\[ \cfrac{\partial c_{s, i}}{\partial t} = \cfrac{D_s}{\delta r ^2} \left[ \cfrac{i-1}{i} c_{s, i-1} – 2 c_{s, i} + \cfrac{i+1}{i}c_{s, i+1} \right] \]境界条件
$r=0$ ($i=1$) での境界条件
上の差分方程式で$i=1$とすれば
\[ \begin{aligned} \cfrac{\partial c_{s, 1}}{\partial t} = & \cfrac{D_s}{\delta r ^2} \left[ \cfrac{0}{1} c_{s, 0} – 2 c_{s, 1} + \cfrac{2}{1}c_{s, 2} \right] \\ = & \cfrac{D_s}{\delta r ^2} \left[ – 2 c_{s, 1} + 2c_{s, 2} \right]\end{aligned} \]となる。
上で書いた
\[ \left. D_s \cfrac{\partial c_s}{\partial r} \right|_{r=0} = 0. \]は考慮されない。これで良いのか、自分は分かっていない。
$r=R_s$ ($i=N$) での境界条件
\[ \left. D_s \cfrac{\partial c_s}{\partial r} \right|_{r=R_s} = \cfrac{I}{S F}. \]を差分方程式にする。
$i$に関して対称形にすることに注意して、$R_s = N \delta r$であるから、
\[ D_s \cfrac{ c_{s, N+1} – c_{s, N-1}}{2\delta r} = \cfrac{I}{S F} \]となる。書き直して
\[ c_{s, N+1} = c_{s, N-1} + \cfrac{2\delta r}{D_s}\cfrac{I}{S F} \]これが境界条件の差分方程式である。
拡散方程式の$r=R_s$に相当する部分は \[ \cfrac{\partial c_{s, N}}{\partial t} = \cfrac{D_s}{\delta r ^2} \left[ \cfrac{N-1}{N} c_{s, N-1} – 2 c_{s, N} + \cfrac{N+1}{N}c_{s, N+1} \right] \]である。ここで$N+1$の項が含まれているから、これを境界条件を使って消去する。
\[ \cfrac{\partial c_{s, N}}{\partial t} = \cfrac{D_s}{\delta r ^2} \left[ \cfrac{N-1}{N} c_{s, N-1} – 2 c_{s, N} + \cfrac{N+1}{N} \left( c_{s, N-1} + \cfrac{2\delta r}{D_s}\cfrac{I}{S F}\right) . \right] \]少しだけ整理すれば、これは最終的に
\[ \cfrac{\partial c_{s, N}}{\partial t} = \cfrac{D_s}{\delta r ^2} \left[ 2 c_{s, N-1} – 2 c_{s, N} \right] + \cfrac{N+1}{N} \cfrac{2I}{S F \delta r} . \]となる。
まとめておくと
拡散方程式そのものは
\[ \cfrac{\partial c_{s, i}}{\partial t} = \cfrac{D_s}{\delta r ^2} \left[ \cfrac{i-1}{i} c_{s, i-1} – 2 c_{s, i} + \cfrac{i+1}{i}c_{s, i+1} \right]. \]$i=1$では
\[ \cfrac{\partial c_{s, 1}}{\partial t} = \cfrac{D_s}{\delta r ^2} \left[ – 2 c_{s, 1} + 2c_{s, 2} \right] \]$i=N$では
\[ \cfrac{\partial c_{s, N}}{\partial t} = \cfrac{D_s}{\delta r ^2} \left[ 2 c_{s, N-1} – 2 c_{s, N} \right] + \cfrac{N+1}{N} \cfrac{2I}{S F \delta r} . \]これで$c_s$ベクトルの微分方程式を数値積分する形に書き直すことができた。
具体的な例
上の差分方程式だけでは、実際にどのように数値積分するのかイメージが湧かない人のために$N=5$の例を書いておく。
\[ \cfrac{\mathrm{d}}{\mathrm{d}t} \begin{bmatrix} c_{s,1} \\ c_{s,2} \\ c_{s,3} \\ c_{s,4} \\ c_{s,5} \end{bmatrix} = \cfrac{2D_s}{\delta r^2} \begin{bmatrix} -2 & 2 & 0 & 0 & 0 \\ 1/2 & -2 & 3/2 & 0 & 0 \\ 0 & 2/3 & -2 & 4/3 & 0 \\ 0 & 0 & 3/4 & -2 & 5/4 \\ 0 & 0 & 0 & 2 & -2 \\ \end{bmatrix} \begin{bmatrix} c_{s,1} \\ c_{s,2} \\ c_{s,3} \\ c_{s,4} \\ c_{s,5} \end{bmatrix} + \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ \cfrac{6}{5}\cfrac{2I}{S F \delta r} \end{bmatrix} \]ここまで来れば、ベクトル$[c_{s,1}, c_{s,2}, c_{s,3}, c_{s,4}, c_{s,5}]^T$を積分するイメージが掴めるだろう。 意外と簡単に数値積分できる形になるし、積分も容易だ。
