等截面圓柱形肋導熱數值模擬研究

劉天昀，劉文波

（1.華北電力大學能源動力與機械工程學院，保定071003；2.廣東生態工程職業學院信息工程系，廣州510520）

0 引言

肋片作為二次傳熱面，擴大傳熱面積，并促進流體介質擾動來強化換熱[1]。在各大領域中，肋片得到廣泛的應用。無論是在空調的冷凝器，還是在汽車發動機的散熱裝置、計算機CPU的散熱器中，都能看到它的身影。本文通過數值模擬的方法分別研究了等截面圓柱形肋片在穩態和非穩態情況下，肋片內部溫度分布，通過肋片的傳熱量以及肋片效率。

1 基于MATLAB穩態導熱問題數值模擬

1.1 計算區域的離散

本課題所研究的等截面圓柱形肋高為0.06m，橫截面直徑為0.01m。材料的比熱為900J/（kg·K），密度為2700kg/m3，導熱系數大小為50W/（m·K）。其中，肋根溫度為100℃，周圍空氣溫度為20℃，不考慮肋片與周圍可見物體的輻射傳熱。模型如圖1所示。

根據等截面圓肋特點，建立沿徑向以及沿軸向二維模型。對該區域離散時，在r方向共有M個節點，空間步長取dr；在長度方向有N個節點，空間步長取dz。

1.2 建立節點方程

本課題根據控制容積熱平衡法分別對底面節點、肋端面節點、內部節點以及邊界節點建立離散方程。

圖1 等截面圓柱形肋模型

沿中心軸節點：

肋端節點：

圓周與肋端交界處節點：

內部節點：

肋根節點：

圓周處節點：

1.3 MATLAB計算結果分析

（1）網格無關性驗證

表1 節點數與散熱量的關系

通過上表可見，當沿高度方向節點數20以及沿半徑方向節點數為20，相對誤差為0.022%，因此本次課題將采用高度方向節點數20以及半徑方向節點數20進行相關計算。

（2）溫度場分析

輸入相關節點數后，可得到等面積圓柱肋片溫度場。如圖所示，對溫度場分析可知肋片溫度主要沿肋高方向變化，沿半徑溫度幾乎不變，肋端溫度為86.48℃。

圖2 圓柱體溫度分布

圖3 肋端溫度分布

（3）肋片散熱量分析

肋片散熱量是肋片邊界與周圍流體對流換熱散熱量。根據肋片特點，將肋片分為三部分進行計算散熱量。第一部分，肋片圓周但不含與肋端交界邊與周圍流體對流散熱量；第二部分，肋端但不含與肋片圓周交界處與周圍流體對流散熱量；第三部分，肋片圓周處與肋端交界處離散區域與周圍流體對流散熱量[1]。根據分塊計算思想，相應程序如圖4所示，其中Q1代表第二部分散熱量，Q2代表第一部分散熱量，q（M,N）代表第三部分散熱量，計算結果為6.6893W。

圖4 散熱量計算程序

（4）肋片效率

在基礎的散熱表面上敷設肋片，主要是想通過增加面積來增大對流傳熱的熱流量。但是肋片表面的溫度是隨肋片的高度變化的，所以肋片的形狀和尺寸（特別是高度）對增強換熱的效果有重要的影響。在計算肋片散熱量時，為簡化計算提出肋片效率ηf一概念，其定義式為：

ηf=（肋片實際散熱量）/（假設整個肋片處于肋根溫度下的散熱量）

根據數值求解散熱量，ηf=85.21%，本課題采用數值求解方法采用二維模型，而理論求解采用一維模型，這是產生誤差的主要原因[2]。

2 基于MATLAB非穩態導熱問題數值模擬

本課題求解非穩態導熱問題，計算區域離散方法以及建立節點方程方法與穩態時相似，在此不再贅述，只給出相應結果。

2.1 建立節點方程

本課題采用顯示格式建立節點方程。

沿中心軸節點：

肋端節點：

圓周與肋端交界處節點：

內部節點：

肋根節點：

圓周處節點：

2.2 肋片溫度

本課題采用顯示格式進行相關計算，通過初始化，將所有節點溫度都設為20℃。該計算中，實際物理所經過的時間等于時間步長數乘以時間步長，同時為了避免迭代結果離散情況發生，時間步長需滿足穩定性條件，本課題根據肋片溫度以及幾何條件，設定最大時間步長為0.004s，迭代次數10000步。

通過對迭代次數為1000、4000、8000以及10000溫度分布圖分析，肋片溫度變化特點：溫度變化位置從肋根開始逐漸至肋端處；隨著時間延長，肋片溫度逐漸不明顯。

圖5 1000步迭代

圖6 4000步迭代

圖7 8000步迭代

圖8 10000步迭代

從迭代時間步長以及實際物理經過時間分析，顯示格式計算非穩態問題時需考慮穩定性條件，在某些情況下出現時間步長過短，迭代次數過長的情況發生。因此，隱式格式計算非穩態情況適用范圍更廣。

2.3 肋片散熱量

非穩態肋片散熱量計算方法與穩態情況相似。肋片散熱量是肋片邊界與周圍流體對流換熱散熱量。其散熱量變化如表2所示，散熱量隨時間變化逐漸層增加，但由于肋片整體溫度逐漸增加，與周圍流體換熱不斷加強，散熱量隨時間變化斜率逐漸增大。

表2 散熱量隨時間變化

由于計算機性能局限，無法得到肋片達到穩態時變化情況，因此對散熱量隨時間變化做線性預測，如圖9所示，通過冪函數預測曲線，當穩態時散熱量6.6893W，所需時間約為85.41s，迭代次數約為21步。

散熱量與時間關系預測曲線：

圖9 散熱量隨時間變化線性曲線

3 基于商業軟件Fluent數值模擬

3.1 幾何模型建立

本課題所研究的等截面圓柱形肋高為0.06m，橫截面直徑為0.01m。

3.2 網格劃分

在確定網格數量時，既需要滿足計算精度要求，又要避免計算機運算負擔過重[3]。以不同的網格劃分數量為變量，以端部溫度為對比參數進行網格無關性分析。

3.3 網格無關性驗證

通過表3可見當網格規模于30萬以上時，肋端溫度相對誤差小于1%，所以取30萬以上規模的網格進行數值計算。

表3 網格總數與肋端溫度的關系

3.4 物理模型構建

本次課題中，使用Gambit 2.4.6建模.使用直接生成體的方法繪制出高為0.06m，橫截面直徑為0.01m的圓柱形。在本例中利用體生成的方式直接生成網格，網格類型設置為map，選取合適的interval size，網格總數要30萬以上。網格劃分完畢后，需要對邊界條件進行設定。在本例中，所有壁面設置為wall。但是需要注意的是，肋端、肋根以及側面的邊界條件需要分開定義。建立好后的幾何模型如圖10所示。

圖10 導入Fluent后的網格文件

3.5 數值模擬分析

在Fluent中設置好邊界條件后，分別對穩態和非穩態條件下等截面圓柱形肋片導熱問題進行迭代計算，然后得到殘差曲線圖、溫度分布云圖、等溫線圖等相關數據。

穩態狀況下時，通過做圓柱體縱向切片觀察溫度分布云圖，發現在距離肋根0.04m以內的范圍內基本只沿肋高方向變化，徑向的溫度梯度較小。而在端部附近的區域，徑向的溫度梯度不可忽略，如圖11所示（右邊為肋根，左邊為肋端）。通過畫plot圖，可以發現隨著與肋根距離的增加，溫度不斷下降，溫度梯度不斷減小，最終在肋端處，溫度達到最小值359.3K即86.15℃，如圖12所示。通過report可以讀取通過肋片的熱流量為6.668 W。

圖11 圓柱形肋溫度 分布云圖

圖12 沿肋高方向的 溫度曲線圖

在非穩態狀況下，改變General選項中Time的設置，選中Time中的Transient選項，選擇合適的時間步長，初始溫度設為20℃，其余邊界條件保持不變。使用同樣的方式做不同時刻下的縱向切片的溫度分布云圖，觀察溫度場隨時間的變化情況，如圖13所示（右邊為肋根，左邊為肋端）。

通過做不同時刻下的縱向切片的溫度分布云圖，觀察到肋片整體溫度隨時間而升高。到一定時刻后，溫度場不再隨時間變化，沿著肋高方向溫度不斷下降，溫度梯度不斷減小。

不同時刻下，還可以通過report讀取對應時刻通過肋片的熱流量。通過讀取可得隨著時間的增加，通過肋片的熱流量不斷減小，最終達到穩態時熱量為6.668W，達到穩態時間所需時間約為116s左右。

圖13 不同時刻下溫度分布云圖

4 Fluent數值模擬與MATLAB模擬對比分析

本課題分別使用MATLAB以及Fluent等軟件對等截面圓柱肋片進行數值模擬工作，為使模擬工作更加精確，對兩種方式得到相關結論進行比較并分析誤差產生原因。

4.1 穩態情況

溫度場分析，如表4所示，通過對比兩種方法得出肋端溫度以及散熱量，兩者相差較小。穩態情況下，兩者模擬結論差別較小。這主要是因為：兩者都是數值解法，基本思想都是把原來在空間與時間域內連續的物理量場用有限個離散點上的值的集合來代替[4]，數值解法的精確度與網格的密度及數量密切相關，網格劃分越小，結果越精確。在網格密度都達到一定標準時，使用編程和商業軟件進行數值模擬得到的結果基本一致，誤差幾乎可以忽略不計。

4.2 非穩態情況

非穩態時，通過對達到穩態情況所需時間分析，兩種模擬方式差別較大，其中MATLAB模擬所需時間為85.61s而Fluent模擬所需時間為116s。

在非穩態時通過商業軟件得到的結果與編程得到的略微存在一定偏差，這主要因為：與編程數值模擬不同的是，Fluent是基于有限元法的商業軟件。有限元分析是用較簡單的問題代替復雜問題后再求解，它將求解域看成是由許多稱為有限元的小的互連子域組成，對每一單元假定一個合適的近似解，然后推導求解這個域總的滿足條件，從而得到問題的解[5]。與編程使用的有限差分法相比，Fluent使用的算法對于不同的算例的求解都有更好的收斂性、穩定性和精度。

使用基于有限差分法的MATLAB編譯處理非穩態問題時，保持空間步長不變，增大時間步長，此時傅里葉數會增大，甚至不再滿足節點離散方程的穩定性條件，因此其結果將進一步發散，與Fluent相比沒那么精確。再者由于計算機性能限制以及程序算法的不足，使用MATLAB方法得到肋片達到穩態所需時間經過數據擬合得到，因此此種方式本就會帶來一定誤差。

5 結語

（1）通過數值模擬，觀察到穩態情況下肋片溫度沿肋高方向不斷下降，沿半徑溫度幾乎不變，肋效率為85.21%。

（2）非穩態情況下，肋片整體溫度隨時間而升高。到一定時刻后，溫度場不再隨時間變化，達到穩態情況。

（3）分別通過MATLAB編譯計算和Fluent的迭代計算，觀察到兩者在穩態情況下得到的結果差異很小，非穩態時有一定誤差且Fluent所得結果更精確。

（4）非穩態情況下，選擇顯式格式進行計算時，時間步長需要滿足穩定性條件，本課題中時間步長為0.4ms。因此，當需得到較完整物理變化過程，迭代步數將特別大，本課題需迭代21萬步左右，對計算機性能有較大考驗。因此，對非穩態工況計算建議采用隱式格式。