オイラー法で非等温反応を解くだけ

背景

　Oiita自体に投稿するのが初めてなので最近始めたGitの練習と合わせて投稿してみるだけ、あとMarkdown記法の練習もかねて。ついでに言うなら化学系の投稿が少ないので少し寂しいので。
　Gitのリポジトリはここ

計算手法

　オイラー法を説明する必要は全くないと思うので、実際の設計方程式を変形した形で下に示す。

$$\frac{dx_A}{dt}=\frac{r_A}{C_{A0}}=k_oexp(-E/RT)(1-x_A)$$

$$\frac{dT}{dt}=\frac{r_AΔH}{ρc_{pm}}-\frac{UA(T_s-T)}{Vρc_{pm}}$$

　上の2式をオイラー法で解く(上：設計方程式、下：熱収支式)。上の反応速度に注目すれば分かるが、反応速度は温度に対して指数関数的に変化する。なので等温反応を仮定して解くと実際の反応時間と反応率に対して大きく違う答えが出てくる。

練習問題

　培風館の『反応工学』の例題を取り上げる。

A→Cで表される液相反応を非等温の回分反応器を用いて行う。すなわち、内径D=0.8 m、高さ1.3 mの撹拌僧にAのみからなる原料を500 ㎏仕込んで、槽の側面と底部をジャケットでつつみ、その中に熱媒体を流して、槽外壁を常に613 Kに保って反応させる。反応率が95%になる反応時間を求めよ。本反応は1次反応であって、反応速度は次式で与えられる。

-ra=3.228×10^13×exp(-E/RT)×Ca

ただし、E=186.2×10^3 J/mol
反応熱 ΔH=62760 J/mol
反応液密度 ρ=900 ㎏/m^3
反応液の平均比熱容量 Cp=2.51×10^3 J/kg/K
Aの分子量 M=385×10^-3 kg/mol
総括伝熱係数 U=523 W/m^2/K

培風館『反応工学』 p170 例題7・4

解法

　基本的にオイラー法では$a<b$に対して区間$[a,b]$を何ステップに分けるかで決めるが、ここではstepで定めた時間数を終了条件まで計算し続ける。反応速度式の計算ならそこまで問題にならない。気になるなら微分方程式を反転させて(逆数を取って)終了条件の$x_A$に対して区間で区切ればいい。

　言語はPython、環境はJupyter notebookで行った。コードはここ。化学系では実行速度や過去の遺産によりいまだにFORTRANやCが活用されているがこれくらいはpythonで問題がない。むしろグラフ描画がやりやすい分、個人的には好みである。設計の段ではExcel+VBAというのも計算とグラフ、表計算の3点セットで使いやすい。

　反応率$x_A=0.95$で終了。$step=0.1 s$で19988回。2000秒で反応終了。これくらいならすぐ。
　反応率の変化と温度変化を下に示す。

　温度の単位はK。613K(330℃)から反応がスタートしたが途中で吸熱反応により温度が590Kまで低下しているのが分かる。

等温と非等温

　例題はここまでだが、実際に温度変化がどれくらい反応に寄与しているのかを計算してみる。つまり熱収支式を無視して613K一定で設計方程式を解く。

　そのあと613Kではなく100℃低い513Kで計算をして結果を見る。コードはここ

　613Kでの等温の結果は以下のようになった。

　オレンジの線が等温での計算結果である。20K程度温度が下がっただけで反応時間が3倍程度になることが分かる。反応速度がいかに温度変化に敏感かよく分かる。

　次に513Kでの非等温反応を解いてみる。が、解いてみた温度変化の式が下になる。

　反応を通して513K一定であることが分かる。100℃低いと反応速度が極端に遅くなるので反応熱より壁面からの熱伝導のほうが圧倒的に大きくなる。

　今度は同じく513Kでの等温ー非等温を比べてみる。

　完全に等温反応と非等温の反応が一致している。

　それにしても613Kで2000sの反応が513Kで800000s、ほとんど9日かかる。この場合、温度が高いほうが有利である。
　しかし、今回は吸熱反応で温度が下がる方向にはたらいたが、発熱反応で温度が上がる場合、上に見たように反応速度が急激に上昇しそれがさらに温度を上げるので熱暴走の危険がある。その辺は熱的安定性の問題になるが、温度が簡単に上昇するような状況は良くない。