高普朗特數(shù)流體圓柱坐標(biāo)下自然對流數(shù)值模擬研究1

邵倩倩，王敏，許諾，梁飛華

(1.廣東石油化工學(xué)院 石油工程學(xué)院，廣東 茂名 525000；2.北京計算科學(xué)研究中心，北京 100193)

在石油儲運行業(yè)中，儲罐內(nèi)油品的傳熱問題是圓柱坐標(biāo)下自然對流在工程上的典型應(yīng)用之一。相較于其他流體，石油產(chǎn)品的特征是普朗特數(shù)(Pr)較高。因此，對高Pr流體在圓柱坐標(biāo)下自然對流規(guī)律進(jìn)行研究是進(jìn)一步深入解決油品儲存過程中所涉及的各種流動與傳熱問題的基礎(chǔ)。如，以圓柱坐標(biāo)下自然對流為基礎(chǔ)，進(jìn)一步研究儲罐和管道內(nèi)油品的溫度分布，并根據(jù)油品閃點、沸點、凝點等指導(dǎo)油品的節(jié)能生產(chǎn)與儲運安全[1-3]；以不同條件(如晝夜變化、天氣變化等)下儲罐內(nèi)油品溫度分布為基礎(chǔ)，可以進(jìn)一步研究油品的蒸發(fā)損耗，同時為減小蒸發(fā)損耗提供寶貴決策[4-6]。

自然對流問題的溫度場與速度場受物理區(qū)域的幾何特征、Pr、瑞利數(shù)(Ra)以及邊界條件等共同影響。兩同心圓筒壁間區(qū)域是常見的用于研究圓柱坐標(biāo)下自然對流問題的幾何形狀。Powe等[7]采用數(shù)值計算方法研究了不同Ra下水平圓柱環(huán)形空腔內(nèi)的自然對流，其中Pr選為0.7。Kuehn等[8]采用實驗和數(shù)值模擬對水平圓柱環(huán)形空間內(nèi)的自然對流進(jìn)行研究，環(huán)形寬度與內(nèi)徑的比值為0.8，Pr選為0.7、5.45，Ra范圍為2.11×104～9.76×105，其實驗結(jié)果與數(shù)值計算結(jié)果吻合良好。Basant等[9]對不對稱加熱或冷卻的兩個立式圓筒之間的自然對流進(jìn)行研究，分析了浮力分布參數(shù)的影響，獲得了瞬態(tài)自然對流的半解析解。Kumar等[10]關(guān)注立式圓柱環(huán)形空間，對Pr為0.7、Ra為10～106、不同高徑比幾何特征下的自然對流傳熱與流動進(jìn)行研究。調(diào)研發(fā)現(xiàn)，以往對圓柱坐標(biāo)下自然對流問題的研究多針對低Pr流體展開。本研究以Pr=1000為例，采用數(shù)值計算方法模擬兩同心圓筒壁間的自然對流。內(nèi)圓筒高溫、外圓筒低溫、頂部和底部絕熱，Ra取為1.0×103、1.0×104、1.0×105及1.0×106。采用加密網(wǎng)格和多重網(wǎng)格算法，得到計算區(qū)域的溫度場、流場以及流函數(shù)分布等，基于數(shù)值計算結(jié)果研究高Pr流體圓柱坐標(biāo)自然對流規(guī)律，分析Ra對溫度場與流場的影響。

1 數(shù)值計算方法

1.1 物理問題及簡化

本研究物理區(qū)域(如圖1a所示)為內(nèi)圓筒與外圓筒之間的區(qū)域?？紤]三維圓柱坐標(biāo)下，物理量沿周向方向(圖1a中θ所在的方向)的變化可近似忽略，本文將三維圓柱坐標(biāo)簡化為二維圓柱坐標(biāo)(見圖1b)。其中，x方向為高度方向，朝上為正方向；r坐標(biāo)指定朝外圓筒為正方向，坐標(biāo)原點位于圓柱底的圓心。對于本研究計算區(qū)域，徑向上長度(即內(nèi)、外圓筒壁徑向上的間距)為l，高度方向長度為h。

1.2 控制方程

圖1 三維圓柱坐標(biāo)物理問題及其二維簡化

圓柱坐標(biāo)下自然對流的控制方程包括連續(xù)性方程、x方向動量方程、r方向動量方程以及能量方程，其分別為

(1)

(2)

(3)

(4)

式中：r、x分別為r方向、x方向坐標(biāo)值，m；ρ為密度，kg/m3；cp為比熱容，J/(kg·K)；p為壓力，Pa；g為重力加速度，取9.81 m/s2；μ為動力黏度，Pa·s；β為熱膨脹系數(shù)，1/K；λ為導(dǎo)熱系數(shù)，W/(m·K)；T為溫度，K；u和v分別為x方向、r方向的速度分量，m/s。

引入無量綱量R=r/l，X=x/l，U=ul/υ，V=vl/υ，P=pl2/ρυ2，Θ=(T-Tc)/(Th-Tc)，Pr=ρcpυ/λ，Ra=Gr·Pr=ρgβcp(Th-Tc)l3/λυ(Gr為格拉曉夫數(shù)，定義為Gr=gβ(Th-Tc)l3/υ2；l為計算區(qū)域的寬度，m；υ為運動黏度，m2/s；Th、Tc分別為高溫邊界和低溫邊界的溫度，K)。

對式(1)～(4)進(jìn)行無量綱化，可得：

(5)

(6)

(7)

(8)

1.3 數(shù)值計算方法及驗證

為了獲得更高的計算精度，采用渦量-流函數(shù)法來計算流函數(shù)，當(dāng)流函數(shù)的殘差小于1.0×10-10時，計算結(jié)果認(rèn)為是收斂的；計算區(qū)域采用1026×1026的非均分網(wǎng)格進(jìn)行離散，網(wǎng)格足夠密，所得結(jié)果具有高精度；采用多重網(wǎng)格算法來提高數(shù)值計算速度[11, 12]。

采用本文計算方法對經(jīng)典文獻(xiàn)[10]中問題進(jìn)行計算，其中Pr=0.7，Ra=2×105。所得結(jié)果與文獻(xiàn)[10]中結(jié)果進(jìn)行比較，如圖2所示。其中R1為無量綱內(nèi)半徑，R2為無量綱外半徑。

圖2 文獻(xiàn)[10]與本研究數(shù)值計算方法的準(zhǔn)確性

由圖2可知，采用本文計算方法所得結(jié)果與文獻(xiàn)[10]中結(jié)果吻合較好，因此，本研究的數(shù)值計算方法被認(rèn)為是精確有效的。

2 計算結(jié)果與分析

2.1 計算條件

本研究中考慮的溫度邊界條件如圖3所示，其中內(nèi)邊界為高溫(Th)邊界，外邊界為低溫(Tc)邊界，上邊界和下邊界均為絕熱邊界 ?T/?x=0。各邊界速度設(shè)置為無滑移條件，即各邊界上u=0，v=0(即U=0，V=0)。Ra設(shè)為1.0×103、1.0×104、1.0×105及1.0×106，Pr設(shè)為1000，代表油品這一類高Pr數(shù)流體。無量綱內(nèi)半徑為R1=0.1；無量綱外半徑為R2=1.1；無量綱高度為H=1。

圖3 二維圓柱坐標(biāo)下自然對流溫度邊界條件

2.2 溫度場與速度場分布

圖4為不同Ra下計算所得的X-R坐標(biāo)下的溫度場與速度場，給出的參照速度矢量單位為m/s。

圖4 不同Ra下的溫度場與速度場

由圖4可知，(1)高溫流體對流到低溫區(qū)域，低溫流體對流到高溫區(qū)域，流體在區(qū)域內(nèi)呈現(xiàn)順時針旋轉(zhuǎn)。(2)中心區(qū)域流速較小，這是由于中心區(qū)域流體溫差小，密度差也小，流體對流較慢。邊界處速度較大則是由于較大的溫差產(chǎn)生的大的密度差，使流體運動劇烈。(3)Ra增大，自然對流增強(qiáng)，流體流動更加劇烈，速度整體增大。(4)當(dāng)Ra增大時，高溫區(qū)域從內(nèi)圓筒向外圓筒及上部轉(zhuǎn)移，一方面，內(nèi)外圓筒溫差的作用，熱量不斷地由內(nèi)圓筒傳到外圓筒；另一方面，對流使得高溫流體向上方移動。

2.3 流函數(shù)場

圖5為不同Ra下計算所得的流函數(shù)場。由圖5可知，(1)不同Ra下流函數(shù)場均清晰地顯示一個主渦。(2)Ra對流函數(shù)的分布和數(shù)值均有很大的影響，當(dāng)Ra較小時，流函數(shù)的等值線形狀較為規(guī)則，不同等值線近似同心均勻分布；當(dāng)Ra增大時，各等值線形狀差別加大，且形狀變得不規(guī)則，流動變得更加復(fù)雜，表現(xiàn)出越來越明顯的湍流特征。(3)當(dāng)Ra增大時，主渦的位置向低溫邊界及上邊界轉(zhuǎn)移，即向圖中的右上角移動，同時，最大流函數(shù)數(shù)值明顯增加。

圖5 不同Ra下的流函數(shù)分布(圖中所標(biāo)數(shù)字代表流函數(shù)數(shù)值)

為了更準(zhǔn)確地展示這一規(guī)律，表1給出了不同Ra下主渦中心的坐標(biāo)以及對應(yīng)的最大流函數(shù)值。從表1可以直觀地看到，當(dāng)Ra增大時，主渦中心的R坐標(biāo)以及X坐標(biāo)均增大，最大流函數(shù)值更是明顯地增大。

表1 不同Ra下計算所得主渦中心坐標(biāo)以及最大流函數(shù)值Ra=1.0×103Ra=1.0×104Ra=1.0×105Ra=1.0×106主渦中心R坐標(biāo)值0.64250.66880.81040.9413主渦中心X坐標(biāo)值0.51420.54540.59910.7446最大流函數(shù)值0.00050.00220.00490.0080

2.4 Ra對溫度場與流場的影響

圖6展示了Ra對計算區(qū)域兩條中心線(X=1/2，R=(R1+R2)/2)上溫度分布的影響，圖7展示相應(yīng)的Ra對速度分布的影響(圖7僅給出垂直于各中心線的速度分量的情況)。

圖6 不同Ra下計算區(qū)域兩條中心線上溫度分布 圖7 不同Ra下計算區(qū)域兩條中心線上速度分布

由圖6可知，對于Θ～R分布曲線，當(dāng)Ra較小時，從內(nèi)圓筒到外圓筒沿著X=1/2中心線，溫度逐漸降低，其分布接近對數(shù)分布，這是由于低Ra下，自然對流較弱，熱傳導(dǎo)占主要地位，傳熱過程類似于純導(dǎo)熱過程。隨著Ra增大，X=1/2中心線上Θ～R分布越來越偏離對數(shù)分布。當(dāng)Ra較大時，在內(nèi)圓筒處溫度急劇下降，遠(yuǎn)離內(nèi)圓筒處，溫度分布近似均勻。對于Θ～X曲線，高溫存在于X較大的區(qū)域，即上部區(qū)域，這是高溫流體向上運動導(dǎo)致的。而Ra增大使得高溫流體向上聚集的現(xiàn)象更為明顯，表現(xiàn)為沿R=(R1+R2)/2這條中心線，上邊界附近溫度梯度隨著Ra的增加而增大。

由圖7可知，由于邊界無滑移速度條件以及自然對流的作用，在各邊界附近，速度均表現(xiàn)為先增大后減小的分布特征，速度梯度隨Ra增大而變大；Ra對速度大小影響亦十分明顯，當(dāng)Ra增大時，流速顯著增大。例如，Ra=1.0×106時，中心線X=1/2上最大U速度約為Ra=1.0×103時的85倍。

為了直觀顯示自然對流的強(qiáng)弱，給出努塞爾數(shù)(Nu)在內(nèi)外邊界上的分布情況，見圖8。其中Nu1為高溫邊界上Nu，Nu2為低溫邊界上Nu。圖8a、圖8b分別為高溫邊界、低溫邊界上Nu沿X方向的分布情況。由圖8可知，在高溫內(nèi)邊界及低溫外邊界上，Nu均隨著Ra增加而顯著增加，即對流換熱增強(qiáng)。

圖8 不同Ra下高溫邊界及低溫邊界上Nu分布

3 結(jié)語

使用數(shù)值計算的方法對高Pr流體在圓柱坐標(biāo)下的自然對流問題進(jìn)行研究，基于精密網(wǎng)格劃分，采用有限容積法并利用多重網(wǎng)格算法來提高計算效率，借助經(jīng)典文獻(xiàn)中的數(shù)據(jù)驗證了本文數(shù)值計算方法的精確可靠。研究得到了不同Ra下自然對流的溫度場、速度矢量場以及流函數(shù)場的分布情況，并對中心線上的溫度分布、速度分布以及高溫邊界和低溫邊界上的Nu分布進(jìn)行了分析，主要結(jié)論如下：(1)Ra對溫度場和流場均有較大的影響，反映自然對流的強(qiáng)弱；(2)當(dāng)Ra增大時，高溫流體對流到低溫區(qū)域、低溫流體對流到高溫區(qū)域的現(xiàn)象更加劇烈，使上下區(qū)域溫度分層現(xiàn)象更加明顯，且對流換熱增強(qiáng)，表現(xiàn)為明顯增大的Nu；(3)Ra增大，速度整體上增大，且最大流函數(shù)值增加，主渦位置向低溫邊界及上邊界移動。