湍流普朗特數模型對水平圓管內sCO2傳熱預測的影響

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

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

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

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

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

1 物理模型與數值方法

1.1 幾何模型

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

表1 計算模型幾何參數

圖1 物理模型

1.2 控制方程

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

連續性方程:

(1)

動量方程:

(2)

能量方程:

(3)

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

湍動能k輸運方程:

(4)

比耗散率ω輸運方程:

(5)

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

更多模型細節可參考文獻[16]。

1.3 典型湍流普朗特數修正模型

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

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

(6)

式中:A為可調節的常數,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 網格劃分與無關性驗證

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

圖2 計算域橫截面網格

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

圖3 網格獨立性驗證

1.5 邊界條件與計算方法

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

表2 驗證工況

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

2 結果與討論

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

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

2.1.1 弱浮升力效應

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

(a) 工況TS1

2.1.2 強浮升力效應

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

(a) 工況TS2

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

由于強浮升力效應下Prt=0.85模型的預測效果不佳,本節嘗試引入典型Prt修正模型,結果如圖6所示。由圖6(a)和圖6(b)可知,傳熱惡化嚴重時,所嵌入的變Prt修正模型并未能改善預測性能。其中,針對圖6(b)中內徑為8 mm的水平圓管,Du模型計算值在趨勢上與Prt=0.85時的近乎一致,上壁溫預測值略偏高(受限于構造形式,Du模型無法應用于4 mm管徑的修正)。然而,Tang模型在溫度預測上的表現比Prt=0.85時更差,尤其是在入口處出現非常尖銳的上壁溫峰值,該峰值隨著A增大而降低。這是因為邊界層緩沖層中Prt=0.85+Pr/A(Pr恒為正值),A越大,緩沖層的Prt越小。可以設想,當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模型在黏性底層、過渡層和湍流核心區的湍流普朗特數均不小于0.85,而根據能量方程(3),湍流普朗特數所在項為湍流擴散項,湍流普朗特數的增大使得湍流摻混效應減弱,傳熱性能下降,繼而壁溫升高,這與文獻[7]和文獻[22]的結論一致。

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

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

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

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

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

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

(a) x=600 mm

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

(a) x=600 mm

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

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

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

(a) 溫度

3 結 論

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

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

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

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