王迪昌,劉澤遠,劉捷?,盧文強

(1 中國科學(xué)院大學(xué)工程科學(xué)學(xué)院, 北京 100049; 2 浙江大華技術(shù)股份有限公司, 杭州 310053)

自然對流換熱是一種生活和工程應(yīng)用中常見的對流換熱方式,其驅(qū)動原理與強迫對流不同,流體運動僅由浮力引起,而浮力通常是由溫度變化引起的密度梯度和重力共同作用的結(jié)果[1]。盡管自然對流的換熱系數(shù)小于強迫對流,但是簡單的驅(qū)動原理和低成本等優(yōu)點使其在許多工業(yè)領(lǐng)域或設(shè)備中得到廣泛使用,如球形電子設(shè)備、循環(huán)流化床、固定床氣化爐[2-3]和核工程設(shè)備中。后者包括加速器驅(qū)動次臨界堆系統(tǒng)(accelerator driven subcritical system, ADS)的流態(tài)固體靶[4]和安全性高的第4代核反應(yīng)堆——球床式高溫氣冷反應(yīng)堆的堆芯[5]。

隨著計算機性能的提高和各種穩(wěn)定的、精確的數(shù)值算法的提出,越來越多的學(xué)者開始通過數(shù)值模擬方法研究球類自然對流換熱問題。Jia和Gogos[11-12]數(shù)值模擬10

1 模型及方法

圖1 “極限”的數(shù)學(xué)思想Fig.1 Idea of the limit in mathematics

在層流自然對流傳熱的數(shù)值研究中,浸沒在廣延和靜止空氣的物體需要足夠大的計算空間,即保證邊界層能自由地發(fā)展而不受外邊界影響。在前人研究中,習(xí)慣上將球體浸入足夠大的外環(huán)境中,此時外內(nèi)邊界的比率d∞/D約為12.5、25和40[14]。然而,巨大的計算域會嚴(yán)重降低計算效率。為此,Bejan等[20]在5～8根水平錯排管束的自然對流優(yōu)化換熱計算中提出一種更經(jīng)濟的二維計算模型。與常規(guī)模型相比,二維Bejan模型在遠端氣體出口附近設(shè)置一個入口來中和換熱過程中由于密度變化產(chǎn)生的壓力差,避免煙囪效應(yīng),同時可以使計算域的寬度減小到d∞/D=3或10從而大幅度提高計算效率。需要指出的是:煙囪效應(yīng)是指在細長區(qū)域內(nèi)流體被施加一個額外的加速度現(xiàn)象,其會影響計算的準(zhǔn)確性。該模型在對稱條件下的應(yīng)用已被證明是可行的,包括在二維條件下的半球(軸對稱結(jié)構(gòu))[19,21]和圓柱(對稱結(jié)構(gòu))[22-24]。然而它不適用于本文研究的順排布置的水平單層球結(jié)構(gòu)。因此本文基于二維Bejan模型,提出拓展后的三維模型,減少計算域和提高計算效率的同時避免煙囪效應(yīng)的產(chǎn)生,如圖2所示。此外,局部結(jié)果的展示采用圖2(b)中的球坐標(biāo)系(r,θ,φ)以便于分析后文數(shù)據(jù)。

圖2 三維拓展計算模型Fig.2 Expanded 3D computational model

這里以3×3排列水平單層球為例解釋三維模型,如圖2(a)所示。坐標(biāo)原點在中心球的球心。由于流動對稱性,僅選取切割后的1/8等腰直角三棱柱(由灰色面組成)的拓展模型進行模擬,這樣可以極大地減少計算量。對于切割后的拓展模型,每個球體和空氣的溫度分別保持恒定值Tw和T∞(Tw>T∞)。這里假設(shè)空氣(Pr=0.72)不可壓縮,其熱物性參數(shù)通過恒定的膜溫Tm(=(Tw+T∞)/2)來獲得,除了Z方向上的動量方程的密度項中,使用布辛涅斯克假設(shè)(Boussinesq approximation)說明密度變化和溫度變化之間的關(guān)系:

(ρ∞-ρ)=ρβ(T-T∞).

(1)

為減少誤差,本文在數(shù)值模擬中使用大約1 K的溫度差來保證布辛涅斯克假設(shè)的有效性。由于溫差極小,空氣可近似為理想氣體,所以其熱膨脹系數(shù)β可簡化為β=1/Tm。另外,由于本文中自然對流是低速層流且采用空氣不可壓縮假設(shè),所以能量方程中忽略黏性耗散和輻射。因此,三維、穩(wěn)態(tài)、常物性等條件下的無量綱控制方程可以簡化成如下:

連續(xù)性方程

(2)

X方向動量方程

(3)

Y方向動量方程

(4)

Z方向動量方程

(5)

能量方程

(6)

其中拉普拉斯算子Δ定義為Δ=?2/?X*2+?2/?Y*2+?2/?Z*2。

上述控制方程中無量綱變量的定義,圖2(a)中計算域的邊界條件的設(shè)定和文獻[14]保持一致。

在開始本文研究內(nèi)容之前,需要進行區(qū)域無關(guān)性驗證和網(wǎng)格無關(guān)性驗證,消除計算域大小和網(wǎng)格數(shù)量對于結(jié)果的影響。入口段和出口段的垂直長度參考文獻[14]使用恒定的3D,如圖2(a)所示。對于計算域?qū)挾萪∞,(2N+1)×(2N+1)排列模型選擇d∞=0.52N+1D+N×0.01D+2D。舉例來說,3×3排列模型選擇d∞=1.5D+0.01D+2D,1.5D+0.01D表示中心球球心到最外圍球的距離,其中0.01D表示3×3排列模型只存在一個間隙,2D表示最外圍球到模型外邊界的距離。間隙的存在是由于模型無法做到任意兩球緊密接觸,因此在網(wǎng)格劃分時采用極小間距(0.01D)來代替緊密接觸,同時避免產(chǎn)生質(zhì)量差的網(wǎng)格。這里要說明的是,本文模擬的多組排列模型中最外圍球到模型外邊界的距離均使用2D,這是由于固定的Gr數(shù)使得溫度和速度邊界層不會因為排列數(shù)的增多而有劇烈變化。

表1 區(qū)域無關(guān)性驗證(3×3排列)Table 1 Domain independence tests for 3×3 array

綜合上述區(qū)域無關(guān)性驗證和網(wǎng)格無關(guān)性驗證,順排布置的水平單層球自然對流換熱模型使用以下參數(shù):縱向長度Z=±3D, 橫向?qū)挾萪∞=0.52N+1D+N×0.01D+2D,另外網(wǎng)格始終使用G1。

表2 網(wǎng)格無關(guān)性驗證(3×3排列)Table 2 Grid independence tests for 3×3 array

2 結(jié)果及討論

2.1 模型方法驗證

2.2 溫度分布和流線分布

為了定性說明排列數(shù)對于水平單層球自然對流換熱的影響,本節(jié)展示3×3和15×15排列的水平單層球自然對流的溫度分布和流線分布。由于球面是定溫的,所以持續(xù)與球周圍的冷空氣發(fā)生熱交換,而這些冷空氣吸熱后產(chǎn)生密度差,進而在浮升力作用下圍繞球體不斷向上運動,在這個過程中還裹挾多球區(qū)域外圍的大量空氣來形成穩(wěn)定運動的卷流。隨著排列數(shù)的增加,最外圍球附近的溫度分布的整體情況沒有改變,但是對于其他位置的大多數(shù)球,15×15排列的溫度分布與3×3排列不相同,特別是對于層球間隙附近的溫度場。從圖4(a1)和圖4(b1)的溫度分布可以看出,隨著排列數(shù)增加,層球間隙附近區(qū)域吸入了更多的冷空氣。這是因為隨著排列數(shù)增加,層球上方的熱空氣被進一步加熱,在強浮力的驅(qū)動下向上運動,由此導(dǎo)致該區(qū)域負壓增大,具體可見2.3節(jié)中的CP分布。雖然這會導(dǎo)致外圍空氣快速補充進來,但中間區(qū)域因為無法直接和這些補充的空氣接觸,從而會通過單層球的間隙從下方吸入大量冷空氣,這使得球和周圍空氣換熱增強,具體可見2.4節(jié)的Nu數(shù)分布。該負壓區(qū)存在范圍較廣,基本涵蓋大多數(shù)的球,且越靠近單層球排列的中心,負壓越大,通過間隙吸入的冷空氣越多。此外,與單球的溫度和流線分布不同,由于外圍空氣向中間位置靠攏,單層球上方的溫度邊界層會受到擠壓,流線產(chǎn)生明顯的彎曲,這從圖4(a2)和圖4(b2)之間的溫度和流線分布的對比可以看出。

表3 三維單球的驗證結(jié)果Table 3 Verification results of three-dimensional single sphere

圖4 不同排列數(shù)下水平單層球的溫度分布和流線分布Fig.4 Isothermal and streamline contours in the vicinity of the horizontal multi-spheres

2.3 局部阻力系數(shù)分布

圖5中展示方位角為φ=π/4位置處中心球的自然對流換熱的CF分布(圖5(a))和CP分布(圖5(b))。總的來說,隨著排列數(shù)增加,CF逐漸增大,CP的絕對值也逐漸增大。這是因為層球上方的卷流隨著排列數(shù)的增加受到的加熱效應(yīng)更強,導(dǎo)致層球附近區(qū)域負壓增大,即壓力梯度增大,由此增強了流體的流動能力。而且層球下游負壓受排列數(shù)的影響更明顯。對于某一排列數(shù)而言,隨著θ的減小(從上游導(dǎo)向下游),中心球壁面附近的流體首先在順壓力梯度下加速,速度梯度也相應(yīng)增大。當(dāng)θ減小到π/2時,雖然流體在順壓力梯度下依然會加速,但因為受到下游卷流的影響,速度邊界層逐漸增厚,速度梯度開始隨著θ的減小而減小。當(dāng)θ減小到3π/7時,流體在逆壓力梯度下開始減速,在下游卷流的共同作用下,速度梯度減小幅度加大,并在θ=0～π/4區(qū)間達到最小且基本保持為零。

圖5 φ=π/4下中心球的局部阻力系數(shù)分布Fig.5 Variation of local drag coefficient along the surface of central sphere from 3×3 to 15×15 for φ=π/4

2.4 局部努塞爾數(shù)分布

圖6展示方位角為φ=π/4位置處不同排列數(shù)(從3×3至15×15)的Nu數(shù)的影響?？偟膩砜?Nu數(shù)隨著排列數(shù)增加而增大。這是由于對流效應(yīng)隨著排列數(shù)增加而增強,從而加大了球體和流體之間的換熱,這種區(qū)別在上游區(qū)域更加明顯。這是因為上游區(qū)域的溫度邊界層相對于下游區(qū)域受排列數(shù)的影響更大。對比13×13和15×15這2種排列,可以發(fā)現(xiàn)當(dāng)排列數(shù)增加到一定限度后,Nu數(shù)分布基本不變,這表明本文引入的“極限”的數(shù)學(xué)思想的合理性。通常來說,Nu數(shù)隨著θ的減小而單調(diào)遞減,即球體底端的Nu數(shù)最大而頂端的Nu數(shù)最小。在θ=5π/8～π區(qū)域,Nu數(shù)變化不大。這是因為雖然流速的增加會增強換熱,但中心球附近球體對流體的預(yù)熱減小了溫差。兩種因素共同作用使得中心球Nu數(shù)變化不大。當(dāng)θ從5π/8開始減小,附近球體對流體的預(yù)熱作用增強,中心球Nu數(shù)逐漸減小。隨著θ進一步減小,流體速度的減小使得Nu數(shù)減小幅度變大。在θ=0～π/4區(qū)間,換熱方式主要是熱傳導(dǎo),因為上方卷流及周圍球體的作用,中心球和流體之間的溫差較小,換熱能力較弱。

圖6 φ=π/4下中心球的局部努賽爾數(shù)分布Fig.6 Variation of local Nusselt number along the surface of two horizontal spheres from 3×3 to 15×15 for φ=π/4

2.5 平均努塞爾數(shù)

圖7 中心球的平均努賽爾數(shù)分布Fig.7 Average Nusselt number distribution of the central sphere

(7)

3 結(jié)論

1)隨著排列數(shù)增加,卷流作用的增強使得單層球上方區(qū)域負壓增大。這不僅導(dǎo)致層球上方溫度邊界層受到擠壓,而且使得通過球體間隙從下方吸入的流體速度升高。越靠近單層球排列的中心,負壓越大,吸入的流體速度越大。

2)總的來看,Nu數(shù)與排列數(shù)呈正相關(guān)。由于上游區(qū)域溫度邊界層受排列數(shù)影響更大,排列數(shù)對該區(qū)域Nu數(shù)的影響更明顯。通常來說,Nu數(shù)隨著θ的減小而單調(diào)遞減。對于θ=5π/8～π區(qū)域,Nu數(shù)隨θ減小幅度較小。隨著流體的減速及周圍球體對中心球預(yù)熱的增強,Nu數(shù)隨θ減小幅度加大。在θ=0～π/4區(qū)間,Nu數(shù)達到最小。