【科学技術計算講座3-9】熱伝導方程式の２次元への拡張

前回はクランク=ニコルソン法を紹介し、計算精度について説明しました。

今回と次回で、２次元の熱伝導の問題を解いてみたいと思います。そして最後には温度が変化する様子をアニメーションにしてみましょう。本日は熱伝導方程式を２次元に拡張し、プログラムでの２次元データの扱いについて説明します。

問題

前回まで１次元の熱伝導の問題を解いてきました。これを２次元に拡張してみます。

問題：　次の図のような０℃に冷却された２次元の銅板がある。端の一部を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ループを回します。

さあこれだけの情報があれば、以前作った完全陰解法のプログラムを２次元に拡張できると思います。ぜひチャレンジしてみてください。

まとめ

本日は２次元熱伝導方程式へ拡張を行い、２次元データをプログラムでどう扱えばよいかを説明しました。

さて次回はいよいよ最終回です。２次元平面で熱が伝わっていく様子をアニメーションにする部分を加えてプログラムを完成させます。

←前回　クランク=ニコルソン法のプログラミング

→次回　２次元熱伝導方程式のシミュレーション

全体の目次