今回のターゲット：１次元熱拡散方程式
前回は1次元の熱拡散方程式を数値計算するために、2階微分に関する中心差分の式を勉強しました。
今回は実際にコードを書いて、１次元の熱拡散方程式の数値計算を実施します。
勉強の際に参考にさせていただいた書籍もご紹介させていただきます。具体的なPythonのコードも載っており、大変勉強になる１冊です。
１次元熱拡散方程式の変形
前回の復習も兼ねて、１次元熱拡散方程式と2階微分に関する中心差分の式を示します。ただし、中心差分の式は $O\left(h^2\right)$ の項を無視しています。
$$
\cfrac{\partial{u\left(x, t\right)}}{\partial{t}} = \alpha\cfrac{\partial{^2}{u\left(x, t\right)}}{\partial{x^2}} \tag{1}
$$
$$
\left(\cfrac{\partial{^{2}f}}{\partial{x^2}}\right)_{i} = \cfrac{f\left(x_{i+1}\right) -2f\left(x_{i+1}\right) + f\left(x_{i-1}\right)}{h^2} \tag{2}
$$
式(1)は２つ以上の変数を持つ微分方程式、すなわち偏微分方程式になっています。もっと言えば、時間 $t$ と位置 $x$ の偏微分方程式となります。
ここで、式(2)を使って、式(1)を変形すれば、以下の式が得られます。
$$
\cfrac{\partial{u\left(x_{i}, t\right)}}{\partial{t}} = \cfrac{u\left(x_{i+1}, t\right) -2u\left(x_{i}, t\right) + u\left(x_{i-1}, t\right)}{h^2} \tag{3}
$$
この式はよく見ると、常微分方程式の形になっています。そのため、1次元熱拡散方程式は常微分方程式として数値計算することができます。
Pythonで計算する
コードは以下のようにしました。solve_ivpで計算してみました。
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
#熱拡散方程式
def diff_eq(t, u):
#uと同じ形で1次元熱拡散方程式の左辺du/dtを全ての要素を0で初期化する
dudt = np.zeros_like(u)
#2階微分の中心差分の式を使った1次元熱拡散方程式（近似）
dudt[1:-1] = alpha * (u[:-2] - 2 * u[1:-1] + u[2:]) / dx**2
#境界の温度（物体の左端と右端）を固定したいので境界の時間変化を0に固定する
dudt[0] = 0
dudt[-1] = 0
return dudt
#各種パラメータを設定
#空間関連
L = 1.0 #物体の長さ
nx = 100 #空間分割数
x = np.linspace(0, L, nx) #座標
dx = x[1] - x[0] #空間分割
alpha = 1.12*10**(-4) #物質の熱拡散率（今回は銅 Cuを想定）
#時間関連
t_start = 0 #シミュレーションの開始時間
t_end = 5000 #シミュレーションの終了時間
t_range = (t_start, t_end) #シミュレーション時間
t_eval = np.linspace(*t_range, 21) #時間の刻み幅を指定
#初期条件を設定
u_ini = 298.15*np.ones_like(x) #今回は１次元物体の両端以外は初期温度298.15K=25 °C
# 初期の境界条件を設定（今回は１次元物体の両端が373.15K=100℃の場合を考える
u_ini[0] = 373.15 # 左端の温度
u_ini[-1] = 373.15 # 右端の温度
#solve_ivpで計算
sol = solve_ivp(diff_eq, t_range, u_ini, method="RK45",dense_output=True, t_eval=t_eval, rtol=1e-8)
#境界条件を100℃で固定する
for i in range(sol.y.shape[1]):
sol.y[0, i] = 373.15 # 左端の温度を固定
sol.y[-1, i] = 373.15 # 右端の温度を固定
#以下は可視化のためにグラフを書くコード
fig, ax = plt.subplots(dpi=240)
color_map = plt.get_cmap('jet')
N = len(sol.t)
for i in range(N):
ax.plot(x, sol.y[:, i], color=color_map(i/N), label = "t="+str(sol.t[i]))
ax.set_xlabel('x [m]')
ax.set_ylabel('Temperature [K]')
ax.legend(fontsize=8, loc='upper left', bbox_to_anchor=(1.0, 1.0))
ax.tick_params(direction='in', axis='x')
ax.tick_params(direction='in', axis='y')
plt.show()実行結果
それでは実行結果を見てみましょう。下のようなグラフが得られると思います。
初期は物体の両端点が373.15 K（=100 °C ）で、その他の部分は298.15 K（=25℃）の状態に設定しました。時間が経過するに従って物体全体の温度が上昇していきます。そして、5000 secの時には物体全体の温度が373.15 Kになっています。
この計算結果を信じれば、（仮に１次元方向しか熱拡散を考えなくても良い場合、）1.0 mの銅Cuを両端から100℃で温め始めると、全体が100℃になるまでざっくり1時間以上がかかるということでしょうか。
間違っていたらすみません。
最後に
今回は偏微分方程式を勉強してみました。
実はこの数値計算で自分が思った通りの結果を出すために2週間ほどかかり、ずっと唸っていました。特に時間がかかったのは、”どうやって物体の両端温度を固定するか”でした。
このあたりのコードをサラッと書けるようになれば良いのですが。
それでは今回はこのあたりで終わりにします。
他にも色々な記事があるのでぜひご覧ください。
コメント