管型反応器のモデリングと数値計算 - モデリングと数値計算（主に化学工学向け）

はじめに

シンプルな条件での管型反応器のモデリング・数値計算を行う。

管型反応器

下図のようなA→Bの反応を行う管型反応器において、初期は管型反応器内に成分A, Bが存在しない状態として、入口から成分Aを供給し始めたあとの反応器内のAの濃度変化を考える。反応熱は小さく等温反応であり、また流体の流れはプラグフローであり管半径方向の流量、濃度は一定であるものとする。

f:id:schemer1341:20200724002116p:plain:w350

物質収支

上図の微小区間⊿z [m]における成分Aの物質収支は次の関係で表される。

(Aの流入) - (Aの流出) + (反応によるAの生成) = (Aの蓄積)

Aの流量を $F_A$ [mol/s]、Aの反応速度を $r_A$ [mol/m3/s]、微小区間内の成分Aの物質量を $n_A$ [mol]として微小時間⊿t [s]における物質収支を数式で表すと下記となる。

$\displaystyle{ F_A|_z - F_A|_{z+\Delta z} + r_A \times S \times \Delta z = \frac{n_A|_{t+\Delta t} - n_A|_t}{\Delta t} }$

$S$ は反応管断面積 [m2] である。ここで両辺を $S\Delta z$ で割り、 $n_A/(S\Delta z)=c_A$ (濃度 [mol/m3])とすると、

$\displaystyle{ \frac{1}{S} \frac{F_A|_z - F_A|_{z+\Delta z}}{\Delta z} + r_A = \frac{c_A|_{t+\Delta t} - c_A|_t}{\Delta t} }$

$\Delta z \rightarrow 0$ 、 $\Delta t \rightarrow 0$ の極限を取ると

$\displaystyle{ -\frac{1}{S} \frac{\partial F_A}{\partial z} + r_A = \frac{\partial c_A}{\partial t} }$

ここで $F_A = vc_A$ ( $v$ は全流体流量 [m3/s])で表され、 $v$ は一定であると仮定し、空塔速度 $u=v/S$ [m/s]を用いると次式となる。

$\displaystyle{ -u \frac{\partial c_A}{\partial z} + r_A = \frac{\partial c_A}{\partial t} }$

反応速度

成分Aの反応速度は濃度に比例するとして次式で表す。

$\displaystyle{ -r_A=kc_A }$

$k$ は反応速度定数 [1/s]、成分Aは反応で消失するのでマイナスが付く。

初期・境界条件

初期は反応管内にAが存在しないのですべての位置で濃度ゼロ。境界条件は入口で供給されるAの濃度 $c_{A0}$ とする。

$\begin{aligned} c_A = 0 \quad &{\rm at} \quad t = 0 \\ c_A = c_{A0} \quad &{\rm at} \quad z = 0 \\ \end{aligned}$

以上で入口から成分Aを供給し始めたあとの反応器内のAの濃度変化を計算できる状態となった。

解析解

このままの形で解析解を得られるか定かでないが、長い時間が経過し定常状態に至った場合については簡単に解析解が得られる。定常状態では濃度の時間変化がないので蓄積を表す項 $\partial c_A/\partial t$ がゼロとなり、独立変数が z だけの常微分方程式になる。この条件では次の解析解が得られる。

$c_A = c_{A0}e^{-kz/u}$

数値解

まず方程式を整理して、

$\displaystyle{ \frac{\partial c_A}{\partial t} = -u\frac{\partial c_A}{\partial z} - kc_A }$

ここで位置 $z$ と時刻 $t$ を次の形で表すこととし( $i$ 、 $j$ は整数)、

$z_i = i \times \Delta z, \quad t_j = j \times \Delta t$

濃度 $c_A$ について次の表記の仕方をすることとする( $A$ は見づらくなるので省略)。

$c_{i, j} = f(z_i, t_j)$

偏微分を次の通り位置について後退差分、時間について前進差分で近似する。

$\begin{aligned} \frac{\partial c_A}{\partial z} &\fallingdotseq \frac{c_{i,j} - c_{i-1,j}}{\Delta z} \\ \frac{\partial c_A}{\partial t} &\fallingdotseq \frac{c_{i, j+1} - c_{i, j}}{\Delta t} \\ \end{aligned}$

また、項 $kc_A$ の $c_A$ は次の通り位置方向の2点の平均で近似する。

$\displaystyle{ c_A \fallingdotseq \frac{c_{i, j} + c_{i-1, j}}{2} }$

これら近似式を微分方程式に代入し整理すると下記の形となり、濃度 $c$ の値を時刻0から順次計算していくことができる。

$\displaystyle{ c_{i, j+1} = c_{i, j} - \frac{u \Delta t}{\Delta z} (c_{i, j} - c_{i-1, j}) - \frac{k \Delta t}{2}(c_{i, j} + c_{i-1, j}) }$

算出イメージは下図のようになる。横軸を $i$ (位置方向)、縦軸を $j$ (時間方向)とした格子であり、各格子点がある位置、ある時刻での濃度を表している。濃度既知の格子点は青色、濃度未知は赤色であるが、緑矢印の方向に値が求まるイメージとなる。位置について後退差分を取ったが、前進差分にすると出口(一番右の格子点)の計算ができない。