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

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

1次元非定常熱伝導の数値計算

はじめに

前回、1次元非定常熱伝導のモデリングを行い偏微分方程式を導いたが、今回は数値計算によりその解を得る。

方程式

加温したステンレス棒を急に冷水につける問題について、前回導いた偏微分方程式および初期・境界条件は下記のとおりである。


\displaystyle{
\frac{\partial T}{\partial t} = \alpha \Bigl( \frac{\partial^2 T}{\partial r^2} + \frac{1}{r} \frac{\partial T}{\partial r} \Bigr)
}


\begin{aligned}
T = T_0 \quad &{\rm at} \quad t = 0 \\
T = T_1 \quad &{\rm at} \quad r = R \\
\frac{\partial T}{\partial r} = 0 \quad &{\rm at} \quad r = 0
\end{aligned}

差分近似による数値解法

微小距離を \Delta r、微小時間を \Delta t とおき、位置 rと時刻 t を次の形で表すこととする。


r_i = i \times \Delta r, \quad t_j = j \times \Delta t

ここで ij は整数とする。温度 T は位置 rと時刻 t の関数であるが、次の記載の仕方をすることとする。


T_{i, j} = f(r_i, t_j)

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


\begin{aligned}
\frac{\partial T}{\partial t} &\fallingdotseq \frac{T_{i, j+1} - T_{i, j}}{\Delta t} \\
\frac{\partial T}{\partial r} &\fallingdotseq \frac{T_{i+1, j} - T_{i, j}}{\Delta r} \\
\frac{\partial^2 T}{\partial r^2} &\fallingdotseq \frac{T_{i+1, j} -2T_{i, j} + T_{i-1, j}}{(\Delta r)^2}
\end{aligned}

3つ目については、次の形で近似していると考えれば導くことができる。


\displaystyle{
\frac{\partial^2 T}{\partial r^2} \fallingdotseq \frac{(\partial T/\partial r)_{i, j} - (\partial T/\partial r)_{i-1, j} }{\Delta r}
}

これら差分近似式を微分方程式に代入して整理すると次の形となる。この式からTの値を順次算出していくことができる。


\displaystyle{
T_{i, j+1} = T_{i, j} + \alpha \Bigl( \frac{T_{i+1, j} -2T_{i, j} + T_{i-1, j}}{(\Delta r)^2} + \frac{1}{r} \frac{T_{i+1, j} - T_{i, j}}{\Delta r} \Bigr) \Delta t
}

算出イメージ

上の式を下図のイメージで捉えてみる。下図は横軸を i (位置方向)、縦軸を j (時間方向)とした格子であり、各格子点がある位置、ある時刻での温度を表している。一番上の行は j=0 における各位置の温度、すなわち初期条件を表しており、それぞれの温度は既知となる(既知の格子点は青色としている)。なお一番右の点は r=R の温度であり、境界条件により初期以外のすべての時刻の値が既知となる。

上の式は、位置 i における次の時刻の温度が、現時刻の位置 iとその前後の位置の温度から求まるという形になっている。下図で言うと緑矢印の方向に値が求まるイメージである。時刻ゼロの温度から次の時刻の温度をすべて求め、それが済んだら時刻を進めて、さらに次の時刻の温度を求める、というふうに順次温度を求めていくことができる。

f:id:schemer1341:20200715000140p:plain:w500

計算

Pythonを用いて計算した。ソースコードは下記リンク先にある。

cemodeling/conduction_1d_unsteady.py at master · nakamura-13/cemodeling · GitHub

ステンレス棒の初期温度は50℃で全体が均一とし、冷水温度は0 ℃とした。ステンレス棒の半径は1.5 cmである。位置(半径)方向の分割数は10とし(よって格子点数は11点)、時刻のステップ幅は0.1秒とした。ステンレス棒の密度、熱伝導率、比熱容量はSUS304の値を参考に決めた(具体的な値はソースコード参照)。

計算結果を下グラフに示す。横軸は半径方向位置、縦軸は温度であり、時刻0秒、5秒、10秒、20秒、40秒の温度の分布を示している。40秒後にはステンレス棒の芯までほぼ冷水温度の0℃まで低下していることがわかる。

f:id:schemer1341:20200716002841p:plain:w500

なお、下記書籍のChapter 12にある1次元非定常熱伝導の解析解グラフを参照して本問題の解析解を確認してみると、棒の芯の温度は10秒で28℃、20秒で10℃ほどになり、数値計算の結果と矛盾がないことを確認できる。

Amazon | Transport Phenomena | Bird, R. Byron, Stewart, Warren E., Lightfoot, Edwin N. | Dynamics

その他の参考書籍

Cによる数値計算法入門(第2版)新装版 | 堀之内 總一, 酒井 幸吉, 榎園 茂 |本 | 通販 | Amazon