SUS316の伝熱を熱回路網法(Thermal-Network)でやってみました

こんちには。

最近の興味は熱の挙動です。実は私は大学時代、機械系で材料系でって感じでFEMとか差分法やってました。

「熱」なんてさらさら興味がなく、エンタルピーやらエントロピーやら熱効率やら対流やら、、、テストでは一夜漬けで対応してぎりぎり単位取得みたいな典型的な大学生でした。

ところが、数年前、設計業務で初めて「熱」を取り扱うことになったのです。もう何すれば良いのか全く分からず、とりあえず大学時代の伝熱工学の教科書を引っ張ってきました。ここまで読んで頂いて分かったと思いますが、自分は見た目と裏腹に結構真面目です。本質的には物理も数学も好きですから。

前置きはここまでにして、熱回路網法(Thermal-Network)は熱の専門家の方々には常識と思います。ちょっとこれを用いて、SUS316の板の伝熱(定常2次元解析)をやってみました。下図のような感じで伝熱モデルを考えます。発熱体の出力は1Wとしますので、5個の要素で分割してみると1要素あたり0.2Wの発熱量です。モデル寸法はx=0.3m,y=0.1m,z=0.1mです。SUS316の熱伝導率はざっくり16.7W/m・K, 右端の大気への放熱による自然対流熱伝達率は10W/(m^2・K)で計算してみます。

f:id:Aobayama-Engineering:20181111191621p:plain

 excel上で入力シートは下図のような感じで作ってみます。熱コンダクタンスを各節点に定義してみます。

f:id:Aobayama-Engineering:20181111220226p:plain

 これらをマクロで読み込んで計算した結果が下記です。セル内の数値は温度[℃]です。

f:id:Aobayama-Engineering:20181111221106p:plain

1Wの発熱だとぬくい感じですね。定性的にも合ってそうですね。

検算でプログラムと数値解が合っていそうか確かめてみます。

伝熱の基礎式、コンダクタンスC[W/K]=伝熱量Q[W]/温度差⊿T[K]で手計算してみます。

私の手計算結果によれば、右端は35℃で左端は36.8℃でした。簡単なモデルなので、そこそこに精度が出ていますね。

とりあえず一安心です。

本計算で用いたプログラムを下記に書いておきますね。マクロにコピペして、入力シートを作れば使えるはずです。動作に不具合あっても補償しません。笑

本日はここまで。

 

↓プログラム↓

Sub thermal_net_steady1()

'型宣言
Dim NODE As Integer '節点数
Dim ELM As Integer '要素数
Dim HNODE As Integer '発熱する節点数
Dim TNODE As Integer '固定温度の節点数
Dim N1 As Integer '節点1 No.
Dim N2 As Integer '節点2 No.
Dim CELM As Single '要素の熱コンダクタンス
Dim A() As Single '熱伝導マトリックス
Dim B() As Single '荷重ベクトル
Dim T() As Single '温度ベクトル

'入力
Worksheets(1).Select 'excel内で1番目のシートを開く

'エクセルシートからのデータ読み込み
NODE = Cells(2, 2) '節点数
ELM = Cells(3, 2) '要素数
HNODE = Cells(4, 2) '発熱する節点数
TNODE = Cells(5, 2) '温度固定箇所の節点数

'マトリックスの配列宣言
ReDim A(NODE, NODE) '熱伝導マトリックス
ReDim B(NODE) '荷重ベクトル
ReDim T(NODE) '温度ベクトル

'計算マトリックスの作成
'熱伝導マトリックス
For i = 1 To ELM
N1 = Cells(2 + i, 4)
N2 = Cells(2 + i, 5)
CELM = Cells(2 + i, 6)
A(N1, N2) = A(N1, N2) - CELM
A(N2, N2) = A(N2, N2) + CELM
A(N2, N1) = A(N2, N1) - CELM
A(N1, N1) = A(N1, N1) + CELM
Next i

'発熱量の処理
For i = 1 To HNODE
num = Cells(2 + i, 7)
B(num) = B(num) + Cells(2 + i, 8)
Next i

'固定温度箇所の処理
For i = 1 To TNODE
num = Cells(2 + i, 9)
A(num, num) = 1
B(num) = Cells(2 + i, 10)

For j = 1 To num - 1
A(num, j) = 0
Next j

For j = num + 1 To NODE
A(num, j) = 0
Next j

Next i

'計算式を解く(熱回路法における熱伝導方程式)

'<前進消去法>
For i = 1 To NODE
For j = i + 1 To NODE
P = A(j, i) / A(i, i)
For k = 1 To NODE
A(j, k) = A(j, k) - P * A(i, k)
Next k
B(j) = B(j) - P * B(i)
Next j
Next i

'<後退代入法>
T(NODE) = B(NODE) / A(NODE, NODE)

For i = NODE - 1 To 1 Step -1 'for文の逆ループ

For j = i + 1 To NODE
B(i) = B(i) - A(i, j) * T(j)
Next j
T(i) = B(i) / A(i, i)

Next i

'計算結果の出力
Worksheets(2).Select 'excel内で2番目のシートを開く
Worksheets(2).Cells.Clear 'excel内で2番目の内容をすべてクリアにする

'出力項目の入力
Cells(2, 2) = "節点番号"
Cells(2, 3) = "温度"

'計算値
For i = 1 To NODE

Cells(i + 2, 2) = i '節点番号の入力
Cells(i + 2, 3) = T(i) '各節点に対する温度の計算値

Next i

End Sub