張 戈，張 昊，顧漢洋，楊燕華

(上海交通大學(xué) 核能科學(xué)與工程學(xué)院，上海 200240)

超臨界壓力下流體物性奇特的變化規(guī)律導(dǎo)致其對流傳熱現(xiàn)象十分獨特。Jackson等[1]和Cheng等[2]分別對超臨界流體對流傳熱進(jìn)行研究，根據(jù)傳熱機(jī)理不同，可分為3種傳熱現(xiàn)象：大流量小熱流下的強(qiáng)迫傳熱、小流量大熱流下的混合對流傳熱和大流量大熱流下的強(qiáng)加速流動傳熱。

當(dāng)流量較小、壁面熱流較大時，壁面附近流體密度隨溫度的升高而急劇減小，導(dǎo)致浮升力的作用加強(qiáng)，并對傳熱產(chǎn)生很大影響，形成傳熱弱化現(xiàn)象。根據(jù)浮升力作用方向不同，混合對流可分為浮升力助推流動(浮升力的方向和流動方向相同)和浮升力阻滯流動(浮升力的方向和流動方向相反)。浮升力作用下的傳熱弱化現(xiàn)象是超臨界流體傳熱的一重要特征，也是人們非常關(guān)心的一個現(xiàn)象。在相關(guān)文獻(xiàn)中，Sharabi等[3-4]、Wilcox[5]、Speziale等[6]用傳統(tǒng)的湍流模型模擬了傳熱弱化現(xiàn)象，其結(jié)果能定性預(yù)測傳熱弱化現(xiàn)象的發(fā)生與恢復(fù)，但在壁面溫度的預(yù)測上均偏高。Petukhov等[7]總結(jié)了浮升力對超臨界流動傳熱的影響，主要歸結(jié)為兩類：第1類是浮升力直接影響平均速度分布，進(jìn)而間接影響湍流流動，稱為外部效應(yīng)；第2類是浮升力直接影響湍流流動，稱為結(jié)構(gòu)效應(yīng)。

本工作通過實驗獲取10 mm單管中超臨界水流動的傳熱弱化數(shù)據(jù)，采用文獻(xiàn)[8]改進(jìn)的k-ε-kt-εt四方程模型對該實驗參數(shù)下的管內(nèi)流動傳熱特征進(jìn)行數(shù)值分析。

1 實驗研究

1.1 實驗

上海交通大學(xué)建成了適用于超臨界流體熱工水力實驗的SWAMUP回路[9]，SWAMUP回路示意圖示于圖1。該回路運行壓力最高30 MPa，出口溫度最高可達(dá)550 ℃，主泵流量5 t/h。預(yù)熱段加熱功率最大300 kW，試驗段加熱功率最大900 kW。穩(wěn)壓器通過頂部的氮氣空間保持回路整體的壓力。試驗段為內(nèi)徑10 mm的單管，其加熱長度為2.5 m。由布置在壁面上的49根間隔為5 cm、直徑為1 mm的熱電偶測量壁面溫度。實驗測量的各參數(shù)的不確定度列于表1。

圖1 SWAMUP回路示意圖

表1 測量參數(shù)的不確定度

1.2 實驗數(shù)據(jù)處理

試驗段的熱平衡實驗表明試驗段的保溫效果很好。因此，加熱的電功率可認(rèn)為全部加入到試驗段中。壁面熱流密度為均勻分布，因此內(nèi)側(cè)壁面熱流密度為：

q=UI/πDiL

(1)

外壁面由49根熱電偶測量溫度。將壁面劃分為100層，在假設(shè)體積發(fā)熱量相同的情況下，內(nèi)壁面溫度由99次的式(2)迭代計算得到。

/4kw)((Di/2)2-(Di+1/2)2)-

(qVl/2kw)(Di/2)2ln(Di/Di+1)

(2)

其中：tw為壁面溫度；i為壁面劃分層數(shù)的編號；kw為管道熱導(dǎo)率；qVl為體積熱流密度。qVl的表達(dá)式為：

(3)

管道熱導(dǎo)率由平均溫度算得，平均溫度為：

/2

(4)

主流溫度由沿流動方向的均勻焓升計算，傳熱系數(shù)由內(nèi)壁溫度和主流溫度計算。

實驗工況列于表2。其中，1#工況為低流量高熱流密度，2#工況為高流量高熱流密度。

表2 實驗工況

2 數(shù)值模擬分析

假設(shè)流動為二維軸對稱流動，其在柱坐標(biāo)下的質(zhì)量、動量方程為：

同樣是金枝玉葉的段譽，第一次來燕子塢吃的那些：“茭白蝦仁”“龍井茶葉雞丁”，看看就教人饞涎欲滴。段譽的當(dāng)時心理評判是這樣的：“魚蝦肉食之中混以花瓣鮮果，色彩既美，自別有天然清香。”

ρρ

(5)

(6)

(7)

其中：r為柱坐標(biāo)的半徑；x為柱坐標(biāo)的長度；ρ為密度；u、v為不同方向的速度；p為壓力；g為重力加速度；μe為有效黏性系數(shù)，μe=μ+μt，μ為動力黏性系數(shù)，μt為湍流黏性系數(shù)。

(8)

其中：Cμ為常數(shù)；fμ為考慮壁面效應(yīng)的衰減函數(shù)；k為湍動能；ε為湍動能耗散率。

能量方程為：

ρρ

(9)

其中：h為焓；σt為湍流普朗特數(shù)，本文取0.9。

k和ε方程為：

μ

ρPk+ρGk-ρε

(10)

(11)

根據(jù)文獻(xiàn)[9]對本文中的k-ε模型進(jìn)行改進(jìn)，改進(jìn)后的k-ε模型的參數(shù)和方程列于表3。

表3 k-ε模型參數(shù)和方程

(12)

(13)

(14)

式中：Ret為湍流雷諾數(shù)；y+為無量綱參數(shù)；Reb為主流雷諾數(shù)；εw為ε方程在壁面處的邊界條件；y為壁面法線方向；ν為運動黏度。

在k方程中，對于垂直管流動，Gk為：

λ

ρCPPt+ρCPε

(15)

(16)

CP2=0.62+1.2exp(-3.2Pr)

(17)

(18)

(19)

根據(jù)文獻(xiàn)[8]給出了改進(jìn)后的kt-εt模型，其模型參數(shù)列于表4。

表4 kt-εt模型參數(shù)和方程

計算中，徑向和軸向網(wǎng)格結(jié)構(gòu)根據(jù)計算工況進(jìn)行相應(yīng)調(diào)整，保證最靠近壁面的網(wǎng)格的y+<0.2。為得到黏性支層區(qū)域的詳細(xì)特征，在y+<5.0的區(qū)域內(nèi)布置40個網(wǎng)格點。對流項采用QUICK離散格式，速度場和壓力場采用SIMPLEC算法。超臨界物性根據(jù)APWS-95程序生成，因管道沿程壓降相對系統(tǒng)壓力很小，物性只是隨流體溫度變化。

3 結(jié)果分析與討論

3.1 壁面溫度的預(yù)測

圖2示出壁面溫度實驗測量與計算結(jié)果的比較。采用的對比模型是FLUENT中常用的低雷諾數(shù)模型，分別為YS模型[10]、AB模型[11]、CH模型[12]、AKN模型[13]、LS模型[14]和常用的兩方程SST模型[15]。這幾種模型均要求近壁面處具有較多的網(wǎng)格。圖2計算了兩種工況：工況1，p=23 MPa，G=596.6 kg/(m2·s)，q=772.6 kW/m2；工況2，p=25 MPa，G=2 021 kg/(m2·s)，q=1 385 kW/m2。由圖2a可見：YS、AB和CH模型所計算的壁面溫度較其他模型計算的均高出很多，其計算結(jié)果不能預(yù)測傳熱弱化現(xiàn)象的恢復(fù)；SST模型能定性預(yù)測到壁面溫度的第1個峰值，但定量上還是差了很多。而本文模型能很好預(yù)測傳熱弱化現(xiàn)象的產(chǎn)生及恢復(fù)，也能定量預(yù)測壁面溫度。由圖2b可見，傳統(tǒng)的低雷諾數(shù)模型也高估了壁面溫度，本文模型能捕捉到沿流動方向上的壁溫逐漸上升趨勢，而YS和AB模型的計算結(jié)果是下降的。

圖2 壁面溫度實驗測量與計算結(jié)果的比較

3.2 傳熱機(jī)理

以工況1為例，為了分析傳熱弱化現(xiàn)象，對5個不同位置處的物性和湍流特征進(jìn)行分析，結(jié)果示于圖3。圖3中，x=0對應(yīng)正常傳熱區(qū)域，x=10R和x=50R分別為傳熱弱化起始位置和弱化過程區(qū)域，x=80R和x=140R分別對應(yīng)傳熱弱化最嚴(yán)重區(qū)域和后期壁面換熱能力恢復(fù)區(qū)域，R為流道半徑。

圖3 工況1的近壁面區(qū)物性和湍流特征

從圖3a可看出，在近壁面區(qū)域出現(xiàn)了溫度的最大梯度。因超臨界流體的物性變化劇烈區(qū)域為擬臨界點附近區(qū)域，將擬臨界溫度±5 ℃稱為大物性變化(LPV)區(qū)。可見，在上游x=0位置的LPV區(qū)范圍十分狹小，而下游x=50R和x=80R位置的LPV區(qū)的范圍明顯擴(kuò)大。對應(yīng)不同位置的大比熱容變化區(qū)的分布特征，從圖3b可十分明顯地得看出這一點，且LPV區(qū)的中心位置開始遠(yuǎn)離壁面，而x=140R位置的LPV區(qū)的中心位置比x=80R更靠近壁面。對于常規(guī)流體，因近壁面區(qū)域的溫度變化最為劇烈，其熱物性也在該區(qū)域變化最劇烈，但超臨界流體的物性變化最劇烈的位置由LPV區(qū)的位置決定，這是超臨界流體的特殊性。由圖3c可見，在x=0和x=10R位置，壁面區(qū)域存在很大的密度變化，該區(qū)域存在顯著的浮升力，并顯著改變了流體速度分布，但其大變化區(qū)域處仍限制在黏性支層區(qū)域(y+<5)，且不可能對湍流特性產(chǎn)生顯著影響。在下游x=50R位置，密度劇烈變化區(qū)域拓寬到整個近壁面區(qū)域，壁面溫度的上升促使更下游x=80R和x=140R位置近壁面區(qū)域的流體密度進(jìn)一步降低。近壁面區(qū)域的低密度流體與中心區(qū)域的高密度流體之間的密度差產(chǎn)生近壁面區(qū)域的浮升力作用，顯著改變了該區(qū)域的流動和湍流特性。圖3d示出速度的分布特征。x=0和x=10R的密度僅局限在黏性支層，其速度表現(xiàn)為典型的強(qiáng)制對流湍流速度分布特征。在x=50R位置，浮升力的作用顯著增強(qiáng)，近壁面區(qū)域的流體速度明顯升高，使得整個速度分布變得平坦。而隨著浮升力作用的進(jìn)一步增強(qiáng)，近壁面區(qū)域的速度進(jìn)一步增加，在x=80R區(qū)域出現(xiàn)典型的混合對流的“M”型速度分布特征，而x=140R位置因湍流和壁面換熱的恢復(fù)，“M”型速度分布的近壁面最大值有所減低。

工況1下不同x軸位置的剪切應(yīng)力分布示于圖4。浮升力對湍流的影響分為兩部分。

1) 浮升力直接影響剪切應(yīng)力分布

由圖4a可見，當(dāng)浮升力作用較小時，總剪切應(yīng)力基本呈線性分布。當(dāng)有浮升力作用時，剪切應(yīng)力發(fā)生改變。因壁面熱流的輸入，近壁面處和主流區(qū)間存在溫差，產(chǎn)生如圖3c所示的密度差，導(dǎo)致近壁面處的速度變大，相對主流處的速度形成趨于平坦的層流化現(xiàn)象(圖3d)。同時，密度差還導(dǎo)致近壁面處的剪切應(yīng)力減小(圖4b、c)。剪切應(yīng)力和平均速度梯度對湍流傳熱的影響很大。平均速度梯度的存在導(dǎo)致較大渦的存在，大渦通過剪切應(yīng)力增強(qiáng)湍流的交混作用，提高熱量傳遞的效率。在浮升力作用下，剪切應(yīng)力和速度梯度變小，使得湍流能力變?nèi)?，并?dǎo)致傳熱能力變?nèi)酰M(jìn)而使近壁面處和主流處溫差進(jìn)一步加大，使得流體密度差進(jìn)一步加大，導(dǎo)致剪切應(yīng)力和速度梯度進(jìn)一步減小，形成一正反饋模式，導(dǎo)致傳熱能力急劇弱化。隨著傳熱能力的弱化，近壁面處和主流處的密度差進(jìn)一步增大，導(dǎo)致負(fù)的剪切應(yīng)力的出現(xiàn)(圖4d、e)，剪切應(yīng)力的絕對值增加。另一方面，近壁面處的速度變大，出現(xiàn)“M”型速度分布，近壁面處與主流處的速度梯度再次增加。這使得湍流能力加劇，傳熱能力開始恢復(fù)。

圖4 工況1不同x軸位置剪切應(yīng)力分布

2) 浮升力直接影響湍流結(jié)構(gòu)

湍流產(chǎn)生項和浮升力項是k-ε方程中的重要項。圖5示出工況1近壁面處的湍流特性。由圖5可見，湍流產(chǎn)生項在x=80R處達(dá)到最小值，而此處恰好是傳熱惡化最劇烈的點。當(dāng)流動狀態(tài)由湍流向?qū)恿鬓D(zhuǎn)變時，湍流產(chǎn)生項急劇變小，當(dāng)流動狀態(tài)由層流向湍流轉(zhuǎn)變時，湍流產(chǎn)生項又急劇增加。浮升力項在x=0位置很小，隨著流動狀態(tài)的變化有一定增加。但因為湍流產(chǎn)生項的減小，導(dǎo)致浮升力項的影響不能被忽略，成為湍動能的一主要源項。

圖5 工況1近壁面處的湍流特性

4 結(jié)論

1) 改進(jìn)的低雷諾數(shù)模型能準(zhǔn)確預(yù)測超臨界流體流動傳熱的傳熱惡化現(xiàn)象發(fā)生的起始、峰值和結(jié)束位置，而傳統(tǒng)模型均高估了壁面溫度。

2) LPV區(qū)域的位置影響浮升力效應(yīng)。當(dāng)LPV區(qū)在過渡區(qū)時，浮升力明顯影響了平均速度和湍流剪切應(yīng)力，導(dǎo)致了傳熱弱化現(xiàn)象的產(chǎn)生與恢復(fù)。

3) 除浮升力的外部效應(yīng)影響，其結(jié)構(gòu)效應(yīng)也會影響湍流產(chǎn)生項和浮升力項，而這兩項均是傳熱弱化情況下湍流方程中重要的項。

參考文獻(xiàn)：

[1] JACKSON J D, HALL W B. Forced convection heat transfer to fluids at supercritical pressure: Turbulence forced convection in channels and bundles[M]. New York: Hemisphere Publishing Corporation, 1978: 563-611.

[2] CHENG X, SCHULENBERG T. Heat transfer at supercritical pressures, literature review and application to an HPLWR[R]. Kalsruhe, Germany: FZKA, 2001.

[3] SHARABI M, AMBROSINI W, FORGIONE N, et al. Prediction of experimental data on heat transfer to supercritical water with two-equation turbulence models[C]∥3rd International Symposium on SCWR-Design and Technology. Shanghai, China: [s. n.], 2007.

[4] SHARABI M, AMBROSINI W. Discussion of heat transfer phenomena in fluids at super critical pressure with the aid of CFD models[J]. Annals of Nuclear Energy, 2009, 36(1): 60-71.

[5] WILCOX D C. Turbulence modeling for CFD[M]. California, USA: DCW Industries Inc., 2000.

[6] SPEZIALE C G, ABID R, ANDERSON E C. Critical evaluation of two-equation models for near wall turbulence[J]. AIAA Journal, 1992, 30(2): 324-331.

[7] PETUKHOV B S, POLYAKOVA H. Heat transfer in turbulent mixed convection[M]. New York: Hemisphere Publishing Corporation,1988.

[8] ZHANG H, XIE Z R. Numerical study on supercritical fluids flow and heat transfer buoyancy[C]∥The 8th International Topical Meeting on Nuclear Thermal-hydraulics, Operation and Safety (NUTHOS-8). Shanghai, China: [s. n.], 2010.

[9] GU H Y, ZHANG G, WEN Q L, et al. Supercritical water heat transfer test in a vertical tube[C]∥The 8th International Topical Meeting on Nuclear Thermal-hydraulics, Operation and Safety (NUTHOS-8). Shanghai, China: [s. n.], 2010.

[10] YANG Z, SHIH T H. New time scale basedk-εmodel for near-wall turbulence[J]. AIAA Journal, 1993, 31(7): 1 191-1 198.

[11] ABID R. Evaluation of two-equation turbulence models for predicting transitional flows[J]. Int Eng Sci, 1993, 31(6): 831-840.

[12] CHIEN K Y. Predictions of channel and boundary layer flows with a low Reynolds number turbulence model[J]. AIAA Journal, 1982, 20(1): 33-38.

[13] ABE K, KONDOH T, NAGANO N. A new turbulence model for predicting fluid flow and heat transfer in separating and reattaching flows, Ⅰ: Flow field calculation[J]. Int J Heat Mass Transfer, 1994, 37(1): 139-151.

[14] LAUNDER B E, SHARMA B I. Application of the energy dissipation model of turbulence to the calculation of flow near a spinning disk[J]. Lett Heat Mass Transfer, 1974, 1(2): 131-138.

[15] MENTER F R. Two equation eddy-viscosity turbulence models for engineering application[J]. AIAA Journal, 1994, 32(8): 1 598-1 605.