環(huán)形肋片導熱問題的數(shù)值模擬

袁炳志，石文，任書衡

（華北電力大學能源動力與機械工程學院，保定071000）

0 引言

數(shù)值模擬是計算機技術在工程上最為廣泛的應用之一，其依靠電子計算機，結合有限元的概念，通過數(shù)值計算和圖像顯示的方法，對工程問題和物理問題的研究與求解。在換熱設備中，通常使用肋片增大傳熱面積、降低熱阻，從而提高傳熱效率。環(huán)形肋片可以布置在圓管上，有著優(yōu)良的導熱性能，但在分析求解中，由于環(huán)肋與普通直肋的區(qū)別，不能直接進行理論計算，需要借助數(shù)值模擬求解[1]。

本文基于顯式和隱式兩種離散方式對傳熱方程進行離散，利用MATLAB 通過高斯-賽德爾迭代和雅各比迭代兩種迭代方式計算了環(huán)肋的穩(wěn)態(tài)傳熱和瞬態(tài)傳熱的問題。通過比較肋效率、傳熱量和溫度場，對離散方式和迭代方式進行比較和分析。

1 穩(wěn)態(tài)問題數(shù)值求解

1.1 數(shù)學模型的建立

1.1.1 問題分析

已知肋根處的溫度t0高于肋片附近的流體溫度tf，熱量由肋根處向肋端處導熱，同時肋片表面和周圍空氣存在對流換熱。

圖1 環(huán)肋的物理模型

嚴格來說，肋片的導熱問題為圓柱坐標系下三維非穩(wěn)態(tài)導熱問題，肋片內部的溫度場為三維的溫度場t=f(r,?,z,τ)，則其一般形式的導熱微分方程為[2]：

表1 符號說明

要分析一般情況下環(huán)肋的導熱問題，需要對該式進行簡化：

（4）由于通過肋端的散熱量很小，肋端可看作絕熱。但是為了減少計算帶來的誤差，對肋高進行修正。修正后肋高

（5）肋片表面的熱量傳遞不能忽略，把表面?zhèn)鬟f的熱量折合成虛擬內熱源Φ?：

導熱微分方程中的內熱源項由微元段表面的對流傳熱量折算得到[3]。沿肋高方向任取dr 作為研究對象，則該微元段表面總散熱量可以表示為：

可得虛擬內熱源強度為：

在工程應用中，Bi＜0.1 則認為肋片的溫度分布符合一維假設。計算可得Bi=0.0005 滿足所需條件。

通過以上的假設分析，可以把導熱微分方程簡化為圓柱坐標系下一維有內熱源的穩(wěn)態(tài)導熱問題。則其導熱微分方程為：

邊界條件為：

1.1.2 離散方程的建立

所研究問題的導熱微分方程是一個二階非齊次常微分方程，為轉化成齊次方程便于求解，引入無量綱過余溫度和無量綱半徑；為便于采用特征方程求解，令。導熱微分方程可進一步化簡為：

對環(huán)肋進行空間區(qū)域離散化，取空間步長Δr，在半徑r 方向上取N 個節(jié)點，節(jié)點編號為n=1，2，3，…，N-1，N。離散化過程中采用二階精度泰勒展開的中心差分格式，建立離散化方程：

也可以推知肋效率表達式：

1.1.3 求解方法

將離散方程簡化整理，可得：

為方便編程求解，將化簡后的離散方程式表示為Aθ=b的形式，其中A 為N×N的系數(shù)矩陣，θ是列向量，表示無量綱過余溫度。在MATLAB 里定義系數(shù)矩陣后，對特殊元素進行賦值，使其滿足離散化方程。并調用Gauss-Seidel 迭代函數(shù)進行數(shù)值計算，求得無量綱過余溫度場，最終得到溫度場。

1.1.4 網格無關性驗證

根據(jù)MATLAB 模擬的程序，在控制r2'/r1以及m 不變的情況下，任意輸入節(jié)點數(shù)N，可以得到不同的肋效率η的值。下表列出了在r2'/r1=2 且m=1.5 的情況下，節(jié)點數(shù)N 分別取20、30、45、60、100、150、200 時肋效率的值：

表2 肋效率的值

由表2 可知，在不同的節(jié)點數(shù)N 的情況下，肋效率的前兩位小數(shù)可以控制不變，但是當N 取到150 及以上時，肋效率的前四位小數(shù)可以控制不變。這時可以認為數(shù)值計算的結果基本與網格數(shù)無關，即實現(xiàn)了網格無關性驗證。

1.2 模型的求解與分析

1.2.1 求解代碼

輸入：r2'/r1、m、節(jié)點數(shù)N

輸出：溫度場矩陣t、肋效率efficiency、傳熱量fi1

k=input（'請輸入r2/r1 的值.'）;

disp（['r2/r1=',num2str（k）]）;

N=input（'請輸入節(jié)點數(shù)N 的值.'）;

disp（['節(jié)點數(shù)N=',num2str（N）]）;

m=input（'請輸入參數(shù)m 的值.'）;

disp（['參數(shù)m=',num2str（m）]）;

dR=1/（N-1）;%定義無量綱半徑的步長

R=zeros（N,1）;%無量綱半徑矩陣

for i=1:1:N

R（i）=1/（k-1）+（i-1）/（N-1）;%為無量綱半徑賦值

end

b=zeros（N,1）;%定義一個N 行1 列的零矩陣，作為迭代求解時等式右邊的矩陣。

b（1）=1;%邊界條件1：Θ1=1。

theta0=ones（N,1）;%設置無量綱過余溫度場Θ的迭代初值。

A=zeros（N,N）;%定義系數(shù)矩陣A。接下來對特殊元素賦值。

A（1,1）=1;%同樣是為了滿足邊界條件1：Θ1=1。

for i=2:1:N-1

A（i,i-1）=（1/（dR^2））-（1/（2*dR*R（i）））;

A（i,i）=-（2/（dR^2））-2*m^2;

A（i,i+1）=（1/（dR^2））+（1/（2*dR*R（i）））;

end%此循環(huán)為系數(shù)矩陣A 賦值。

A（N,N-1）=-1;

A（N,N）=1;%以上兩個賦值是為了滿足邊界條件2：ΘN-ΘN-1=0。（即頂端絕熱條件）

theta=gaussseidel（A,b,theta0）;%調用高斯—塞德爾迭代函數(shù)。

%以下開始計算肋效率

a1=0;

a2=0;

for i=2:1:N-1

a1=a1+R（i）*theta（i）;

a2=a2+R（i）;

end

a1=a1+R（1）/2*theta（1）+R（N）/2*theta（N）;

a2=a2+R（1）/2+R（N）/2;

efficiency=a1/a2;%肋效率

disp（['以下為計算結果：']）

disp（['肋效率為：',num2str（efficiency）]）;

%以下為計算傳熱量和溫度場

h=50;

r1=0.01;

r2=0.04;

t0=100;

tf=20;

fi1=a1*2*pi*（r2-r1）^2*dR*（t0-tf）*h;

t=theta*（t0-tf）+tf+273.15;

1.2.2 求解結果

表3 具體參數(shù)

（1）溫度分布和傳熱量

在MATLAB 中利用高斯賽德爾迭代法進行計算，在r2'/r1=4.1 以及m=0.49015，節(jié)點數(shù)N 取150 的情況下對各點處的溫度進行求解，得到肋片傳熱量為15.135W，同時對比Fluent 各點處的溫度分布如表4。

將MATLAB 中獲得溫度場導出為dat 文件，導入Tecplot 里分別繪制出MATLAB 以及Fluent 的溫度分布圖像如圖2a、圖2b 所示。

表4 溫度在不同節(jié)點位置的分布情況

圖2

得到整體的溫度分布曲線為圖3。

圖3 溫度分布曲線

由上述圖表可以看出：

隨節(jié)點坐標的增加，越靠近肋端處溫度越低，并且在傳熱過程中溫度變化的斜率逐漸減小，說明溫差降低使傳熱速率降低，傳熱量逐漸減少，與定性分析的結果相符。因此推測：在一定范圍內，增加肋片高度可以使肋端溫度進一步降低，增加肋片的傳熱量，但同時肋效率也降低。

MATLAB 和Fluent 得到的溫度分布情況十分接近。兩種方式獲得的肋端處溫度分別下降到了350.0562K 和357.1077K，相對誤差約為2%，誤差較小，結果較為精確，與定性分析的結果相符。進一步證明了MATLAB 模擬結果的可信度與精確度。

（2）肋效率

根據(jù)網格無關性的驗證結果，取N=150 分別計算了當r2'/r1為2、3、4.1、5，m 為0.1、0.5、1、1.5、2、2.5 時肋效率的值：

表5 不同r2'/r1 和m 情況下肋效率的值

①根據(jù)表中數(shù)據(jù)，繪制r2'/r1=2 的圖線與教材中給出的圖線做對比。

圖4 r2'/r1=2 時的η-m 曲線圖

曲線actual（教材中的圖線）與simulated（模擬的曲線）高度擬合，說明上述編譯的解決方案精確度很高，可以用來求解所述案例。

②根據(jù)表中數(shù)據(jù)繪制r2'/r1=4.1 時的η-m曲線圖。

肋效率的物理意義為實際散熱量與假設整個肋表面處于肋基溫度下的散熱量之比。根據(jù)r2'/r1=4.1的圖線可得：m 越小即肋高越小，肋片表面的平均溫度越接近肋基溫度，肋效率就會越大。對比所有的曲線可得：r2'/r1的值越小，在m 相同的情況下肋效率的值越大。所以，適當減小。r2'/r1和m 可以增大肋效率，但會減小散熱量。

圖5 不同r2'/r1 時的η-m 曲線圖

2 非穩(wěn)態(tài)問題數(shù)值求解

2.1 數(shù)學模型的建立與求解方法

由于對于穩(wěn)態(tài)問題的求解，采取泰勒級數(shù)展開法進行離散化，為了與穩(wěn)態(tài)的結果相印證對比，并比較兩種離散方式的差異，采用控制容積熱平衡法來對非穩(wěn)態(tài)問題進行離散化處理。該方法是根據(jù)每個節(jié)點所代表的控制容積的能量守恒關系，得到該節(jié)點與相鄰節(jié)點溫度間的關系式。在離散方程中，空間步長為Δr，時間步長為Δτ。為滿足穩(wěn)定性條件，時間步長取0.0005s。

（1）采用顯式格式離散方程得：

其中：

化簡得：

其中：

至此得到非穩(wěn)態(tài)問題的顯式離散方程和相應邊界條件。采用高斯—賽德爾方法對方程進行迭代處理，可以求得非穩(wěn)態(tài)問題的溫度場變化。

（2）同理可求得隱式格式離散方程

為方便MATLAB 求解，將方程組轉化為如下形式：

即：

這樣一來，就把求解溫度分布的問題轉化為求解X 列向量的問題，便于編程計算。在穩(wěn)態(tài)問題的求解中，運用高斯賽德爾迭代，運用兩種迭代方法求解方程，并比較溫度場的差異。對于高斯賽德爾迭代，利用先前構造的高斯賽德爾迭代函數(shù)即可求解。

對于雅各比迭代，由于初始溫度已知，則可求出等式右邊列向量B，利用X=A/B 可以求得下一時層中X的值，將其賦值給B，形成迭代計算，最終求得溫度場。

%參數(shù)設定

h=50;

a=0.00004115;

N=10;

r1=0.01;

r2=0.041;

tf=20;

delta=0.002;

%定義空間步長，并通過迭代得出半徑矩陣

dr=0.031/（N-1）;

r=zeros（N,1）;

for i=1:1:N

r（i）=r1+（i-1）*dr;

end

dtime=1/20;%定義時間步長

tt=20*ones（N,1）;%初始溫度

tt（1）=100;%邊界條件

t=ones（N,1）;

counter=0;

while counter＜=1000000

t（1）=100;

counter=counter+1;

for i=2:1:N-1

t（i）=（a*（r（i）-dr/2）*dtime/（r（i）*dr*dr））*t（i-1）+（1-（a*（r（i）-dr/2）*...

dtime/（r（i）*dr*dr））-（a*（r（i）+dr/2）*dtime/（r（i）*dr*dr））-（2*h*dtime/（0.002*270...0*900）））*tt（i）+（a*（r（i）+dr/2）*dtime/（r（i）*dr*dr））*tt（i+1）+（（2*h*dtime/（0.002...*2700*900））*tf）;

end

t（N）=（a*（r2-dr/2）*dtime/（dr*r2*（dr/2）））*t（N-1）+（1-（a*（r2-dr/2）*dtime/...

（dr*r2*（dr/2）））-（2*h*r2*dtime/（0.002*2700*900）））*tt（N）+（2*h*r2*dtime/（0.002*...2700*900））*tf;

tt=t;

end

t=t+273.15;

2.2 結果分析

2.2.1 溫度隨時間變化

在MATLAB 中顯式格式采用高斯—賽德爾迭代法進行計算，在r2'/r1=4.1，節(jié)點數(shù)N 取100 的情況下對各點處的溫度進行求解，同時對比Fluent 各點處的溫度分布可得溫度隨時間的變化圖像（圖6）。

圖6

從溫度曲線可以看出靠近肋基的地方溫度下降的較快，靠近肋端的地方溫度變化較為平緩，越靠近肋端溫度越低，說明在溫差大的情況下，傳熱的效果強。

顯式格式與Fluent 對比如圖7，以5s,25s,50s 時刻為例，溫度場的相對誤差均在1%以下，誤差較小，并且隨著時間的推移，各個節(jié)點的溫度不斷升高，直到50s左右溫度分布基本不變，非穩(wěn)態(tài)過程趨于穩(wěn)定，此時溫度分布情況與穩(wěn)態(tài)的溫度分布情況近乎一致，肋端處的溫度都約為350℃。

圖7 不同時刻MATLAB與Fluent溫度分布曲線

2.2.2 顯式差分、隱式差分比較

在求解顯式差分格式的問題中，從MATLAB 中導出溫度分布隨時間變化的視頻，可以看出溫度上升的過程先快后慢。這是由于最一開始溫差較大，傳熱量大，溫度變化大，后來隨著溫度逐漸升高，溫差越來越小，傳熱量減少，溫度場的變化也越來越慢。到50s 以后溫度變化很小，溫度場幾乎不變，此時穩(wěn)態(tài)和非穩(wěn)態(tài)的溫度場分布幾乎一致，相對誤差小于1%。

溫度分布在肋基處溫度最高，越靠近肋端處溫度越低，并且在顯隱式情況下求解得到的結果十分接近。而顯式與隱式的差別就是對時間的向前向后差分，所以對同一問題所得結果極為接近，這說明了上述模型符合預期，能達到模擬該案例的水平。

圖8 50秒時，顯式和隱式差分溫度分布對比圖

2.2.3 雅各比迭代和高斯賽德爾迭代比較

為比較兩種迭代方法，取10s 時的溫度場進行對比，如圖9 所示。可以看出采用兩種迭代方法對該問題的影響很小。

圖9 10秒時，雅各比迭代和高斯賽德爾迭代溫度分布對比圖

2.2.4 傳熱量隨時間變化

圖9 為MATLAB 中傳熱量隨時間變化的曲線，從圖中可看出開始階段溫差大，導致熱流量變化大，之后溫差逐漸減小，傳熱量趨于平緩。為對比驗證，做一條傳熱量為15.135W 的直線，50 秒以后非穩(wěn)態(tài)的傳熱量與穩(wěn)態(tài)的傳熱量完全重合，相對誤差為0.2%。

在50s 之前傳熱量不斷增加，起初傳熱量變化較大，之后逐漸平穩(wěn)，趨近于一個穩(wěn)定值，可以認為近乎達到了穩(wěn)態(tài)。通過觀察溫度變化曲線可知，初始時刻溫度變化較大，50s 后溫度幾乎不變，可以看做穩(wěn)態(tài)情況下的溫度分布。溫度分布情況與穩(wěn)態(tài)時求解結果相近，肋端處溫度均為350℃左右。同時在非穩(wěn)態(tài)傳熱過程中，由于肋基溫度不變，熱量逐漸從肋基傳到肋端，并散失到外界，離肋基越遠，溫度越低，溫差越小，溫度分布曲線趨于平緩，符合理論分析。

3 結果分析與總結

對于穩(wěn)態(tài)問題，通過簡化和合理假設，先推導出了環(huán)肋一維有內熱源的導熱微分方程，進一步進行無量綱化并采用泰勒級數(shù)展開法中心差分格式進行離散，利用MATLAB 求解溫度場和傳熱量，并計算相應的肋效率。對于非穩(wěn)態(tài)問題，利用控制容積熱平衡法進行空間區(qū)域離散，并分別運用顯式差分和隱式差分對時間區(qū)域離散，對于隱式差分分別采用高斯—賽德爾迭代和雅各比兩種迭代方式進行計算和比較。結果總結如下：

圖10 MATLAB傳熱量變化曲線

●MATLAB、Fluent 模擬所得溫度場分布一致，肋端溫度在350K 左右，相對誤差不超過2%?？拷呋帨囟认陆递^快，靠近肋端處溫度變化較為平緩，越靠近肋端溫度越低。增大肋高會增大散熱量，降低肋效率。

●非穩(wěn)態(tài)狀況，MATLAB 所得顯式與隱式結果一致，不同節(jié)點在5s,25s,50s 的溫度數(shù)值與Fluent 模擬結果相對誤差不超過1%。取10s 時溫度場對比，雅各比迭代和高斯賽德爾迭代結果也無差別。

●非穩(wěn)態(tài)下傳熱量隨時間增加，增長速率逐漸降低，50s 左右后與穩(wěn)態(tài)散熱量保持一致，相對誤差為0.2%。