管式換熱器的VisualBasic模擬
一.引言
列管式換熱器是工業生產特別是熱原料熱設備溫度調節套用普適單元，其模擬原理是基於熱傳導三種基本方式：傳導（conduction），對流（convection）以及輻射（rayonnement）的算法方程來構建的。該模擬可用於在確定出，如口溫度和列管規格的前提下，計算所需長度，或根據已知列管長度，規格，入口溫度推算出口溫度。
二.問題闡述
舉例：氮氣冷卻試驗
列管管路規格：管外徑 De，管壁厚度 e，環境空氣溫度 Tamb，入口溫度Tentrée，出口溫度 Tsortie，氣體流量Debit_M(Kg/h）——直觀條件
管內：強制對流
管外：自然對流
標準大氣壓下空氣及氮氣比熱Cp，導熱係數lamda，黏度µ，密度rho，beta隨溫度變化關係——間接條件須檢索
所用管材外壁黑度與溫度變化關係——間接條件須檢索
目標：
求得各分段管內壁溫度Tpi
求得各分段管外壁溫度Tpe
求得達到各分段所需溫度列管長度
求得達到T sortie 所需總管路長度
三.模擬步驟與算法
1.分段處理
由於氣體物性（比熱Cp，導熱係數lamda，黏度µ，密度rho）在管路中隨溫度變化顯著，為了模擬的精確併兼顧程式計算速度要求，即以總溫度範圍均勻有限分段，取每段平均溫度代表物性變化特徵溫度，把每段求得長度疊加就可得到總管長。
2.單段求解方法
i. 在溫度取值確定以後，氣體物性（比熱Cp，導熱係數lamda，黏度µ，密度rho）可由已知方程得到，氣體在該段流速可由質量流量，密度以及列管規格求得。
ii. 計算管路特性條件需要管路長度L這一必要條件，這裡採用試差法，Visual Basic 中先給L賦初值，添加循環，進行數值比較的語句為 Do until......Loop，解得管內管路特性係數Re，Pr，Gr。
iii. 由Re值範圍確定傳熱方程。其中，對於層流傳熱（Re<2300情況）中需要的管壁溫度Tpi下的氮氣黏度亦可利用試差法解決。最終由求得的Nusselt準數算出對流傳熱係數alpha.
iv. 由每段進，出口溫度求得應散熱量Q0。
以下是該模擬理論核心關係式：
Q0=Q convection d'azote=Q conduction=Q convection extérieure+ Q rayonnement
管內氮氣對流傳熱量 = 管壁傳導熱量 = 管外壁空氣對流傳熱量+管壁輻射傳熱量
v. 由以上關係式，可由對流傳熱方程求出Tpi，由熱傳導方程求出管外壁溫度Tpe。
已知Tpe：
利用管外空氣自然對流方程求得Q convection extérieure
關於輻射傳導，可以根據管束與環境角度情況改變套用方程，本例中設角係數為1（管路被無限大空間包圍情況）
以及最終管外壁總散熱量（Q convection extérieure+ Q rayonnement)
Q0，Q convection extérieure+ Q rayonnement的數值差異就是前面步驟ii 提到的邏輯比較。
vi. 由分段長度L疊加就得到列管總長度。
四.假設與擴展
本模型是列管式換熱器中最簡單的情況，其中分段熱傳導未考慮氮氣傳導方向的管壁傳導傳熱，並且在以下舉例中，管材黑度為定值。
其他複雜列管式換熱器可以根據實際情況改動套用方程或算法，也可有限度地改變已知條件種類，這是其可重建性優勢。
五.程式源碼
Dim i As Integer
Dim j As Integer
Dim k As Integer
Dim SurS As Single '截面積
Dim Nu As Single '紐塞爾準數
Dim Surpi As Single '管內表面積
Dim Surpe As Single '管外表面積
Dim L As Single '分段長度
Dim TpiH As Single '假設管內壁溫度
Dim mu2 As Double '管內壁溫度下氮氣黏度
Public Const n As Integer = 20
Public Const p As Integer = 10
Public Const m As Integer = 1000
Public ro(1 To n) As Double '氮氣密度
Public mu(1 To n) As Double '氮氣黏度
Public Cp(1 To n) As Double '氮氣比熱
Public lambda(1 To n) As Double '氮氣導熱係數
Public beta(1 To n) As Double '溫度倒數
Public T(1 To n) As Single '分段溫度（在Excel 表格中寫入各分段入口溫度及出口溫度中間值，由VB讀入）
Public Re(1 To n) As Single '雷諾準數
Public Pr(1 To n) As Single '普朗特準數
Public Gr(1 To n) As Single '格拉紹夫準數
Public h(1 To n) As Single '對流傳熱係數 同“alpha”
Public V(1 To n) As Single '氮氣流量
Public Tpi(1 To n) As Single '管內壁溫度
Public Tpe(1 To n) As Single '管外壁溫度
Public TAIR(1 To n) As Single '空氣溫度
Public roa(1 To n) As Double '空氣密度
Public mua(1 To n) As Double '空氣黏度
Public Cpa(1 To n) As Double '空氣比熱
Public lambdaa(1 To n) As Double '空氣導熱係數
Public betaa(1 To n) As Double
Public Rea(1 To n) As Single
Public Pra(1 To n) As Single
Public Gra(1 To n) As Single
Public ha(1 To n) As Single
Public Flux1(1 To n) As Double 'Q0
Public Flux2(1 To n) As Double 'Q convection d'azote
Public Flux3(1 To n) As Double
Public ReThP(1 To n) As Single '管壁熱阻
'Constants
Public Const Debit_M As Single = 0.005138 '質量流量
Public Const De As Single = 0.1143 '管外徑
Public Const Di As Single = 0.1053 '管內徑
Public Const delta_T As Single = 76 '氮氣從1600攝氏度降溫到80攝氏度分為20段，每段降溫76度
Public Const lamdap As Single = 15 '管壁導熱係數
Public Const Cpp As Single = 500 '管壁熱容
Public Const Tamb As Single = 20 '環境溫度
Public Const EMIS As Single = 0.6 '管壁黑度
Public Sub evaluation()
SurS = 3.141592653 * (Di / 2) ^ 2
For i = 1 To n
T(i) = Sheets("Simul").Cells(2 + i, 2)
ro(i) = -0.000000001 * T(i) ^ 3 + 0.000003 * T(i) ^ 2 - 0.0026 * T(i) + 1.1279
mu(i) = 0.000000000000007 * T(i) ^ 3 - 0.00000000002 * T(i) ^ 2 + 0.00000004 * T(i) + 0.00002
Cp(i) = -0.0000002 * T(i) ^ 3 + 0.0004 * T(i) ^ 2 + 0.002 * T(i) + 1038.3
lambda(i) = 0.00000000002 * T(i) ^ 3 - 0.00000004 * T(i) ^ 2 + 0.0000735 * T(i) + 0.023934
V(i) = Debit_M / (SurS * ro(i))
beta(i) = 1 / T(i)
Re(i) = ro(i) * Di * V(i) / mu(i)
Pr(i) = Cp(i) * mu(i) / lambda(i)
Flux1(i) = Cp(i) * Debit_M * delta_T
'parametres physiqes et conditions donnees
L = 8
Do Until Abs(Flux1(i) - (Flux2(i) + Flux3(i))) < 60
If Re(i) <= 2100 Then
TpiH = 70
Do Until Abs(Tpi(i) - TpiH) < 20
mu2 = 0.000000000000007 * TpiH ^ 3 - 0.00000000003 * TpiH ^ 2 + 0.00000006 * TpiH + 0.000003
Nu = 1.86 * ((Re(i) * Pr(i)) / (L / Di)) ^ (1 / 3) * ((mu(i) / mu2) ^ 0.14)
h(i) = Nu * lambda(i) / Di
Surpi = 3.141592653 * Di * L
Tpi(i) = T(i) - Flux1(i) / (h(i) * Surpi)
TpiH = TpiH + 5
Loop
ElseIf 2100 < Re(i) And Re(i) < 10000 Then
Nu = (0.023 * Re(i) ^ (4 / 5) * Pr(i) ^ 0.3) * (1 - 600000 / (Re(i) ^ 1.8))
h(i) = Nu * lambda(i) / Di
Surpi = 3.141592653 * Di * L
Tpi(i) = T(i) - Flux1(i) / (h(i) * Surpi)
Else
End If
ReThP(i) = Application.WorksheetFunction.Ln(De / Di) / (2 * 3.141592653 * L * lamdap)
Tpe(i) = Tpi(i) - Flux1(i) * ReThP(i)
Tair(i) = (Tpe(i) - Tamb) / 2 + Tamb
roa(i) = 0.000000000002 * Tair(i) ^ 4 - 0.000000005 * Tair(i) ^ 3 + 0.000006 * Tair(i) ^ 2 - 0.0036 * Tair(i) + 1.2594
Cpa(i) = 0.0000000002 * Tair(i) ^ 4 - 0.0000006 * Tair(i) ^ 3 + 0.0006 * Tair(i) ^ 2 + 0.0189 * Tair(i) + 1002.7
lambdaa(i) = 0.00000000000004 * Tair(i) ^ 4 - 0.00000000007 * Tair(i) ^ 3 + 0.00000003 * Tair(i) ^ 2 + 0.00006 * Tair(i) + 0.0255
mua(i) = 9E-18 * Tair(i) ^ 4 - 0.00000000000002 * Tair(i) ^ 3 - 0.000000000001 * Tair(i) ^ 2 + 0.00000004 * Tair(i) + 0.00002
betaa(i) = 1 / Tair(i)
Gra(i) = betaa(i) * 9.81 * (Tpe(i) - Tamb) * De ^ 3 * roa(i) / mua(i) ^ 2
Pra(i) = Cpa(i) * mua(i) / lambdaa(i)
Surpe = 3.141592653 * De * L
If 1 < Gra(i) * Pra(i) < 10 ^ 4 Then
Nu2 = 1.09 * (Gra(i) * Pra(i)) ^ 0.2
ha(i) = Nu2 * lambdaa(i) / De
Flux2(i) = ha(i) * (Tpe(i) - Tamb) * Surpe
ElseIf 10 ^ 4 < Gra(i) * Pra(i) < 10 ^ 9 Then
Nu2 = 0.53 * (Gra(i) * Pra(i)) ^ 0.25
ha(i) = Nu2 * lambdaa(i) / De
Flux2(i) = ha(i) * (Tpe(i) - Tamb) * Surpe
ElseIf 10 ^ 9 < Gra(i) * Pra(i) < 10 ^ 12 Then
Nu2 = 0.13 * (Gra(i) * Pra(i)) ^ 0.33
ha(i) = Nu2 * lambdaa(i) / De
Flux2(i) = ha(i) * (Tpe(i) - Tamb) * Surpe
End If
Flux3(i) = emis * 5.669 * 1 * Surpe * (((Tpe(i) + 273) / 100) ^ 4 - ((Tamb + 273) / 100) ^ 4)
L = L - 0.01
Loop
ThisWorkbook.Worksheets("Simul").Cells(2 + i, 4) = L ' 結果寫入與定位
Next
End Sub
