#はじめに
私は大学で半導体物性を研究しています。
半導体材料にレーザーを照射し、光吸収による熱膨張を圧電素子(膨張や収縮を電圧に変換する素子)で計測する実験などを行っており、実験結果を計算で再現することで熱伝導率などの物性値を算出しようとしています。

その計算のために熱拡散方程式を数値計算で解きたいのです。
実験結果を再現するには1次元では不十分だったので、2次元の熱拡散方程式を解きました。

この記事では、2次元の熱拡散方程式を
ADI法(Alternating direction implicit method,交互方向陰解法)で解きます。
間違いやより良いプログラムの書き方があればどうか教えて下さい。

解きたい方程式は以下の通りです。

\frac{\partial T}{\partial t}=\lambda \left(\frac{\partial^2 T}{\partial x^2}+\frac{\partial^2 T}{\partial y^2}\right)

ここで、$T$は温度、$t$は時間、$\lambda$は熱拡散率、$x,y$は座標です。

熱拡散方程式については、熱伝導方程式の導出などで勉強しました。

また、1次元の熱拡散方程式の数値計算は流体シミュレーション : 拡散方程式を実装するなどで勉強しました。

1次元では棒状、2次元では面状、3次元では立体の温度を計算できます。

1次元の熱拡散方程式の数値計算でよく使われるクランクニコルソン法を2次元に適用すると計算量が膨大になります。一方、ADI法を用いると計算量を抑えつつ2次元の計算が可能になります。

#ADI法のイメージ
ADI法ではある時刻$k$での温度がすべて求められているとして次の$k+1$時刻を求める時、図のように$y$座標$j=1,2,...m$のそれぞれで$x$方向に$i=1,2,...n$を計算することで面全体の温度を求めます。
その次の$k+2$時刻を求める時は$x$座標$i=1,2,...n$のそれぞれで$y$方向に$j=1,2...m$を計算することで全体の温度を求めます。
$k+3$時刻は$x$方向に、$k+4$時刻は$y$方向に...と繰り返します。
$x$方向と$y$方向を交互に計算するので交互方向陰解法と言われているのだと思います。

(普通は左下に原点があり、横方向にx軸、上方向にy軸を書きますが、Spyderの変数エクスプローラで解のndarrayを開いたときに左上に0があって下方向にx、右方向にyの解があるので私の中でそのイメージになってしまっています。わかりにくくてすみません。)

#ADI法の計算

第$k$時刻での$T_{i,j}^{k}$がすべて計算されたとして、次の第$k+1$時刻での$T_{i,j}^{k+1}$を計算するために「はじめに」で示した熱拡散方程式を次のように差分化する。
$$
\frac{T_{i,j}^{k+1} - T_{i,j}^k}{\Delta t} = \lambda \left(\frac{T_{i-1,j}^{k+1} - 2 T_{i,j}^{k+1} + T_{i+1,j}^{k+1} }{\Delta x^2} + \frac{T_{i,j-1}^k - 2 T_{i,j}^k + T_{i,j+1}^k }{\Delta y^2}\right)
$$

これを左辺に$k+1$右辺に$k$となるように式変形すると次のようになる。
$$
T_{i-1,j}^{k+1}-\left(2+\frac{\Delta x^2}{\lambda \Delta t}\right)T_{i,j}^{k+1}+T_{i+1,j}^{k+1} = \left[2 \left( \frac{\Delta x^2}{\Delta y^2}\right)-\frac{\Delta x^2}{\lambda \Delta t}\right]T_{i,j}^k-\frac{\Delta x^2}{\Delta y^2} \left( T_{i,j-1}^k + T_{i,j+1}^k \right)
$$
これを次のように置く
$$
T_{i-1,j}^{k+1}+dT_{i,j}^{k+1}+T_{i+1,j}^{k+1} = BT_{i,j}^k-\frac{\Delta x^2}{\Delta y^2} \left( T_{i,j-1}^k + T_{i,j+1}^k \right)
$$

ここで、
$$
d = -\left(2+\frac{\Delta x^2}{\lambda \Delta t}\right), B=\left[2 \left( \frac{\Delta x^2}{\Delta y^2}\right)-\frac{\Delta x^2}{\lambda \Delta t}\right]
$$

$i=1,2,...m$とおき、これを行列で表すと、


\left(
\begin{array}{cccc}
      d & 1 &  0 & \ldots  & \ldots & 0 \\
      1 & d & 1 & 0 & \ldots & 0 \\
      0  &1 & d & 1 & 0 & \ldots  \\
      \vdots & \ddots & \ddots & \ddots & \ddots & \ddots\\
      0  & 0 & \ddots & 1 & d & 1 \\
      0  & 0 & \ldots & 0 & 1 & d
    \end{array}
  \right)
  \left(
    \begin{array}{c}
      T_{1,j}^{k+1}  \\
      T_{2,j}^{k+1}  \\
      T_{3,j}^{k+1}    \\
      \vdots \\
      T_{m-1,j}^{k+1} \\
      T_{m,j}^{k+1} 
\end{array}
\right)


= \left( \begin{array}{c}
    B T_{1,j}^{k} -\frac{\Delta x^2}{\Delta y^2} \left( T_{1,j-1}^k + T_{1,j+1}^k \right)- T_{0,j}^k \\
      B T_{2,j}^{k} -\frac{\Delta x^2}{\Delta y^2} \left( T_{2,j-1}^k + T_{2,j+1}^k \right)  \\
      B T_{3,j}^{k} -\frac{\Delta x^2}{\Delta y^2} \left( T_{3,j-1}^k + T_{3,j+1}^k \right)  \\
      \vdots \\
      B T_{m-1,j}^{k} -\frac{\Delta x^2}{\Delta y^2} \left( T_{m-1,j-1}^k + T_{m-1,j+1}^k \right)  \\
      B T_{m,j}^{k} -\frac{\Delta x^2}{\Delta y^2} \left( T_{m,j-1}^k + T_{m,j+1}^k \right)-T_{m+1,j}^k
    \end{array} \right)

$m$元連立一次方程式が得られる。

この式変形ではディリクレ境界条件(境界の温度が一定であるという条件)使っている。
$T_{0,j}^k$や$T_{m+1,j}^k$は境界の温度であり既知数であるから未知数は$m$個。

この連立方程式を$j=1$から$j=n$のそれぞれの場合に対して解くと第$k+1$時刻での温度が求まる.

次に、第$k+1$時刻での$T_{i,j}^{k+1}$がすべて計算されたとして、次の第$k+2$時刻での$T_{i,j}^{k+2}$を計算するために次のように差分化する。
$$
\frac{T_{i,j}^{k+2} - T_{i,j}^{k+1}}{\Delta t} = \lambda \left(\frac{T_{i-1,j}^{k+1} - 2 T_{i,j}^{k+1} + T_{i+1,j}^{k+1} }{\Delta x^2} + \frac{T_{i,j-1}^{k+2} - 2 T_{i,j}^{k+2} + T_{i,j+1}^{k+2} }{\Delta y^2}\right)
$$

同様にして
$$
T_{i,j-1}^{k+2}-\left(2+\frac{\Delta y^2}{\lambda \Delta t}\right)T_{i,j}^{k+2}+T_{i,j+1}^{k+2} = \left[2 \left( \frac{\Delta y^2}{\Delta x^2}\right)-\frac{\Delta y^2}{\lambda \Delta t}\right]T_{i,j}^{k+1}-\frac{\Delta y^2}{\Delta x^2} \left( T_{i-1,j}^{k+1} + T_{i+1,j}^{k+1} \right)
$$
これを次のように置く
$$
T_{i,j-1}^{k+2}-eT_{i,j}^{k+2}+T_{i,j+1}^{k+2} = CT_{i,j}^{k+1}-\frac{\Delta y^2}{\Delta x^2} \left( T_{i-1,j}^{k+1} + T_{i+1,j}^{k+1} \right)
$$

ここで、
$$
e=-\left(2+\frac{\Delta y^2}{\lambda \Delta t}\right), C=\left[2 \left( \frac{\Delta y^2}{\Delta x^2}\right)-\frac{\Delta y^2}{\lambda \Delta t}\right]
$$

第$k+1$時刻から次の第$k+2$時刻の温度を計算する場合も同様にして行列を書ける。(省略します)

$j=1,2,...n$とおくと、$n$個の未知数に対する$n$元連立1次方程式が得られる。
その連立方程式を$i=1$から$i=m$までの各場合にそれぞれ解くと第$k+2$時刻の温度が求まる。

ADI法では1回に解くべき連立方程式の元数は$m$元か$n$元で良いことになる。

クランクニコルソン法で同じように計算しようとすると$m\times n$元連立1次方程式を解かなければならない。

#プログラム

計算の条件は

縦10、横10の大きさ、熱拡散率が0.1の物体
$dx=dy=0.1$, $dt=0.1$として計算
時刻0の時物体全体の温度が0度
物体の4辺の温度は100度で固定されている


# -*- coding: utf-8 -*-
"""
交互方向陰解法による2次元の熱拡散方程式の計算
"""

import numpy as np
import numpy.matlib
from decimal import Decimal, ROUND_HALF_UP
import scipy.sparse.linalg as spla
from scipy.sparse import csr_matrix
import matplotlib.pyplot as plt
import seaborn as sns
import time


dtnum = 5001 #時間をいくつに分割して計算するか 一の位を1とした方が良い 時刻tのfor文の終わりの数字に関係あり
dxnum = 101 #xとyをいくつに分割して計算するか
dynum = 101

thick = 10 #x方向の大きさ
width = 10 #y方向の大きさ
time_calc = 500 #計算する時間
beta = 0.1 #上の式中ではラムダ　温度拡散率 lambda = 熱伝導率/(密度*比熱)

"""計算準備"""
#空の解を用意
solution_k1 = np.zeros([dxnum,dynum])
solution_k2 = np.zeros([dxnum,dynum]) #solution_k2の初期値が初期条件になる


#境界条件
irr_boundary = 100 #表面の境界条件 (0,t)における温度 
rear_boundary = 100 #裏面と反対側の境界条件　(x,t)における温度 0を基準となる温度とした
side_0 = 100 #y=0の温度
side_y = 100 #y=0とは反対側の境界の温度

dt = time_calc / (dtnum - 1)
dx = thick / (dxnum - 1)
dy = width / (dynum - 1)

d = -(2+(dx**2)/(beta*dt))
B = (2*((dx/dy)**2)-((dx**2)/(beta*dt)))

e = -(2+(dy**2)/(beta*dt))
C = (2*((dy/dx)**2)-((dy**2)/(beta*dt)))


"""Ax=bのAを用意"""
a_matrix_1 = np.identity(dxnum) * d \
            + np.eye(dxnum, k=1) \
            + np.eye(dxnum, k=-1)
            
a_matrix_2 = np.identity(dynum) * e \
            + np.eye(dynum, k=1) \
            + np.eye(dynum, k=-1)      
    
#疎行列を格納　csr方式
a_matrix_csr1 = csr_matrix(a_matrix_1)
a_matrix_csr2 = csr_matrix(a_matrix_2)


#ADI法ではk+1時刻とk+2時刻を計算するので、for文の回数は半分になる
number = Decimal(str(dtnum/2)).quantize(Decimal('0'), rounding=ROUND_HALF_UP)
number = int(number)

solution = np.zeros([dxnum,dynum,number]) #解を代入する行列を作成

#temp_temperature_arrayは、ずれた列を作成して、足すためのもの
temp_temperature_array1 = np.zeros([dxnum,dynum+2])
temp_temperature_array1[:,0] = side_0
temp_temperature_array1[:,-1] = side_y

temp_temperature_array2 = np.zeros([dxnum+2,dynum])
temp_temperature_array2[0,:] = irr_boundary
temp_temperature_array2[-1,:] = rear_boundary

#ADI法の計算
for k in range(number): #時刻tについて
    for j in range(dynum):#xを計算し、yの数繰り返すことで、時刻k+1のTijを計算       
        temp_temperature_array1[:,1:dynum+1] = solution_k2
        b_array_1 = B * temp_temperature_array1[:,j+1] \
            - ((dx/dy)**2) * (temp_temperature_array1[:,j] + temp_temperature_array1[:,j+2])
        
        # bの最初と最後に境界条件
        b_array_1[0] -= irr_boundary
        b_array_1[-1] -= rear_boundary
        
        #解を求める
        temperature_solve1 = spla.dsolve.spsolve(a_matrix_csr1, b_array_1)#xについての解
        solution_k1[:,j] = temperature_solve1
        
    for i in range(dxnum):#yを計算し、xの数繰り返すことで、時刻k+2のTijを計算
        temp_temperature_array2[1:dxnum+1,:] = solution_k1
        b_array_2 = C * temp_temperature_array2[i+1,:] \
            - ((dy/dx)**2) * (temp_temperature_array2[i,:] + temp_temperature_array2[i+2,:])

        # bの最初と最後に境界条件
        b_array_2[0] -= side_0
        b_array_2[-1] -= side_y
        
        #解を求める
        temperature_solve2 = spla.dsolve.spsolve(a_matrix_csr2, b_array_2)#yについての解
        solution_k2[i,:] = temperature_solve2
    solution[:,:,k] = solution_k2


ax = sns.heatmap(solution[:,:,10], linewidth=0, vmin=0, vmax=100)
plt.show()

ax = sns.heatmap(solution[:,:,100], linewidth=0, vmin=0, vmax=100)
plt.show()

ax = sns.heatmap(solution[:,:,500], linewidth=0, vmin=0, vmax=100)
plt.show()

ax = sns.heatmap(solution[:,:,2000], linewidth=0, vmin=0, vmax=100)
plt.show()

実行していただくと私のパソコンだと50秒ぐらいかかって計算結果が出てきます。
計算結果は以下の通りです。

初期条件の0度から時間の経過とともに境界付近から次第に熱が伝わり100度に近づいていく様子が計算できています。

#終わりに
適切な刻み幅ってあるんでしょうか。
物体の大きさがマイクロやナノメートルオーダーで計算する時刻がミリ秒オーダーの時はどの程度の刻み幅が必要でしょうか。
陰解法でも$dx$と$dy$は$dt$より小さい方が良いのでしょうか。
教えていただけますと幸いです。

#参考文献

山崎郭滋, 偏微分方程式の数値解法入門, 森北出版株式会社 (1993).
http://nalab.mind.meiji.ac.jp/~mk/labo/text/heat-fdm-1.pdf
https://qiita.com/KQTS/items/97daa509991bb9777a4a