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

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

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

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

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

1 實驗研究

1.1 實驗

上海交通大學建成了適用于超臨界流體熱工水力實驗的SWAMUP回路[9]，SWAMUP回路示意圖示于圖1。該回路運行壓力最高30 MPa，出口溫度最高可達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ù)處理

試驗段的熱平衡實驗表明試驗段的保溫效果很好。因此，加熱的電功率可認為全部加入到試驗段中。壁面熱流密度為均勻分布，因此內(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的表達式為：

(3)

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

/2

(4)

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

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

表2 實驗工況

2 數(shù)值模擬分析

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

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

ρρ

(5)

(6)

(7)

其中：r為柱坐標的半徑；x為柱坐標的長度；ρ為密度；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ù)文獻[9]對本文中的k-ε模型進行改進，改進后的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ù)文獻[8]給出了改進后的kt-εt模型，其模型參數(shù)列于表4。

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

計算中，徑向和軸向網(wǎng)格結(jié)構(gòu)根據(jù)計算工況進行相應(yīng)調(diào)整，保證最靠近壁面的網(wǎng)格的y+<0.2。為得到黏性支層區(qū)域的詳細特征，在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 傳熱機理

以工況1為例，為了分析傳熱弱化現(xiàn)象，對5個不同位置處的物性和湍流特征進行分析，結(jié)果示于圖3。圖3中，x=0對應(yīng)正常傳熱區(qū)域，x=10R和x=50R分別為傳熱弱化起始位置和弱化過程區(qū)域，x=80R和x=140R分別對應(yīng)傳熱弱化最嚴重區(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ū)的范圍明顯擴大。對應(yīng)不同位置的大比熱容變化區(qū)的分布特征，從圖3b可十分明顯地得看出這一點，且LPV區(qū)的中心位置開始遠離壁面，而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ū)域的流體密度進一步降低。近壁面區(qū)域的低密度流體與中心區(qū)域的高密度流體之間的密度差產(chǎn)生近壁面區(qū)域的浮升力作用，顯著改變了該區(qū)域的流動和湍流特性。圖3d示出速度的分布特征。x=0和x=10R的密度僅局限在黏性支層，其速度表現(xiàn)為典型的強制對流湍流速度分布特征。在x=50R位置，浮升力的作用顯著增強，近壁面區(qū)域的流體速度明顯升高，使得整個速度分布變得平坦。而隨著浮升力作用的進一步增強，近壁面區(qū)域的速度進一步增加，在x=80R區(qū)域出現(xiàn)典型的混合對流的“M”型速度分布特征，而x=140R位置因湍流和壁面換熱的恢復(fù)，“M”型速度分布的近壁面最大值有所減低。

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

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

由圖4a可見，當浮升力作用較小時，總剪切應(yī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)力增強湍流的交混作用，提高熱量傳遞的效率。在浮升力作用下，剪切應(yīng)力和速度梯度變小，使得湍流能力變?nèi)?，并?dǎo)致傳熱能力變?nèi)?，進而使近壁面處和主流處溫差進一步加大，使得流體密度差進一步加大，導(dǎo)致剪切應(yīng)力和速度梯度進一步減小，形成一正反饋模式，導(dǎo)致傳熱能力急劇弱化。隨著傳熱能力的弱化，近壁面處和主流處的密度差進一步增大，導(dǎo)致負的剪切應(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處達到最小值，而此處恰好是傳熱惡化最劇烈的點。當流動狀態(tài)由湍流向?qū)恿鬓D(zhuǎn)變時，湍流產(chǎn)生項急劇變小，當流動狀態(tài)由層流向湍流轉(zhuǎn)變時，湍流產(chǎn)生項又急劇增加。浮升力項在x=0位置很小，隨著流動狀態(tài)的變化有一定增加。但因為湍流產(chǎn)生項的減小，導(dǎo)致浮升力項的影響不能被忽略，成為湍動能的一主要源項。

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

4 結(jié)論

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

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

3) 除浮升力的外部效應(yīng)影響，其結(jié)構(gòu)效應(yīng)也會影響湍流產(chǎ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.