管内の1次元拡散方程式 - モデリングと数値計算（主に化学工学向け）

はじめに

管内の拡散を対象として次の拡散方程式を導出し、数値計算する。

$\displaystyle{ \frac{\partial c}{\partial t} = D \frac{\partial^2 c}{\partial x^2} }$

導出

f:id:schemer1341:20201219203616p:plain:w300

管内に液体が満たされていて、その中に溶解した物質が左から右へ拡散しているとする。液体中のある成分Aについて、微小区間 $\Delta x$ における微小時間 $\Delta t$ の間の物質収支は、

(拡散による微小区間へのAの流入量) - (拡散による微小区間からのAの流出量) = (微小区間内のAの蓄積量)

と考えられる。ここでAの質量流束を $j$ [kg/m2/s]、微小区間内の成分Aの濃度を $c$ [kg/m3]として物質収支を数式で表すと下記となる（流束は管の半径方向に対し一定と考える)。

$\displaystyle{ j_x \times S \times \Delta t - j_{x+\Delta x} \times S \times \Delta t = c_{t+\Delta t} \times S \times \Delta x - c_t \times S \times \Delta x }$

両辺を $S \Delta x \Delta t$ で割ると、

$\displaystyle{ \frac{j_x - j_{x+\Delta x}}{\Delta x} = \frac{c_{t+\Delta t} - c_t }{\Delta t} }$

$\Delta x \rightarrow 0, \quad \Delta t \rightarrow 0$ の極限を取ると、

$\displaystyle{ -\frac{\partial j}{\partial x} = \frac{\partial c}{\partial t} }$

ここで、フィックの法則より、

$\displaystyle{ j = - D \frac{\partial c}{\partial x} }$

これを上の式に代入すると、次の1次元拡散方程式が得られる( $D$ は拡散係数 [m2/s]で、ここでは一定値で不変とする)。

$\displaystyle{ \frac{\partial c}{\partial t}=D \frac{\partial^2 c}{\partial x^2} }$

差分近似による数値解法

まず微小距離を $\Delta x$ 、微小時間を $\Delta t$ とし、位置 $x$ と時刻 $t$ を次の形で表すこととする。

$x_i = i \times \Delta x, \quad t_j = j \times \Delta t$

ここで $i$ 、 $j$ は整数とする。濃度 $c$ は位置 $x$ と時刻 $t$ の関数であるが、 $c_{i, j}$ は次を意味するものとする。

$c_{i, j} = f(x_i, t_j) = f(i\Delta x, j\Delta t)$

このとき偏微分の差分法による近似式は次の形で表される。

$\begin{aligned} \frac{\partial c}{\partial t} &\fallingdotseq \frac{c_{i, j+1} - c_{i, j}}{\Delta t} \\ \frac{\partial^2 c}{\partial x^2} &\fallingdotseq \frac{c_{i+1, j} -2c_{i, j} + c_{i-1, j}}{(\Delta x)^2} \end{aligned}$