謝雄州, 葉道銘, 黃俊輝, 王建勇

(中山大學·深圳 航空航天學院,深圳 518107)

超臨界二氧化碳(supercritical carbon dioxide, sCO2)是一種在動力循環(huán)和能量轉(zhuǎn)換中能顯著提高熱效率的工質(zhì)[1],在動力工程領域有著廣闊的應用前景[2-3]。不同于一般的常物性流體,sCO2的熱物理性質(zhì)在擬臨界點附近會發(fā)生劇烈變化。由于物理性質(zhì)的劇烈變化而誘發(fā)的浮升力或熱加速效應導致流場和溫度場之間強烈耦合,從而影響流場和湍流結(jié)構(gòu)[4-5],表現(xiàn)為傳熱強化(HTE)和傳熱惡化(HTD)等異常傳熱現(xiàn)象。傳熱惡化的發(fā)生會導致傳熱系數(shù)降低,局部壁面溫度急劇升高,傳熱效率降低,并對超臨界換熱裝置的安全、穩(wěn)定運行帶來嚴峻挑戰(zhàn)。因此,深刻且準確了解sCO2湍流流動傳熱行為對于相關(guān)超臨界流體換熱部件的設計及研發(fā)至關(guān)重要。

為了可靠地預測超臨界流體的流動傳熱特性,學者們嘗試了各種不同的數(shù)值方法和模型。Reynolds Average Navier-Stokes(RANS)湍流模型在計算精度與成本之間達到了較好的平衡,因此在工程中得到廣泛應用,但傳熱惡化時其對超臨界換熱的預測性能通常較差。對于RANS湍流模型,湍流普朗特數(shù)Prt是一個重要參數(shù),且通常設定為定常值(Prt=0.85)。此舉對常物性流體而言,其合理性得到普遍接受,但對于物性變化劇烈的超臨界流體,其適用性遭到廣泛質(zhì)疑。因此,一些學者嘗試提出湍流普朗特數(shù)修正模型以改善RANS湍流模型對超臨界湍流傳熱的預測能力。Bae[6]基于混合長度理論提出了變Prt模型,并將其應用于多種超臨界工質(zhì)在豎直管內(nèi)的湍流傳熱計算中,結(jié)果顯示修正模型的計算結(jié)果與實驗數(shù)據(jù)吻合良好。Tang等[7]在Kays[8]的總結(jié)推導上,基于SSTk-ω模型提出了分段函數(shù)式的變Prt模型;Du等[9]又在此模型基礎上進一步引入壓力修正因子和管徑修正因子,形成了適用性更廣的變模型。此外,還有眾多研究者(Tian等[10]、Kong等[11])提出了一系列變湍流普朗特數(shù)模型。然而,上述變Prt模型的提出背景及修正應用均是面向豎直超臨界管流的。水平布局同樣在超臨界換熱設備中被廣泛采用,但不同于豎直管道換熱所呈現(xiàn)的周向均勻性,水平換熱過程中密度不同的sCO2受浮升力的驅(qū)使會產(chǎn)生顯著的流體分層現(xiàn)象,繼而導致巨大的上、下壁溫差。而RANS湍流模型應用于水平管道超臨界傳熱計算時同樣存在較大偏差,此情形下湍流普朗特數(shù)的適用性及其對計算可信度的影響有待進一步剖析,同時在豎直管流動預測上卓有成效且普適性較好的變Prt修正模型應用于水平超臨界流動傳熱的預測亦值得嘗試。

鑒于此,筆者分析了定常湍流普朗特數(shù)(Prt=0.85)在強、弱2種浮升力效應下對水平sCO2湍流傳熱的預測性能,而后引入現(xiàn)有典型的Prt修正模型進行數(shù)值模擬,并將模擬結(jié)果與實驗數(shù)據(jù)進行詳細對比,以探究傳熱預測的可靠性,并系統(tǒng)分析Prt及其相關(guān)修正模型對水平sCO2流動傳熱計算的影響。研究成果對水平超臨界湍流傳熱RANS湍流模型模擬及其深層修正理論具有重要意義。

1 物理模型與數(shù)值方法

1.1 幾何模型

本文開展數(shù)值模擬時所采用的物理模型如圖1所示,圖中:g為重力加速度;din為圓管內(nèi)徑;φ為圓管周向角;L1、L2和L3分別表示圓管發(fā)展段、加熱段和緩沖段的長度。該模型的幾何尺寸與Theologou等[12]和Adebiyi等[13]的實驗圓管參數(shù)相同,其針對水平sCO2的湍流傳熱進行了寬工況范圍下系統(tǒng)而細致的測量,具體參數(shù)見表1。加熱段前、后分別設置了絕熱段以降低進、出口效應。圓管周向角度φ為0°、90°(270°)和180°的位置分別對應圓管的頂母線、中母線和底母線。

表1 計算模型幾何參數(shù)

圖1 物理模型

1.2 控制方程

采用ANSYS Fluent 19.2軟件求解器開展模擬計算,笛卡爾坐標系下連續(xù)性方程、動量方程、能量方程的穩(wěn)態(tài)形式分別如式(1)～式(3)所示。

連續(xù)性方程:

(1)

動量方程:

(2)

能量方程:

(3)

SSTk-ω模型由于近壁區(qū)域作了低雷諾數(shù)處理,可以滿足邊界層高精度的模擬需求,其在超臨界流體傳熱計算中表現(xiàn)出色,受到眾多研究者的青睞[14-15],因此也被選用至本文所探究的sCO2水平流動傳熱模擬中。SSTk-ω模型中湍流控制方程如下:

湍動能k輸運方程:

(4)

比耗散率ω輸運方程:

(5)

式中:Γk為湍動能k的有效擴散系數(shù);Γω為比耗散率ω的有效擴散系數(shù);Gk為湍動能生成項;Gω為耗散率生成項;Yk為湍動能耗散項;Yω為比耗散率ω的耗散項;Dω為交叉擴散項;Sk為湍動能k的自定義源項;Sω為比耗散率ω的自定義源項。

更多模型細節(jié)可參考文獻[16]。

1.3 典型湍流普朗特數(shù)修正模型

現(xiàn)有Prt修正模型大多面向豎直管內(nèi)超臨界流動傳熱計算,且主要分為2類:第一類基于豎直管內(nèi)超臨界湍流傳熱狀態(tài),依靠相關(guān)理論推導及無量綱分析而得出,如Bae[6]基于混合長度理論提出的變Prt修正模型;第二類主要基于管內(nèi)超臨界湍流傳熱特征,對流動進行分層(尤其是邊界層內(nèi))繼而相應合理賦值湍流普朗特數(shù)。第一類湍流普朗特數(shù)修正模型的提出背景及依據(jù)與水平流動情形明顯不同,故不能直接應用。因此,本文引入并測試第二類Prt修正模型針對水平sCO2湍流傳熱的計算預測。

此處選取普適性較好的Tang等[7]和Du等[9]提出的變Prt修正模型。其中,Tang等[7]提出的模型描述為:

(6)

式中:A為可調(diào)節(jié)的常數(shù),Tang等[7]推薦A取值為15。

而Du等[9]提出的模型為:

(7)

式中:fd、α1、α2均為管徑修正因子,fd=6.2(din-4)0.24,α1=0.24din-0.7,α2=0.05din-0.07;fp為壓力修正因子,fp=[1+0.019(p/pcr)29],pcr為臨界壓力。

為方便敘述,Du等[9]、Tang等[7]提出的模型在后文分別稱為Du模型和Tang模型。

1.4 網(wǎng)格劃分與無關(guān)性驗證

流體計算域網(wǎng)格由ANSYS ICEM軟件生成,采用結(jié)構(gòu)化六面體網(wǎng)格,圓管橫截面采用O形網(wǎng)格劃分。由于黏性底層和緩沖層的網(wǎng)格質(zhì)量對SSTk-ω低雷諾數(shù)湍流模型求解近壁處流動傳熱細節(jié)至關(guān)重要,因此在近壁面處對網(wǎng)格進行局部加密,劃分的網(wǎng)格橫截面如圖2所示。

圖2 計算域橫截面網(wǎng)格

為檢驗計算網(wǎng)格的獨立性,分別對內(nèi)徑為4 mm、8 mm和22.14 mm的水平圓管進行網(wǎng)格無關(guān)性驗證。以內(nèi)徑為4 mm的水平圓管為例,在保持第一層網(wǎng)格無量綱高度y+<1的情況下,對流體域近壁面處網(wǎng)格進行不同程度的加密處理,分別得到84.4萬、115.8萬、280.6萬和445.3萬等4種不同密度的計算網(wǎng)格。這些網(wǎng)格在質(zhì)量流量G=800 kg/(m2·s)、均勻加熱熱流密度q=30 kW/m2、入口溫度Tin=30 K的工況下計算得到的加熱壁面平均溫度Tw隨軸線位置的變化分布如圖3所示。由圖3可以看出,當網(wǎng)格數(shù)大于280.6萬時,加熱壁面平均溫度隨網(wǎng)格數(shù)增加已無明顯變化。為節(jié)省計算成本,4 mm水平圓管的計算網(wǎng)格數(shù)選定為280.6萬。同理,分別對內(nèi)徑為8 mm和22.14 mm水平圓管的計算模型進行網(wǎng)格獨立性分析,選用網(wǎng)格數(shù)分別為360.4萬和437.6萬。

圖3 網(wǎng)格獨立性驗證

1.5 邊界條件與計算方法

為與相應實驗測量條件保持一致,表1中“TS”“TL”系列工況入口設置為質(zhì)量流速入口,“A”系列工況入口設置為質(zhì)量流量入口,出口均設置為壓力出口,入口發(fā)展段和出口緩沖段的壁面設置為絕熱。加熱管壁設置為均勻熱流密度且為無滑移邊界。邊界參數(shù)與實驗工況相同,具體見表2。本文所選擇的驗證工況均處于臨界點附近,CO2物性變化明顯,典型熱物性參數(shù)具體范圍為:1.14 kJ/(kg·K)

表2 驗證工況

控制方程采用有限體積法進行離散,壓力-速度耦合運用Simplec算法求解,動量方程和能量方程采用Quick格式,動量方程中的壓力項應用Body Force Weighted離散,除密度采用二階迎風格式外,湍動能和耗散項均采用一階迎風格式以獲取更優(yōu)的收斂性能。計算過程中,sCO2真實氣體效應物性參數(shù)調(diào)用自基于求解器內(nèi)嵌的NIST物性庫所生成的查找表,變Prt修正模型通過UDF(User-Defined Functions)載入Fluent軟件求解器。

2 結(jié)果與討論

2.1 恒定湍流普朗特數(shù)(0.85)模型對壁溫的預測

為了驗證SSTk-ω湍流模型采用定常Prt=0.85模型開展預測的準確性和可靠性,基于Theologou等[12]和Adebiyi等[13]的實驗結(jié)果,開展了強、弱浮升力效應下的仿真驗證。

2.1.1 弱浮升力效應

圖4給出了弱浮升力效應下(工況 TS1、TL1、A1和A2)定常湍流普朗特數(shù)(Prt=0.85)模型對圓管內(nèi)壁溫Tw,i的預測值與實驗值隨主流焓值ib增加的對比結(jié)果。圖4(a)的結(jié)果表明該工況下的浮升力還不足以使上下管壁產(chǎn)生明顯的溫度差,3條母線處的模擬值與實驗值在趨勢和數(shù)值上均有著較高的吻合度。圖4(b)～圖4(d)中的工況下受熱流密度q和管徑d增大的影響,sCO2頂部和底部溫差略微增大,但采用Prt=0.85開展仿真的計算結(jié)果仍能較好地預測壁溫的變化趨勢,且頂母線處壁溫預測偏差小于5 K。以上驗證結(jié)果表明,浮升力較小時,采用Prt=0.85時模型具有較好的預測可信度。

(a) 工況TS1

2.1.2 強浮升力效應

浮升力效應隨管徑d的增大和熱流密度q與質(zhì)量流速G比值(q/G)的增加愈發(fā)顯著,直觀反映于上下壁溫差[19-21]。圖5給出了強浮升力效應下(工況 TS2、TS3、TL2、TL3、A3和A4)定常湍流普朗特數(shù)(Prt=0.85)模型對壁溫預測值與實驗值的對比結(jié)果。圖5(a)、圖5(c)、圖5(e)的結(jié)果顯示,強浮升力效應下,當q/G相對較小時,頂母線處壁溫預測值較實驗值存在明顯偏差,存在高估或者低估上壁溫的現(xiàn)象。而隨著q/G增大,浮升力進一步加強,如圖5(b)、圖5(d)和圖5(f)所示,兩者偏差進一步加大,上述現(xiàn)象在頂母線處壁溫預測時尤為突出。由此可見,強浮升力效應下仍采用Prt=0.85開展sCO2湍流換熱預測將變得“力不從心”。

(a) 工況TS2

2.2 變湍流普朗特數(shù)修正模型對壁溫的預測

由于強浮升力效應下Prt=0.85模型的預測效果不佳,本節(jié)嘗試引入典型Prt修正模型,結(jié)果如圖6所示。由圖6(a)和圖6(b)可知,傳熱惡化嚴重時,所嵌入的變Prt修正模型并未能改善預測性能。其中,針對圖6(b)中內(nèi)徑為8 mm的水平圓管,Du模型計算值在趨勢上與Prt=0.85時的近乎一致,上壁溫預測值略偏高(受限于構(gòu)造形式,Du模型無法應用于4 mm管徑的修正)。然而,Tang模型在溫度預測上的表現(xiàn)比Prt=0.85時更差,尤其是在入口處出現(xiàn)非常尖銳的上壁溫峰值,該峰值隨著A增大而降低。這是因為邊界層緩沖層中Prt=0.85+Pr/A(Pr恒為正值),A越大,緩沖層的Prt越小?？梢栽O想,當A趨于無窮大時,Tang模型將會無限逼近Prt=0.85時的壁溫預測值。

(a) 工況TS3

圖6(c)和圖6(d)所反映的傳熱惡化程度較圖6(a)和圖6(b)的相對較輕。其中,如圖6(c)所示,載入Du和Tang等提出的Prt模型修正后,在加熱段下游,上壁溫預測值更加接近于實驗值,但變化趨勢并無明顯改善,同時下壁溫預測值與實驗值的差異略微擴大。當q/G進一步增大使得浮升力效應加強時,Du模型和Tang模型的溫度預測值均更加偏離實驗值。綜合圖6可知,利用變Prt模型開展計算所得的壁溫預測值均大于Prt=0.85時的預測值。這是由于變Prt模型在黏性底層、過渡層和湍流核心區(qū)的湍流普朗特數(shù)均不小于0.85,而根據(jù)能量方程(3),湍流普朗特數(shù)所在項為湍流擴散項,湍流普朗特數(shù)的增大使得湍流摻混效應減弱,傳熱性能下降,繼而壁溫升高,這與文獻[7]和文獻[22]的結(jié)論一致。

2.3 湍流普朗特數(shù)模型對水平sCO2流動傳熱的影響

上述結(jié)果顯示在強浮升力效應工況下,現(xiàn)有典型的分段賦值變Prt修正模型對水平sCO2流動傳熱的預測改善效果不佳。本節(jié)基于工況 TL3這一典型強浮升力效應工況進行計算,從主流(“宏觀”)、邊界層(“微觀”)2個層面深層次探究Prt模型對水平sCO2湍流傳熱模擬的影響。為了更加全面地揭示Prt的影響規(guī)律及作用機制,此處直接改變定常Prt值進行計算分析。

2.3.1 湍流普朗特數(shù)模型對軸向主流的影響

圖7在圖6的基礎上添加了不同定常Prt的計算結(jié)果,據(jù)前文可發(fā)現(xiàn)Prt的影響更突出體現(xiàn)在對上壁溫的預測上,因此圖7重點展示水平圓管上壁溫的分布。由圖7可以看出,Prt的減小會使上壁溫大幅降低,其預測值的整體趨勢與實驗值更加吻合。

圖7 不同湍流普朗特數(shù)模型上壁溫預測值與實驗值

圖8給出了沿軸向方向不同截面處各Prt時計算得到的中垂線處速度輪廓,其中,圖8(a)為x=600 mm截面處的速度分布,r為徑向距離(某點到圓心的距離),R為半徑。從圖中可以發(fā)現(xiàn),不同Prt模型獲得的速度分布曲線底部幾乎重疊,但頂部已呈現(xiàn)較為明顯的差異。隨著加熱的持續(xù)進行,Prt對主流速度的影響進一步凸顯,頂母線近壁區(qū)的軸向速度差異繼續(xù)擴大。在x=700 mm截面處的頂母線近壁區(qū),速度隨Prt的增大而增加。其中,Tang模型所得速度輪廓在近上壁處出現(xiàn)第二速度峰值(半M型),這是Prt的增大使得流體傳熱受阻,浮升力效應影響進一步加劇的結(jié)果,即頂部積聚的“類氣態(tài)”sCO2愈發(fā)增多,二次流向下動量輸送中止,使得上半?yún)^(qū)域形成第二速度峰值。當sCO2流經(jīng)x=1 000 mm截面處時,除Prt=0.5時外,其他模型所得速度剖面均于近壁處變平,直至加熱段下游x=1 500 mm處,頂部流體速度已高于底部流體。

(a) x=600 mm

圖9給出了不同Prt模型所計算得到的各截面處湍動能分布。由圖9可以看出,Prt=0.5時在4個截面處的湍動能均呈現(xiàn)出底母線的比頂母線的高的特點,但Prt=0.85以及Tang模型、Du模型的載入使得湍動能在各截面呈現(xiàn)出一系列特殊現(xiàn)象。與圖8中截面的上母線近壁區(qū)速度分布相對應,Prt的增大使得該區(qū)域速度剖面變平,甚至在r/R=0.85處形成第二速度峰值。這雖使緊貼上壁處的速度梯度增大,湍動能強度繼而增大,但r/R=0.85處速度峰值帶來的零速度梯度直接導致該徑向湍動能最低值的出現(xiàn),該湍動能值對流動傳熱起著更為決定性的負面影響(這一點反映在圖9(a)中),因此傳熱惡化加劇。

(a) x=600 mm

2.3.2 湍流普朗特數(shù)模型對徑向邊界層的影響

由于動量傳遞和對流換熱主要受速度邊界層和熱邊界層的影響,以及Prt的影響集中體現(xiàn)在加熱段上游的“類液區(qū)”,因此選取加熱初始段x=600 mm處典型截面進行邊界層層面的深入分析。鑒于不同Prt對底母線近壁區(qū)的各參數(shù)影響并不顯著,因此下述分析將基于該截面中垂線處頂母線近壁區(qū)徑向參數(shù)的變化展開。

圖10(a)給出了近壁區(qū)的溫度輪廓。從圖10可以發(fā)現(xiàn),對于定常和非定常Prt,Prt的減小使得壁溫下降,同時熱邊界層厚度δt亦變薄,對流換熱得到增強。圖10(b)則展示了不同Prt模型計算下Pr的變化,由圖10(b)可以看出,隨著定常Prt的減小,Pr峰值變化不大(Pr與比熱容cp的變化近乎同步,峰值所在處為擬臨界點,Pr峰值附近為大比熱區(qū)),但大比熱區(qū)會朝緩沖層與黏性底層移動并靠近,邊界層內(nèi)(占據(jù)著絕大部分的熱阻)傳熱熱阻下降,這再次佐證了sCO2傳熱性能的上升。從圖10還可以看出,變Prt模型時大比熱區(qū)相較于定常Prt模型時更窄,Pr峰值亦更低,且均位于對數(shù)律區(qū)。而從圖10(c)中近壁區(qū)Prt分布可以發(fā)現(xiàn),針對Tang模型,A值的調(diào)整不僅會影響Prt峰值大小,還會影響Prt的分布。A值增大即模型緩沖層的Prt減小,會使得Prt整體輪廓沿無量綱距離y+左移,但峰值仍位于對數(shù)律區(qū)內(nèi),且與Pr的分布較為一致。但對于Du模型,Prt峰值與Pr峰值錯位分布,Prt峰值已進入緩沖層。圖10(d)展示了湍動能的分布,可以發(fā)現(xiàn)在y+<60緊貼壁面處,湍動能隨著Prt的增大而減小,3個定常Prt和3個變Prt模型均在y+≈60處完成了數(shù)值高低的轉(zhuǎn)換,而該處附近(對應圖8中近上壁速度剖面的峰值處)的湍動能大小對sCO2傳熱起著更為主導的作用。

(a) 溫度

3 結(jié) 論

(1) SSTk-ω低雷諾數(shù)湍流模型中Prt=0.85的設置在弱浮升力效應情形下能夠比較可靠且準確地預測sCO2在水平管內(nèi)的流動換熱,而在強浮升力效應下則會產(chǎn)生較大偏差。

(2) 現(xiàn)有典型的分層賦值變湍流普朗特數(shù)修正模型應用于水平sCO2在弱浮升力效應工況下的傳熱模擬收效甚微,其并不能顯著提升SSTk-ω模型的預測性能。

(3)Prt模型對sCO2湍流傳熱RANS計算具有深層次的影響,且對上壁區(qū)的影響顯著大于下壁區(qū)。在軸向主流方面,Prt的增大會加劇浮升力效應,這使得上壁區(qū)出現(xiàn)第二速度峰值,進而影響頂母線近壁區(qū)的湍動能分布。邊界層方面,Prt越小,上壁區(qū)熱邊界層越薄,大比熱區(qū)往壁面靠近,傳熱熱阻相應減小,壁溫預測值隨之下降。

(4) 后續(xù)針對水平sCO2傳熱計算Prt修正時,應更多聚焦于浮升力影響顯著的上壁區(qū),同時可在管徑修正因子和壓力修正因子的基礎上,加入q/G修正因子,以更準確地反映浮升力影響。