モデリングと数値計算（主に化学工学向け）

化学反応、移動現象のモデリングと数値計算の話

バッチ反応のモデリングと数値計算

はじめに

バッチ反応(液相、一次反応)の非常に簡単な例でモデリング、数値計算を行う。

モデリング

下図に示すとおり成分Aから成分Bを生じる液相バッチ反応について、時刻tにおける成分Aの量を求めることを考える。

f:id:schemer1341:20200709232018p:plain:w300

物質収支

(成分A物質量の変化速度) = (成分Aの反応速度)であり、時刻をt [s]として次式で表される。

$\displaystyle{\frac{dn_A}{dt}=r_AV}$

両辺を定数である体積Vで割り、 $n_A/V=c_A$ (濃度)とすると

$\displaystyle{\frac{dc_A}{dt}=r_A}$

反応速度

成分Aの反応速度は成分Aの濃度の１次に比例とすると、反応速度定数を $k$ [1/s]として次式で表される。

$-r_A = kc_A$

(Aは消失するのでマイナスが付く)

解析解

物質収支と反応速度から次の形となる。

$\displaystyle{\frac{dc_A}{dt}=-kc_A}$

ここで初期条件を $c_A=c_{A0}$ at $t=0$ として方程式を解くと $c_A$ と $t$ の関係が求まる。

$\begin{aligned} \frac{1}{c_A}dc_A &= -kdt \\ \int^{c_A}_{c_{A0}}\frac{1}{c_A}dc_A &= -k\int^{t}_{0}dt \\ \ln \frac{c_A}{c_{A0}} &= -kt \\ c_A &= c_{A0} e^{-kt} \end{aligned}$

数値解

オイラー法を用いると時刻 $t$ における $c_A|_t$ とその微小時間 $\Delta t$ 後の $c_A|_{t+\Delta t}$ は次式で関連付けられる。

$\displaystyle{c_A|_{t+\Delta t} = c_A|_t + \frac{dc_A|_t}{dt} \Delta t}$

微分部分を代入すると、

$\displaystyle{c_A|_{t+\Delta t} = c_A|_t - k c_A|_t \Delta t}$

具体的な微小時間値を決めたら初期条件から順次数値を計算していけばよい。

解析解と数値解の比較

Pythonを用いて計算した。ソースコードは下記リンク先にある。
cemodeling/batch_reaction.py at master · nakamura-13/cemodeling · GitHub

初期濃度 $c_{A0}=1$ [mol/m3]、反応速度定数 $k=0.01$ [1/s]、微小時間を30秒として計算した結果を下グラフに示す。横軸は時刻、縦軸は濃度である。600秒くらいで成分Aは反応してほぼ無くなる。青い線が数値解である。

f:id:schemer1341:20200710223918p:plain:w500

微小時間を1秒にすると、解析解と数値解はほぼ一致する。

f:id:schemer1341:20200710223957p:plain:w500

おわりに

ここでは常微分方程式の問題を取り扱ったが、次回は偏微分方程式の問題に取り組み予定。