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

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

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

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

1 數(shù)值計(jì)算方法

1.1 物理問(wèn)題及簡(jiǎn)化

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

1.2 控制方程

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

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

(1)

(2)

(3)

(4)

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

引入無(wú)量綱量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為計(jì)算區(qū)域的寬度，m；υ為運(yùn)動(dòng)黏度，m2/s；Th、Tc分別為高溫邊界和低溫邊界的溫度，K)。

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

(5)

(6)

(7)

(8)

1.3 數(shù)值計(jì)算方法及驗(yàn)證

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

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

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

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

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

2.1 計(jì)算條件

本研究中考慮的溫度邊界條件如圖3所示，其中內(nèi)邊界為高溫(Th)邊界，外邊界為低溫(Tc)邊界，上邊界和下邊界均為絕熱邊界 ?T/?x=0。各邊界速度設(shè)置為無(wú)滑移條件，即各邊界上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，代表油品這一類(lèi)高Pr數(shù)流體。無(wú)量綱內(nèi)半徑為R1=0.1；無(wú)量綱外半徑為R2=1.1；無(wú)量綱高度為H=1。

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

2.2 溫度場(chǎng)與速度場(chǎng)分布

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

圖4 不同Ra下的溫度場(chǎng)與速度場(chǎng)

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

2.3 流函數(shù)場(chǎng)

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

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

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

表1 不同Ra下計(jì)算所得主渦中心坐標(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對(duì)溫度場(chǎng)與流場(chǎng)的影響

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

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

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

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

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

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

3 結(jié)語(yǔ)

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