非定常の熱平衡方程式の基本

比熱と熱伝導率を定数としたときの1次元の非定常熱平衡方程式は以下のように表される。

\frac{\partial u}{\partial t} = \frac{k}{\lambda c} \frac{\partial^2 u}{\partial x^2}

ここで、 $k$ : 熱伝導率 (W/m K)、 [ $\lambda$ ]: 密度 (g/m $^3$ )、 $c$ : 比熱 (J/g K)、 $u$ : 温度 (K)、 $t$ : 時刻 (s)、 $x$ : 距離 (m)する。

陰解法やクランクニコルソン法などあるが、まずは簡単な陽解法と呼ばれる解き方で数値計算してみる。

まず左辺を差分法を用いて表すと、

\frac{\partial u}{\partial t} = \frac{u_{i, j+1} - u_{i, j}}{\Delta t}

同様に右辺については

\frac{k}{\lambda c} \frac{\partial^2 u}{\partial x^2} = \frac{k}{\lambda c} \frac{ u_{i-1, j} - 2 u_{i, j} + u_{i+1, j}}{(\Delta x)^2}

と表されます。ここで $i$ と $j$ は距離 $x$ と時間 $t$ を離散化したときのそれぞれの格子点を表す。

右辺と左辺をまとめると以下のようになる。

u_{i, j+1} - u_{i, j} = r (u_{i-1, j} - 2 u_{i, j} + u_{i+1, j})

ここで $r = \frac{\Delta t k}{\lambda c (\Delta x)\^2}$ としています。

この式を行列式で書くと以下のようになる。

\mathbf{u^{ j+1}} = A \mathbf{u^{j}} + C

このように陽解法では既知の温度値から次の時間ステップの温度値を求める方法。これをJuliaで書いてみる。

Juliaで数値計算する

まずは解析モデルについて以下のようなモデルを想定する。

材料は銅
室温（290 K）中
片端部に600 Kを設定（片端部にハンダごてをあてた）
もう片方の端部は室温温度一定とする
$\Delta x = 0.01$ , $\Delta t = 1$ E-6

以下、クソコード実装。


using LinearAlgebra # 今回は行列計算で使用

# 定数
dx = 0.01; # delta x
dt = 1e-6; # delta t
xnum = 100; # xの要素数（格子点の数）
tnum = 1e6; # tの要素数
Ti = 290; # (K) 初期導体温度（初期条件）
Tleft = 600; # (K) 左側温度（境界条件）
Tright = 290; # (K) 右側温度（境界条件）
c = 0.379 # (J/ g K) 銅の比熱
k = 398 # (W/m K) 銅の熱伝導率
density = 8.96 * 1e6 # (g/m^3) 銅の密度

r = dt * k / c / density / dx^2 # 係数の計算

matA = diagm(0 => (1-2*r) * ones(xnum), -1 => r * ones(xnum-1), 1 => r * ones(xnum-1)) # 係数行列Aの作成
matc = [r *Tl; zeros(xnum-2); r *Tr];
for t in 1:tnum
    newx = matA * x + matC;
    x = newx;
end
p = plot(newx)
savefig(p, "data.png")

このようにして1秒後の銅の1次元温度分布を求めることができる。ただし陽解法の場合は収束条件があり

r \leq \frac{1}{2}

の場合の解が収束する。