1次元熱平衡式をJuliaで数値計算してみる（陽解法）

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

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

\displaystyle \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)とします。

まずは第一のステップとして陽解法と呼ばれる解き方で数値計算してみましょう。

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

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

同様に右辺については

\displaystyle \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を離散化したときのそれぞれの格子点を示します。

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

\displaystyle 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} としています。

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

f:id:morita0079:20210714181933p:plain
equation 1

まとめると

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

このように陽解法では既知の温度値から次の時間ステップの温度値を求める方法です。これをプログラミング言語Juliaで書いてみましょう。

Juliaで数値計算する

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

  • 材料は銅
  • 室温（290 K）中
  • 片端部に600 Kを設定（片端部にハンダごてをあてた）
  • もう片方の端部は室温温度一定とする
  • \Delta x = 0.01, \Delta t = 1E-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")

以上のコードを実行すると以下のような結果が得られます。 f:id:morita0079:20210715022846p:plain

単に温度配列をプロットしただけなので横軸は xの格子点になってます。縦軸は温度です。 このようにして1秒後の銅の1次元温度分布を求めることができました。 ただし陽解法の場合は収束条件があり

\displaystyle r \leq \frac{1}{2}

の場合の解が収束します。