Pythonでやってみた(Engineering)：1次元の非定常伝熱シミュレーション
１．概要
非定常の熱移動による温度上昇のシミュレーションをPythonを使用して計算しました。
良い文献が手に入らず合っているかわかりませんので自分用備忘録です。
２．伝熱論
２－１．Fourierの法則
1次元の熱流量はFourierの法則に従い計算式は下記の通りです。
Q：熱流量[W], k：熱伝導度(率)[W/(m・K)], \frac{dT}{dx}：温度勾配[K/m], A：伝熱面積[m^2]
$$
\begin{array}{}Q = - Ak\frac{dT}{dx}\end{array}
$$
厚みb, 熱伝導率k, 伝熱面積Aの平板壁の両面がそれぞれ一定の温度$${T_{1}}$$, $${T_{2}}$$に保たれている条件において、定常状態ではどの断面でもQが等しくなる(すなわち温度勾配$${\frac{dT}{dx}}$$が等しい)。
$$
平板壁の熱流量Q [W] = -Ak\frac{(T_{2}-T_{1})}{b} = Ak\frac{(T_{1}-T_{2})}{b}
$$
上記は定常状態ですが熱流量が[W]=[J/s]のため時間刻みごとの熱量を計算して温度上昇を計算すれば非定常状態での計算ができるはずです。次項でPythonを利用して非定常での熱移動計算を実施しました。
３．シミュレーション１：端が高温と低温
片端が高温で片端が低温であり、熱が片方向に向けて流れていく条件で計算しました。
３ー１．計算条件および要領
計算の条件および要領を記載しました。
【計算条件】
●常温25℃の鉄(棒)の両端を420℃、0℃にして、両端温度は常に維持される。
●鉄棒長さx=100mmとして、それを10分割にする（メッシュで区切る）。
●鉄の比熱c=0.444[kJ/(kg・K)]、密度ρ=7870[kg/m3]、熱伝導率k=80.3 [W/(m・K)]を使用(kの出典：はじめての化学工学)
●熱源は固定（ヒーターのように熱量の制限がない）ため、伝熱面積Aは計算しやすいように1m2と仮定
●時間ステップは100ms=0.1sで計算
(10msだと11MBくらいになりnoteで埋め込みできなかった)
【計算要領】
１．NumpyでMesh分の配列を作成して初期状態を作成する。今回は常温の配列を作成して頭とお尻だけ高温・低温に変更した。
２．配列を一つずらして計算することで温度差を配列dTで取得した。
３．各セルは熱の流入と流出の２つがある。熱収支の計算として、配列dTを用いて流入Q_highと流出Q_lowを計算してその合計を伝熱量とした。
４．温度上昇=熱量÷(比熱[kJ/(kg・K)]×重量[kg])で計算した。重量はメッシュの区切りから体積を求めて密度と合わせて計算した(A=1m2と仮定)。
５．温度上昇計算->配列へ反映をループ計算することで時間ごとの温度変化を計算した。
６．Gif化するためにFuncAnimationを使用した。
ー＞ArtistAnimationだとステップ時間を表示させるためにplt.title()を使用しても更新されず、plt.annotate()だと文字が重なったため断念
【温度勾配ΔTの計算方法の概念】
【参考：モデル低次元化／モデル縮退（モデル縮約）】
熱が伝わる100mmの物体を10分割しており計算としては粗いと感じると思います。参考までにモデル低次元化に関しては下記記事が参考になります。
３－２．コード：FuncAnimation
計算コードの注意点は下記の通りです。
●厚みdxを小さくしたり（温度勾配：$${\frac{dT}{dx}}$$増加）刻み時間dtを長くしたりすると熱流量Qが増加しすぎて値が発散する。よって時間やmeshを適切な値に調整が必要である。
●Meshを切りすぎると温度勾配：$${\frac{dT}{dx}}$$増加＋1Meshあたりの体積(重量減少)よりdxの2乗でQが大きくなるため発散しやすくなる。よってdtを小さくする必要があり計算が重くなる。
３－２－１．線形プロット：plt.plot()
コード及び結果は下記の通りです。時間ごとに熱が伝わり最終的に温度勾配は線形に近づくことが確認できました（あってるのかな？）。
[IN]
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython.display import HTML
from tqdm import tqdm
import japanize_matplotlib
np.set_printoptions(precision=3, floatmode="maxprec") #小数点3桁に変更
def calc_dx(length: int, mesh: int) ->float:
len_meter = length/mesh/1000 #総長さ[mm]をmeshで分割してmeterに変換
return len_meter
#温度条件
T_steady = 25 #定常状態(開始時)の温度
T_High, T_Low = 420, 0 #両端温度[℃]
#物性条件 ※鉄
k=80.3 #熱伝導率[W/(m・K)]
heat_spec = 0.444 #比熱[kJ/(kg・K)]
density = 7870 #[kg/m3]
#配列作成
array = np.ones(10)*T_steady #配列の長さおよび常温を設定
array[0], array[-1] = T_High, T_Low #配列の先頭と末尾に温度を設定
#環境条件
area = 1 #1m2
length_bar = 100 #条件の長さ[mm]
mesh = len(array) #分割数
dx = calc_dx(length_bar, mesh) #meshの長さ[m]
V = area * dx #メッシュの体積 [m3]
weight = V * density #メッシュ重量[kg]->比熱計算用
dt = 0.1 #時間ステップ[s] 100ms
#図形描写
fig = plt.figure(figsize=(10,10)) #土台作成
x = np.linspace(0, length_bar, mesh) #x軸のメッシュを長さに変換
def heating_iron(i, array, x, mesh, dt):
total_time = i*dt #総時間ループ回数×時間ステップ
plt.cla() #描画をクリア
array_copy = array.copy() #配列をコピー
dT = array_copy[1:] - array_copy[:-1] #高温側 - 低温側の温度差
Q = -area * k * dT/dx /1000 #熱流量Q[kW]
Q_high, Q_low = Q[:-1], Q[1:] #1次元で見た時、左と右から移動する方向
Q_total = Q_high - Q_low #熱流量の合計
dT_Total = Q_total*dt/(heat_spec * weight)
array[1:-1] = array_copy[1:-1] + dT_Total
#plot用の条件記載
im = plt.plot(x, array, c='red')
plt.grid(linestyle='dotted') #グリッド
plt.xlabel('長さ [mm]'); plt.ylabel('温度 [℃]')
plt.xlim(0, max(x)); plt.ylim(0, 450) #x,y軸の上限・下限
plt.xticks(np.arange(0, max(x)*1.1, max(x)/mesh)); plt.yticks(np.arange(0, 500, 50)) #x,y軸の刻み
plt.title(f'1次元の非定常伝熱計算_運転時間{round(total_time,2)}[s]', y= -0.1) #タイトル
ani = animation.FuncAnimation(fig, heating_iron,
fargs = [array,x, mesh,dt],
frames=500,
interval=10)
ani.save('output.gif', writer='pillow')
plt.show()
[OUT]
※Gif画像はnoteへのUpload時間のためframes=200に変更したものを使用(500:6MB、200:2.3MB)
３－２－２．heatmap：sns.heatmap()
出力を視覚的にするためにheatmap()で描写しました。前のコードとの違いは下記の通りです。。
【前コードとの違いおよび注意点】
●描写以外の計算式はすべて同じ
●Seabornのheatmapを使用するが2次元配列である必要があるためreshape()を使用して配列の次元を変更
●Colarbar（図右にある数値の色の相関図）がptl.cla()でリセットされないため重ね書きされる。よってFuncAnimationの前に一度Cbar=Trueで宣言して、関数内(ループ処理)ではFalseとすることで綺麗に描写した。
[IN]
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython.display import HTML
from tqdm import tqdm
import japanize_matplotlib
np.set_printoptions(precision=3, floatmode="maxprec") #小数点3桁に変更
def calc_dx(length: int, mesh: int) ->float:
len_meter = length/mesh/1000 #総長さ[mm]をmeshで分割してmeterに変換
return len_meter
#温度条件
T_steady = 25 #定常状態(開始時)の温度
T_High, T_Low = 420, 0 #両端温度[℃]
#物性条件 ※鉄
k=80.3 #熱伝導率[W/(m・K)]
heat_spec = 0.444 #比熱[kJ/(kg・K)]
density = 7870 #[kg/m3]
#配列作成
array = np.ones(10)*T_steady #配列の長さおよび常温を設定
array[0], array[-1] = T_High, T_Low #配列の先頭と末尾に温度を設定
#環境条件
area = 1 #1m2
length_bar = 100 #条件の長さ[mm]
mesh = len(array) #分割数
dx = calc_dx(length_bar, mesh) #meshの長さ[m]
V = area * dx #メッシュの体積 [m3]
weight = V * density #メッシュ重量[kg]->比熱計算用
dt = 0.1 #時間ステップ[s] 100ms
#図形描写
fig = plt.figure(figsize=(10,10)) #土台作成
x = np.linspace(0, length_bar, mesh) #x軸のメッシュを長さに変換
#ColarBarを一回だけ描画するために一度宣言
im = sns.heatmap(array.reshape(1,-1),
cmap='coolwarm',
annot=True,
fmt='1.0f',
cbar = True,
vmin=0, vmax=450)
def heating_iron(i, array, x, mesh, dt):
total_time = i*dt #総時間ループ回数×時間ステップ
plt.cla() #描画をクリア
array_copy = array.copy() #配列をコピー
dT = array_copy[1:] - array_copy[:-1] #高温側 - 低温側の温度差
Q = -area * k * dT/dx /1000 #熱流量Q[kW]
Q_high, Q_low = Q[:-1], Q[1:] #1次元で見た時、左と右から移動する方向
Q_total = Q_high - Q_low #熱流量の合計
dT_Total = Q_total*dt/(heat_spec * weight)
array[1:-1] = array_copy[1:-1] + dT_Total
#plot用の条件記載
im = sns.heatmap(array.reshape(1,-1),
cmap='coolwarm', #カラーマップ
annot=True, #数値を表示
fmt='1.0f', #数値フォーマット
cbar = False, #Colarbarを非表示
vmin=0, vmax=450)
plt.axis('off')
plt.title(f'1次元の非定常伝熱計算_運転時間{round(total_time,2)}[s]', y= -0.1) #タイトル
ani = animation.FuncAnimation(fig, heating_iron,
fargs = [array,x, mesh,dt],
frames=200,
interval=10)
ani.save('output_heatmap.gif', writer='pillow')
plt.show()
[OUT]
４．シミュレーション２：両端が高温
両端が高温であり中心に向かって熱が流れるパターンを計算しました。
４－１．計算条件および要領
計算要領は基本的に全章と同じです。条件に関しては下記の通りです。
【計算条件】
●常温25℃の鉄(棒)の両端を420℃にしており両端温度は常に維持される。
●鉄棒長さを前より短くx=20mmとして、それを10分割にする（メッシュで区切る）。
●鉄の比熱c=0.444[kJ/(kg・K)]、密度ρ=7870[kg/m3]、熱伝導率k=80.3 [W/(m・K)]を使用(kの出典：はじめての化学工学)
●熱源は固定（ヒーターのように熱量の制限がない）ため、伝熱面積Aは計算しやすいように1m2と仮定
●時間ステップは発散しない値で調整
４－２．コード
４－２－１．失敗例：数値の発散(FuncAnimation)
基本的な値は3章と同じで厚みを100->20mmに変更して計算しました。
[IN]
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython.display import HTML
from tqdm import tqdm
import japanize_matplotlib
def calc_dx(length: int, mesh: int) ->float:
len_meter = length/mesh/1000 #総長さ[mm]をmeshで分割してmeterに変換
return len_meter
#温度条件
T_steady = 25 #定常状態(開始時)の温度
T_High, T_Low = (420, 420) #両端温度[℃]
#物性条件 ※鉄
k=80.3 #熱伝導率[W/(m・K)]
heat_spec = 0.444 #比熱[kJ/(kg・K)]
density = 7870 #[kg/m3]
#配列作成
array = np.ones(10)*T_steady #配列の長さおよび常温を設定
array[0], array[-1] = T_High, T_Low #配列の先頭と末尾に温度を設定
#環境条件
area = 1 #1m2
thickness = 20 #条件の長さ[mm]
mesh = len(array) #分割数
dx = calc_dx(thickness, mesh) #meshの長さ[m]
V = area * dx #メッシュの体積 [m3]
weight = V * density #メッシュ重量[kg]->比熱計算用
dt = 0.1 #時間ステップ[s] 100ms
#図形描写
fig = plt.figure(figsize=(10,10)) #土台作成
x = np.linspace(0, thickness, mesh) #x軸のメッシュを長さに変換
def heating_iron(i, array, x, mesh, dt):
total_time = i*dt #総時間ループ回数×時間ステップ
plt.cla() #描画をクリア
array_copy = array.copy() #配列をコピー
dT = array_copy[1:] - array_copy[:-1] #高温側 - 低温側の温度差
Q = -area * k * dT/dx /1000 #熱流量Q[kW]
Q_high, Q_low = Q[:-1], Q[1:] #1次元で見た時、左と右から移動する方向
Q_total = Q_high - Q_low #熱流量の合計
dT_Total = Q_total*dt/(heat_spec * weight)
array[1:-1] = array_copy[1:-1] + dT_Total
#plot用の条件記載
im = plt.plot(x, array, c='red')
plt.grid(linestyle='dotted') #グリッド
plt.xlabel('長さ [mm]'); plt.ylabel('温度 [℃]')
plt.xlim(0, max(x)); plt.ylim(0, 450) #x,y軸の上限・下限
plt.xticks(np.arange(0, max(x)*1.1, max(x)/mesh)); plt.yticks(np.arange(0, 500, 50)) #x,y軸の刻み
plt.title(f'1次元の非定常伝熱計算_運転時間{round(total_time,2)}[s]', y= -0.1) #タイトル
ani = animation.FuncAnimation(fig, heating_iron,
fargs = [array,x, mesh,dt],
frames=30,
interval=500)
ani.save('output2.gif', writer='pillow')
plt.show()
[OUT]
結果として瞬間的(各ステップ)で隣のメッシュに過剰に熱移動が生じて温度の逆転が生じたため数値が発散しました。
対策は過剰な熱移動[kJ]がないように温度勾配$${\frac{dT}{dx}}$$を小さく（Meshを荒くする、厚みを大きくする）するか、時間の刻みをdtを小さくする必要があります。厚みは固有値でMeshもすでに荒いため、今回はdtを小さくしますが同じ時間計算すると計算回数が大きくなりメモリを圧迫します。
そこでArtistAnimationを利用して、表示するデータ保存を間引いてやることでメモリ負担を小さくしてGifを作成しました。
４－２－２．成功例（ArtistAnimation）
ArtistAnimationを使用したコードのポイントは下記の通りです。なお本当はFuncAnimationのように画像の中に動くTotal時間を表示したかったのですが、私では無理だったので諦めました。
【コードのポイント】
●時間の刻みdtを0.1->0.001に変更して1回あたりの温度変化を小さくすることで発散を防止した。
●dtを小さくする代わりに処理回数(FuncAnimationでのframes)を増やした。(同じ運転時間をみるならdtが1/100のため100倍にする)
●処理回数が多く全データをGif用リストに入れるとメモリを圧迫するため、特定回数ごとに処理結果を保存する
[IN]
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython.display import HTML
from tqdm import tqdm
import japanize_matplotlib
def calc_dx(length: int, mesh: int) ->float:
len_meter = length/mesh/1000 #総長さ[mm]をmeshで分割してmeterに変換
return len_meter
#温度条件
T_steady = 25 #定常状態(開始時)の温度
T_High, T_Low = (420, 420) #両端温度[℃]
#物性条件 ※鉄
k=80.3 #熱伝導率[W/(m・K)]
heat_spec = 0.444 #比熱[kJ/(kg・K)]
density = 7870 #[kg/m3]
#配列作成
mesh =10 #分割数
array = np.ones(mesh)*T_steady #配列の長さおよび常温を設定
array[0], array[-1] = T_High, T_Low #配列の先頭と末尾に温度を設定
#環境条件
area = 1 #1m2
thickness = 20 #条件の長さ[mm]
dx = calc_dx(thickness, mesh) #meshの長さ[m]
V = area * dx #メッシュの体積 [m3]
weight = V * density #メッシュ重量[kg]->比熱計算用
dt = 0.001 #時間ステップ[s] 1ms
#図形描写
fig = plt.figure(figsize=(10,10)) #土台作成
x = np.linspace(0, thickness, mesh) #x軸のメッシュを長さに変換
plt.grid(linestyle='dotted') #グリッド
plt.xlabel('厚み [mm]'); plt.ylabel('温度 [℃]')
plt.xlim(0, max(x)); plt.ylim(0, 450) #x,y軸の上限・下限
plt.xticks(np.arange(0, max(x)*1.1, max(x)/mesh)); plt.yticks(np.arange(0, 500, 50)) #x,y軸の刻み
anims = [] #データ格納用
for i in tqdm(range(10000)):
total_time = i*dt #総時間ループ回数×時間ステップ
array_copy = array.copy() #配列をコピー
dT = array_copy[1:] - array_copy[:-1] #高温側 - 低温側の温度差
Q = -area * k * dT/dx /1000 #熱流量Q[kW]
Q_high, Q_low = Q[:-1], Q[1:] #1次元で見た時、左と右から移動する方向
Q_total = Q_high - Q_low #熱流量の合計
dT_Total = Q_total/(heat_spec * weight) *dt #各セルの温度変化
array[1:-1] = array_copy[1:-1] + dT_Total #温度の更新
#plot用の条件記載
if i%100 == 0:
im = plt.plot(x, array, c='red')
plt.title(f'1次元の非定常伝熱計算_運転時間{total_time}[s]', y= -0.15) #タイトル
anims.append(im) #データ格納
ani = animation.ArtistAnimation(fig, anims, interval=200)
ani.save('output2.gif', writer='pillow')
plt.show()
[OUT]
各Frameで100ms毎の動きを表示しているが実際は1ms毎に計算している
参考記事
あとがき
大学の時に非定常計算をちゃんと理解できてなかったからかなり難しい・・・・・・
あとこういう時に変数名のつけ方にセンスがでるな・・・・