袁炳志，石文，任書衡

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

0 引言

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

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

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

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

1.1.1 問題分析

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

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

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

表1 符號說明

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

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

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

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

可得虛擬內(nèi)熱源強度為：

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

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

邊界條件為：

1.1.2 離散方程的建立

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

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

也可以推知肋效率表達(dá)式：

1.1.3 求解方法

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

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

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

根據(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ù)可以控制不變，但是當(dāng)N 取到150 及以上時，肋效率的前四位小數(shù)可以控制不變。這時可以認(rèn)為數(shù)值計算的結(jié)果基本與網(wǎng)格數(shù)無關(guān)，即實現(xiàn)了網(wǎng)格無關(guā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（'請輸入?yún)?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）;%設(shè)置無量綱過余溫度場Θ的迭代初值。

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）;%調(diào)用高斯—塞德爾迭代函數(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（['以下為計算結(jié)果：']）

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 求解結(jié)果

表3 具體參數(shù)

（1）溫度分布和傳熱量

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

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

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

圖2

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

圖3 溫度分布曲線

由上述圖表可以看出：

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

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

（2）肋效率

根據(jù)網(wǎng)格無關(guān)性的驗證結(jié)果，取N=150 分別計算了當(dāng)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曲線圖。

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

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

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

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

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

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

其中：

化簡得：

其中：

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

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

為方便MATLAB 求解，將方程組轉(zhuǎn)化為如下形式：

即：

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

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

%參數(shù)設(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 結(jié)果分析

2.2.1 溫度隨時間變化

在MATLAB 中顯式格式采用高斯—賽德爾迭代法進(jìn)行計算，在r2'/r1=4.1，節(jié)點數(shù)N 取100 的情況下對各點處的溫度進(jìn)行求解，同時對比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 中導(dǎo)出溫度分布隨時間變化的視頻，可以看出溫度上升的過程先快后慢。這是由于最一開始溫差較大，傳熱量大，溫度變化大，后來隨著溫度逐漸升高，溫差越來越小，傳熱量減少，溫度場的變化也越來越慢。到50s 以后溫度變化很小，溫度場幾乎不變，此時穩(wěn)態(tài)和非穩(wěn)態(tài)的溫度場分布幾乎一致，相對誤差小于1%。

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

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

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

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

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

2.2.4 傳熱量隨時間變化

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

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

3 結(jié)果分析與總結(jié)

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

圖10 MATLAB傳熱量變化曲線

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

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

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