前回はクランク=ニコルソン法を紹介し、計算精度について説明しました。
熱伝導方程式をクランク=ニコルソン法で解くPythonプログラムを作ります。スキームの精度についても説明します。クランク=ニコルソン法は２次精度で無条件安定なスキームです。科学技術計算講座３「熱伝導方程式のシミュレーション」の第８回目です。
cattech-lab.com
2020-05-10 14:30
今回と次回で、２次元の熱伝導の問題を解いてみたいと思います。そして最後には温度が変化する様子をアニメーションにしてみましょう。本日は熱伝導方程式を２次元に拡張し、プログラムでの２次元データの扱いについて説明します。
問題
前回まで１次元の熱伝導の問題を解いてきました。これを２次元に拡張してみます。
問題： 次の図のような０℃に冷却された２次元の銅板がある。端の一部を100℃に固定したとき温度の変化するようすを求めたい。物性値などはこれまでの問題と同じとする。また、100℃に固定した以外の面は全て断熱である。
今度は２次元の板の熱伝導の問題です。
２次元熱伝導方程式
まずは数式を考えましょう。きちんと導出からやろうとすると、前と同じように微小体積の熱の出入りを考えてやればよいです。方程式としては、１次元の熱伝導方程式の空間項に $y$ 方向の空間微分が追加された式となります。
$$\frac{\partial T}{\partial t} = \alpha (\frac{\partial^2 T}{\partial x^2} +\frac{\partial^2 T}{\partial y^2}) \tag{1}$$
これが２次元熱伝導方程式です。
離散化
次に(1)式を数値的に解くために離散化を考えます。今度は空間に $y$ 方向が加わったので、この方向にも離散化してみます。
$y$ 方向を $\Delta y$ の幅で分割して離散化すると図のようになります。
ここで、$x$ 方向の点の添字を $i$、$y$ 方向の点の添字を $j$ として、位置 $(x_i, y_j)$ の温度を $T_{i,j}$ と書くことにします。
時間方向にはこれまでと同じように上付きの添字 $n$ をつけて、$T_{i,j}^n$ とします。
完全陰解法
(1)式を解く数値計算法ですが、ここでは完全陰解法を使うことにしましょう。(1)式を完全陰解法で表記すると、
$$\frac{T_{i,j}^{n+1}-T_{i,j}^n}{\Delta t} = \alpha (\frac{T_{i+1,j}^{n+1} - 2 T_{i,j}^{n+1} + T_{i-1,j}^{n+1}}{\Delta x^2}+\frac{T_{i,j+1}^{n+1} - 2 T_{i,j}^{n+1} + T_{i,j-1}^{n+1}}{\Delta y^2}) \tag{2}$$
となります。空間の差分をとるときは、微分する方向の点に対して差分をとっていけばよいです。
(2)式を $T_{i,j}^{n+1} $ について整理すると最終的に次式のようになります。
$$T_{i,j}^{n+1} = \frac{1}{1+2c_x+2c_y} ( T_{i,j}^n + c_x T_{i+1,j}^{n+1} + c_x T_{i-1,j}^{n+1}+ c_y T_{i,j+1}^{n+1} + c_y T_{i,j-1}^{n+1}) \tag{3}$$
ここで、$c_x = \alpha \Delta t/ \Delta x^2$、$c_y = \alpha \Delta t/ \Delta y^2$。
あとは(3)式をガウス=ザイデル法で解いていけば解が求まります。
どうでしょうか、２次元への拡張はそれほど難しくはないと思います。これまでの式に機械的に当てはめていけばよいのです。
プログラミング
ここまで出来たらあとはプログラミングしていくだけですが、今回のように２次元のデータはどのように扱えばよいでしょうか。
これまでは一次元のデータを配列temp[i] を使って表してきました。２次元のデータも配列で表します。この配列は２次元配列と呼ばれます。
２次元配列の要素はtemp[i][j] のように表します。$T_{i,j}$ と全く同じように考えればよいです。２次元配列は次のように箱が並んでいるとイメージしてください。
Pythonでは、要素数がnx個×ny個で、全ての値が0.0 の２次元配列を作るためには、
temp = np.zeros((nx, ny))
値が100.0 の2次元配列を作るには、
temp = np.full((nx, ny), 100.0)
と書きます。
また、全ての要素に対して処理をしようとすると、
for i in range(nx): for j in range(ny): temp_old[i][j] = temp[i][j]
のように、i とj でforループを回します。
さあこれだけの情報があれば、以前作った完全陰解法のプログラムを２次元に拡張できると思います。ぜひチャレンジしてみてください。
まとめ
本日は２次元熱伝導方程式へ拡張を行い、２次元データをプログラムでどう扱えばよいかを説明しました。
さて次回はいよいよ最終回です。２次元平面で熱が伝わっていく様子をアニメーションにする部分を加えてプログラムを完成させます。
←前回 クランク=ニコルソン法のプログラミング
→次回 ２次元熱伝導方程式のシミュレーション
全体の目次
