差分法で熱伝導方程式を解く

差分法では、偏導関数を次のように近似する。

　（前進差分）もしくは、
　（後退差分）
2階の偏微分は1階の偏導関数の差分であるから、前進差分と後退差分を用いれば、次式が得られる。

　（中心差分）

---

熱伝導方程式

差分法を用いて熱伝導方程式の数値解を求めてみる。一次元の熱伝導方程式は次式で表される。


Tは温度 [K]、αは温度伝導率（温度拡散率）であり、次式より求まる。

ここで、λは熱伝導率 [W/(m・K)]、ρは密度 [kg/m³]、c_pは定圧比熱 [J/(kg・K)]である。

---

陽解法

前進差分を用いて数値解を求める方法は陽解法と呼ばれる。T_x,tとして、陽解法による差分近似式は次式となる。


T_i,j+1について解けば、


が得られる。ここで、γ=Δt/Δx²である。この式より、T_i-1,j, T_i,j, T_i+1,jの3点から、次の時刻のT_i,j+1を計算することができる。

Pythonでプログラムを書くとこうなる。

#coding: utf-8

L = 0.1  #長さ[m]
M = 10   #分割数
dx = L/M #空間刻み[m]

dt = 1.0 #時間刻み[s]
N = 100  #時間ステップ数

#温度伝導率(鉄)
alpha = 80.2/(7874.0*440.0)

gamma = dt/dx**2
a = alpha*gamma

#初期条件(温度273.15K)
T = [273.15 for i in range(M+1)]

#境界条件
T[0] = 300.0  #温度固定
T[M] = T[M-1] #断熱

f = open('output','w')

for j in range(1,N):
    for i in range(1,M):
        preT = T

        #陽解法による差分式
        T[i] = a*preT[i+1]+(1-2*a)*preT[i]+a*preT[i-1]
        T[0] = 300.0
        T[M] = T[M-1]

    #計算結果をファイルへ出力
    f.write('# t=%ss\n'%(j*dt))
    for i in range(M+1):
        f.write('%s, %s\n'%(i*dx,T[i]))
    f.write('\n\n')

f.close()

グラフにプロットするとこんな感じ。時間とともに、左側から温度が上昇していく様子がわかる。

---

陰解法

後退差分を用いて数値解を求める方法は陰解法と呼ばれる。差分近似式は次式となる。


整理すると、


陽解法とは異なり、T_i,jから、次の時刻でのT_i-1,j+1, T_i,j+1, T_i+1,j+1の3点を求めなければならない。前式は行列で書くと次のようになる。


ただし、行列のT_1,j+1とT_M,j+1には境界条件を設定する。従って、T_{(1,2,...,M),j+1}は次式より求まる。


Python + numarray でプログラムを書くとこうなる。

#coding: utf-8

import numarray
import numarray.linear_algebra as la

L = 0.1  #長さ[m]
M = 10   #分割数
dx = L/M #空間刻み[m]

dt = 1.0 #時間刻み[s]
N = 100  #時間ステップ数

#温度伝導率(鉄)
alpha = 80.2/(7874.0*440.0)

gamma = dt/dx**2
a = alpha*gamma

#初期条件(温度273.15K)
T = numarray.array([[273.15,] for i in range(M+1)])
#境界条件
T[0][0] = 300.0  #温度固定
T[M][0] = T[M-1][0] #断熱

#行列の作成
A = numarray.array([[0. for i in range(M+1)] for j in range(M+1)])
for i in range(1,M):
    A[i][i-1] = -a
    A[i][i]   = 1+2*a
    A[i][i+1] = -a
#境界条件
A[0][0] = 1.   #一定温度
A[M][M]   = 1. #断熱

#Aの逆行列
A_inv = la.inverse(A)

f = open('output','w')

for j in range(1,N):
    preT = T

    T = numarray.dot(A_inv, preT)
    T[0][0] = 300.0
    T[M][0] = T[M-1][0]

    #計算結果をファイルへ出力
    f.write('# t=%ss\n'%(j*dt))
    for i in range(M+1):
        f.write('%s, %s\n'%(i*dx,T[i][0]))
    f.write('\n\n')

f.close()

グラフにプロットするとこんな感じ。

---

クランク・ニコルソン法