劉天昀，劉文波

（1.華北電力大學(xué)能源動力與機(jī)械工程學(xué)院，保定071003；2.廣東生態(tài)工程職業(yè)學(xué)院信息工程系，廣州510520）

0 引言

肋片作為二次傳熱面，擴(kuò)大傳熱面積，并促進(jìn)流體介質(zhì)擾動來強(qiáng)化換熱[1]。在各大領(lǐng)域中，肋片得到廣泛的應(yīng)用。無論是在空調(diào)的冷凝器，還是在汽車發(fā)動機(jī)的散熱裝置、計算機(jī)CPU的散熱器中，都能看到它的身影。本文通過數(shù)值模擬的方法分別研究了等截面圓柱形肋片在穩(wěn)態(tài)和非穩(wěn)態(tài)情況下，肋片內(nèi)部溫度分布，通過肋片的傳熱量以及肋片效率。

1 基于MATLAB穩(wěn)態(tài)導(dǎo)熱問題數(shù)值模擬

1.1 計算區(qū)域的離散

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

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

1.2 建立節(jié)點方程

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

圖1 等截面圓柱形肋模型

沿中心軸節(jié)點：

肋端節(jié)點：

圓周與肋端交界處節(jié)點：

內(nèi)部節(jié)點：

肋根節(jié)點：

圓周處節(jié)點：

1.3 MATLAB計算結(jié)果分析

（1）網(wǎng)格無關(guān)性驗證

表1 節(jié)點數(shù)與散熱量的關(guān)系

通過上表可見，當(dāng)沿高度方向節(jié)點數(shù)20以及沿半徑方向節(jié)點數(shù)為20，相對誤差為0.022%，因此本次課題將采用高度方向節(jié)點數(shù)20以及半徑方向節(jié)點數(shù)20進(jìn)行相關(guān)計算。

（2）溫度場分析

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

圖2 圓柱體溫度分布

圖3 肋端溫度分布

（3）肋片散熱量分析

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

圖4 散熱量計算程序

（4）肋片效率

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

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

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

2 基于MATLAB非穩(wěn)態(tài)導(dǎo)熱問題數(shù)值模擬

本課題求解非穩(wěn)態(tài)導(dǎo)熱問題，計算區(qū)域離散方法以及建立節(jié)點方程方法與穩(wěn)態(tài)時相似，在此不再贅述，只給出相應(yīng)結(jié)果。

2.1 建立節(jié)點方程

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

沿中心軸節(jié)點：

肋端節(jié)點：

圓周與肋端交界處節(jié)點：

內(nèi)部節(jié)點：

肋根節(jié)點：

圓周處節(jié)點：

2.2 肋片溫度

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

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

圖5 1000步迭代

圖6 4000步迭代

圖7 8000步迭代

圖8 10000步迭代

從迭代時間步長以及實際物理經(jīng)過時間分析，顯示格式計算非穩(wěn)態(tài)問題時需考慮穩(wěn)定性條件，在某些情況下出現(xiàn)時間步長過短，迭代次數(shù)過長的情況發(fā)生。因此，隱式格式計算非穩(wěn)態(tài)情況適用范圍更廣。

2.3 肋片散熱量

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

表2 散熱量隨時間變化

由于計算機(jī)性能局限，無法得到肋片達(dá)到穩(wěn)態(tài)時變化情況，因此對散熱量隨時間變化做線性預(yù)測，如圖9所示，通過冪函數(shù)預(yù)測曲線，當(dāng)穩(wěn)態(tài)時散熱量6.6893W，所需時間約為85.41s，迭代次數(shù)約為21步。

散熱量與時間關(guān)系預(yù)測曲線：

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

3 基于商業(yè)軟件Fluent數(shù)值模擬

3.1 幾何模型建立

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

3.2 網(wǎng)格劃分

在確定網(wǎng)格數(shù)量時，既需要滿足計算精度要求，又要避免計算機(jī)運算負(fù)擔(dān)過重[3]。以不同的網(wǎng)格劃分?jǐn)?shù)量為變量，以端部溫度為對比參數(shù)進(jìn)行網(wǎng)格無關(guān)性分析。

3.3 網(wǎng)格無關(guān)性驗證

通過表3可見當(dāng)網(wǎng)格規(guī)模于30萬以上時，肋端溫度相對誤差小于1%，所以取30萬以上規(guī)模的網(wǎng)格進(jìn)行數(shù)值計算。

表3 網(wǎng)格總數(shù)與肋端溫度的關(guān)系

3.4 物理模型構(gòu)建

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

圖10 導(dǎo)入Fluent后的網(wǎng)格文件

3.5 數(shù)值模擬分析

在Fluent中設(shè)置好邊界條件后，分別對穩(wěn)態(tài)和非穩(wěn)態(tài)條件下等截面圓柱形肋片導(dǎo)熱問題進(jìn)行迭代計算，然后得到殘差曲線圖、溫度分布云圖、等溫線圖等相關(guān)數(shù)據(jù)。

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

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

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

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

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

不同時刻下，還可以通過report讀取對應(yīng)時刻通過肋片的熱流量。通過讀取可得隨著時間的增加，通過肋片的熱流量不斷減小，最終達(dá)到穩(wěn)態(tài)時熱量為6.668W，達(dá)到穩(wěn)態(tài)時間所需時間約為116s左右。

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

4 Fluent數(shù)值模擬與MATLAB模擬對比分析

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

4.1 穩(wěn)態(tài)情況

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

4.2 非穩(wěn)態(tài)情況

非穩(wěn)態(tài)時，通過對達(dá)到穩(wěn)態(tài)情況所需時間分析，兩種模擬方式差別較大，其中MATLAB模擬所需時間為85.61s而Fluent模擬所需時間為116s。

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

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

5 結(jié)語

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

（2）非穩(wěn)態(tài)情況下，肋片整體溫度隨時間而升高。到一定時刻后，溫度場不再隨時間變化，達(dá)到穩(wěn)態(tài)情況。

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

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