差分法(陽解法，陰解法，Crank-Nicolson法)による1次元非定常伝熱解析

ht_fdm.py

#!/usr/bin/env python
# coding: utf-8

# 差分法(陽解法，陰解法，Crank-Nicolson法)による1次元非定常伝熱解析

import numpy as np
from numpy import matlib
from numpy import linalg
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d

# 解析条件
init_temp = 273. + 20.  # 初期温度 [K]
total_time = 600.  # 解析時間 [s]

span = 0.5  # 全長 [m]
diameter = 0.01  # 直径 [m]
heat_source = 100.  # 熱源の伝熱量 [W]

rho = 2700.  # 密度 [kg/m^3]
c = 9170.  # 比熱 [J/(kg･K)]
k = 238.  # 熱伝導率 [W/(m･K)]

# 時間分割数，X方向分割数
# M, N = (10, 10)
# M, N = (60, 50)
M, N = (600, 100)

# 図の保存
save_enable = True

# 解析定数の計算
dx = span / N  # X方向刻み幅
dt = total_time / M  # 時間刻み幅
vol_element = dx * 0.25 * np.pi * diameter ** 2  # 要素あたりの体積
qv = heat_source / vol_element  # 生成熱
Gv = np.matlib.zeros((N, 1))  # 生成項ベクトル
Gv[0] = dt * qv / (rho * c)  # x=0に生成項を追加
D = k * dt / (rho * c * dx ** 2)  # 拡散数

print("拡散数 D = {0}".format(D))

# 温度配列の初期化
T_em = np.matlib.zeros((M, N))
T_em[0, :] = init_temp  # θ(0,0)～θ(0,N)に初期温度を格納
T_im = np.matrix(T_em)
T_cnm = np.matrix(T_em)

# 陽解法の係数行列(C行列)
C = np.matlib.zeros((N, N))
C[0, 0:2] = np.matrix([[1. - D, D]])  # x=0の係数
C[N - 1, N - 2:N] = np.matrix([[D, 1. - D]])  # x=Nの係数
for i in range(0, N - 2):
    C[i + 1, i:i + 3] = np.matrix([[D, 1. - 2 * D, D]])  # x=1~N-1の係数

# 陽解法の実行(行列とベクトルの演算)
for t in range(0, M - 1):
    T_em[t + 1] = (C * T_em[t].T + Gv).T

# 陰解法の係数行列(A行列)
A = np.matlib.zeros((N, N))
A[0, 0:2] = np.matrix([[1. + D, -D]])  # x=0の係数
A[N - 1, N - 2:N] = np.matrix([[-D, 1. + D]])  # x=Nの係数
for i in range(0, N - 2):
    A[i + 1, i:i + 3] = np.matrix([[-D, 1. + 2 * D, -D]])  # x=1~N-1の係数
invA = A.I  # Aの逆行列をあらかじめ計算

# 陰解法の実行(連立1次方程式の計算)
for t in range(0, M - 1):
    T_im[t + 1] = (invA * (T_im[t].T + Gv)).T

# A行列とC行列の再定義(Crank-Nicolson法の計算用)
A[0, 0:2] = np.matrix([[1. + 0.5 * D, -0.5 * D]])  # x=0の係数
C[0, 0:2] = np.matrix([[1. - 0.5 * D, 0.5 * D]])
A[N - 1, N - 2:N] = np.matrix([[-0.5 * D, 1. + 0.5 * D]])  # x=Nの係数
C[N - 1, N - 2:N] = np.matrix([[0.5 * D, 1. - 0.5 * D]])
for i in range(0, N - 2):
    A[i + 1, i:i + 3] = np.matrix([[-0.5 * D, 1. + D, -0.5 * D]])  # x=1~N-1の係数
    C[i + 1, i:i + 3] = np.matrix([[0.5 * D, 1. - D, 0.5 * D]])
invA = A.I  # Aの逆行列をあらかじめ計算

# Crank-Nicolson法の実行
for t in range(0, M - 1):
    T_cnm[t + 1] = (invA * (C * T_cnm[t].T + Gv)).T

# 3次元プロット
fig = plt.figure(1, figsize=(17.5, 5))
fig.patch.set_facecolor("white")
# plt.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0.05, hspace=0)
plt.subplots_adjust(left=0.01, bottom=0.05, right=0.99, top=0.95, wspace=0.01, hspace=0.)

ax_em = fig.add_subplot(131, projection="3d")
ax_im = fig.add_subplot(132, projection="3d")
ax_cnm = fig.add_subplot(133, projection="3d")
x, time = np.meshgrid(np.arange(0, span, dx), np.arange(0, total_time, dt))

if N >= 20 and M >= 20:
    ax_em.plot_wireframe(x, time, T_em, cstride=int(N / 20), rstride=int(M / 20))
    ax_im.plot_wireframe(x, time, T_im, cstride=int(N / 20), rstride=int(M / 20))
    ax_cnm.plot_wireframe(x, time, T_cnm, cstride=int(N / 20), rstride=int(M / 20))
else:
    ax_em.plot_wireframe(x, time, T_em)
    ax_im.plot_wireframe(x, time, T_im)
    ax_cnm.plot_wireframe(x, time, T_cnm)

ax_em.set_xlabel(r"$x$ [m]")
ax_em.set_ylabel(r"$t$ [s]")
ax_em.set_zlabel(r"$T$ [K]")
ax_em.set_title("Explicit")

ax_im.set_xlabel(r"$x$ [m]")
ax_im.set_ylabel(r"$t$ [s]")
ax_im.set_zlabel(r"$T$ [K]")
ax_im.set_title("Implicit")

ax_cnm.set_xlabel(r"$x$ [m]")
ax_cnm.set_ylabel(r"$t$ [s]")
ax_cnm.set_zlabel(r"$T$ [K]")
ax_cnm.set_title("Crank-Nicolson")

if save_enable:
    plt.savefig("fdm_3d_x{0}t{1}_N{2}M{3}_D{4:.3f}.png".format(span, total_time, N, M, D))

# 2次元プロット
fig = plt.figure(2, figsize=(12.5, 5))
fig.patch.set_facecolor("white")
# plt.subplots_adjust(left=0.05, bottom=0.05, right=0.95, top=0.95, wspace=0.15, hspace=0.15)

plt.subplot(1, 2, 1)
plt.plot(np.arange(0, span, dx), T_em[-1].T, label="Explicit")
plt.plot(np.arange(0, span, dx), T_im[-1].T, label="Implicit")
plt.plot(np.arange(0, span, dx), T_cnm[-1].T, label="Crank-Nicolson")
plt.xlabel(r"$x$ [m]")
plt.ylabel(r"$T$ [K]")
plt.legend()
plt.title(r"Temperature $(t = M)$")

plt.subplot(1, 2, 2)
plt.plot(np.arange(0, total_time, dt), T_em[:, 0], label="Explicit")
plt.plot(np.arange(0, total_time, dt), T_im[:, 0], label="Implicit")
plt.plot(np.arange(0, total_time, dt), T_cnm[:, 0], label="Crank-Nicolson")
plt.xlabel(r"$t$ [s]")
plt.ylabel(r"$T$ [K]")
plt.legend(loc='lower right')
plt.title(r"Temperature $(x = 0)$")

if save_enable:
    plt.savefig("fdm_2d_x{0}t{1}_M{2}N{3}_D{4:.3f}.png".format(span, total_time, M, N, D))

plt.show()