基于CFD方法的矩形窄縫通道內(nèi)過冷流動沸騰模型研究

劉 凱，李林峰，王明軍,*，章 靜，田文喜，秋穗正，蘇光輝

(1.西安交通大學(xué) 動力工程多相流國家重點實驗室，陜西 西安 710049；2.西安交通大學(xué) 核科學(xué)與技術(shù)學(xué)院，陜西 西安 710049)

板型燃料元件具有活性區(qū)比功率大、結(jié)構(gòu)緊湊等優(yōu)勢，在眾多研究堆、材料輻照堆中獲得了應(yīng)用。其堆芯冷卻劑通道截面形狀為矩形窄縫，具有寬高比大的幾何特征，在發(fā)生過冷沸騰時，汽泡運動行為將受到高度方向上更陡峭的速度分布的影響而發(fā)生改變。

在加熱壁面上的汽泡行為中，汽泡滑移現(xiàn)象受到很多學(xué)者的關(guān)注。Thorncroft和Klausner[1]對豎直通道中冷卻劑FC-87的沸騰換熱系數(shù)進行了測量，結(jié)果表明，向上流動沸騰的換熱系數(shù)顯著高于向下流動工況，分析認為向上流動時壁面附近汽泡的滑移運動增強了換熱。Li等[2]進行了運行壓力0.1 MPa、高2 mm的窄矩形通道內(nèi)豎直向上的流動沸騰實驗，同樣觀察到了頻繁的汽泡滑移。Yuan等[3-4]在運行壓力0.3～1 MPa工況的2 mm高窄縫通道實驗中觀察到了汽泡滑移現(xiàn)象，并對滑移汽泡的瞬態(tài)換熱項建模，結(jié)果表明滑移汽泡對換熱的貢獻可達靜止汽泡的10倍。Yoo等[5]使用熱成像記錄了滑移汽泡的影響面積，結(jié)果表明汽泡滑移過程中影響面積與汽泡直徑相關(guān)。

計算流體力學(xué)(CFD)目前已廣泛應(yīng)用于工程模擬中[6-7]，包括棒束通道的兩相沸騰模擬[8-9]，但對窄矩形通道內(nèi)沸騰分析時往往不能采用已有模型和關(guān)系式，因此，開發(fā)適用于矩形窄縫通道的沸騰模型具有相當?shù)闹匾浴?/p>

本文針對矩形窄縫通道中的汽泡滑移對沸騰換熱的影響，構(gòu)建包含滑移熱流的壁面熱流分配模型，并建立機理性的汽泡受力模型和滑移模型計算汽泡脫離直徑、浮升直徑和滑移距離等輔助參數(shù)，開發(fā)一套適用于矩形窄縫通道內(nèi)豎直向上流動沸騰的壁面沸騰模型。選用Nuthel窄縫通道沸騰實驗進行數(shù)值模擬，通過對比壁面過熱度的計算結(jié)果與實驗結(jié)果，驗證本文模型對矩形通道內(nèi)中低壓過冷沸騰流動換熱特性模擬的有效性。

1 數(shù)學(xué)物理模型

本文基于歐拉兩流體六方程模型，分別對氣、液兩相建立質(zhì)量、動量、能量守恒方程，并考慮相間相互作用，引入相間質(zhì)量、動量、能量交換模型[10]。相間動量交換中的作用力包括曳力和非曳力，非曳力又包含升力、湍流耗散力、壁面潤滑力，相間質(zhì)量和能量交換包括近壁面沸騰換熱和相界面上的蒸發(fā)、冷凝過程。

此外，本文考慮汽泡滑移現(xiàn)象對沸騰換熱的影響，構(gòu)建包括滑移熱流的壁面沸騰模型，并建立機理性的輔助模型計算相關(guān)汽泡關(guān)鍵參數(shù)，以對模型整體進行封閉，形成壁面熱流密度和壁面溫度的計算關(guān)系。

1.1 考慮汽泡滑移的壁面沸騰模型

Kurul和Podowski[11]提出的RPI壁面沸騰模型廣泛應(yīng)用于流動沸騰計算，該模型僅考慮了汽化核心當?shù)氐钠菪袨?，認為壁面熱流包括對流換熱、蒸發(fā)、淬滅熱流3部分。

本文針對窄矩形通道中的汽泡滑移行為，采用以下過程描述汽泡生長的周期，如圖1所示。汽泡在脫離汽化核心后將垂直于壁面浮升或平行于壁面滑移。滑移過程中由于其與壁面接觸處的微液膜蒸發(fā)，并與其他汽化核心處汽泡融合，汽泡直徑持續(xù)增大，當達到一定值后汽泡浮升進入主流區(qū)域。該過程中，汽泡離開汽化核心處、開始滑移時的直徑為脫離直徑dd，汽泡垂直壁面運動、遠離壁面并進入主流時的直徑為浮升直徑dl。

圖1 壁面沸騰模型對汽泡生長周期的描述Fig.1 Description of bubble period by wall boiling model

圖2 壁面沸騰模型對總熱流的分配Fig.2 Distribution of total heat flux by wall boiling model

1) 汽泡生長階段(0～tg，其中tg為汽泡生長時間)中，汽化核心區(qū)域的相變和滑移影響區(qū)域的微液層蒸發(fā)計入蒸發(fā)熱流。

2) 滑移影響階段(tg～t*，其中t*為換熱增強效果消失時刻)中，在汽化核心區(qū)域，固體側(cè)由淬滅熱流主導(dǎo)，而流體側(cè)為由于汽泡擾動的滑移熱流。在滑移影響區(qū)域，由汽泡滑移的擾動帶來的額外換熱增強計入滑移熱流。

3) 在流場恢復(fù)階段(t*～1/f，其中f為汽泡脫離頻率)中，單相對流換熱重新占據(jù)主導(dǎo)，而其他區(qū)域在整個汽泡周期中僅存在單相對流換熱過程。

據(jù)此，本文所開發(fā)的壁面沸騰模型將壁面總熱流qt分為4項熱流，分別為單相對流熱流qc、蒸發(fā)熱流qe、淬滅熱流qq和滑移熱流qs，其表達式為：

qt=qc+qe+qq+qs

(1)

蒸發(fā)熱流的表達式為：

(2)

式中：hfg為汽化潛熱；f為汽泡脫離頻率；N*為汽化核心密度。

由于汽泡滑移擾動，來自主流區(qū)域的較低溫度流體補充汽泡原先占據(jù)的位置，并與過熱壁面接觸換熱。將其簡化為半無限大區(qū)域冷卻平板瞬態(tài)導(dǎo)熱問題，則壁面瞬態(tài)熱流表達式為：

(3)

式中：kl為液相導(dǎo)熱系數(shù)；Tw為壁面溫度；Tl為流體溫度，此處取液相近壁面溫度；ηl為液相熱擴散率；τ為瞬態(tài)過程時間，即滑移熱流的影響時間，從汽泡滑移開始至換熱增強效果消失，即瞬態(tài)導(dǎo)熱系數(shù)與單相對流換熱系數(shù)相等時為止，記作0～τ*。令瞬態(tài)導(dǎo)熱系數(shù)與單相對流換熱系數(shù)相等，得到滑移影響時間表達式：

(4)

式中，hfc為當?shù)亓鲃踊謴?fù)時的單相對流換熱系數(shù)。

考慮汽泡影響時間和影響面積，積分可得滑移過程中的瞬態(tài)導(dǎo)熱熱流，即滑移熱流：

(5)

式中，asl為汽泡的影響面積。

單相對流熱流根據(jù)近壁面對流換熱系數(shù)hfc計算：

qc=hfc(Tw-Tl)·(1-aslN*+aslN*(1-τ*f))

(6)

Gilman等[12]認為，汽泡生長時底部存在高溫干斑；當汽泡脫離后這部分壁面顯熱會釋放到流體中，組成新的淬滅熱流。具體地，認為汽泡底部干斑呈圓形，其直徑為汽泡脫離直徑的1/2。干斑下的半球形高溫區(qū)較其他區(qū)域溫度高ΔTh，淬滅點與周圍壁面的溫差約為3 K。因此淬滅熱流表達式為：

(7)

1.2 輔助模型

1) 汽泡受力模型

由于現(xiàn)有經(jīng)驗關(guān)系的適用范圍有限或未區(qū)分汽泡脫離直徑和浮升直徑，本文建立汽泡受力模型計算汽泡脫離直徑和浮升直徑。采用Sugrue和Buongiorno[13]推薦使用的受力表達式和系數(shù)組合列于表1，其在反應(yīng)堆典型應(yīng)用場景內(nèi)得到了一定的驗證。1個傾斜壁面上汽泡受力情況如圖3[13]所示，圖中x方向平行于壁面為流動方向，y方向垂直于壁面為浮升方向，θ為壁面水平傾角。

表1 汽泡所受作用力表達式Table 1 Expression for force acting on bubble

圖3 單個靜止汽泡受力示意圖Fig.3 Diagram of force on single static bubble

對于靜止汽泡，兩方向上的合力分別為：

(8)

式中：θ為加熱壁面的傾角；φ為汽泡偏向流動方向的傾角，一般假設(shè)汽泡生長方向與壁面法向的夾角為π/15。當ΣFx>0時，汽泡達到dd并平行于壁面運動；當ΣFy>0時，汽泡達到dl并垂直于壁面運動。

2) 汽泡滑移模型

通過滑移過程中汽泡生長方程計算汽泡滑移時間，進而得到汽泡滑移距離，本文使用Maity[14]擬合的汽泡直徑關(guān)系式計算，如下：

(9)

式中：d1和d2分別為汽泡滑移起始和終止時的直徑；τs為滑移時間；Jasup和Jasub分別為壁面過熱度和液相過冷度雅各比數(shù)；Reb為以汽泡平均直徑、汽泡滑移速度計算的雷諾數(shù)。

根據(jù)Basu[15]對滑移汽泡的建模，在汽泡由上一個汽化核心與當?shù)仄萑诤喜⑾蛳乱粋€汽化核心滑移過程中，汽泡直徑由dN-1不斷增大，在該過程中判斷汽泡是否脫離即可得到汽泡滑移距離：

l=Nms+ldN-1～dl

(10)

滑移過程中單個汽泡的影響面積為：

asl=l(dd+dl)/2

(11)

此外，當浮升直徑小于汽泡脫離直徑時，汽泡將不發(fā)生滑移而直接浮升。此時，汽泡行為實際上與RPI模型的描述相同，汽泡影響面積按下式計算：

(12)

式中，K為影響范圍因子。根據(jù)Valle和Kenning等[16]對汽化核心處汽泡的觀測分析使用以下公式計算：

(13)

3) 汽泡脫離頻率模型

汽泡脫離頻率使用Cole[17]關(guān)系式計算，其表達式如下：

(14)

式中，g為重力加速度。

4) 汽化核心密度模型

汽化核心密度使用Lemmert-Chawla[18]公式計算：

Nw=[n0(Tw-Tsat)]m

(15)

式中，參考密度n0和指數(shù)m分別取210和1.805。根據(jù)汽泡滑移距離，汽泡融合個數(shù)可以按l/s估計，則修正后的汽化核心密度為：

(16)

2 實驗?zāi)M與結(jié)果分析

2.1 實驗介紹

選取西安交通大學(xué)反應(yīng)堆熱工水力研究室(Nuthel)梁振輝[19]的過冷沸騰實驗進行模型驗證。實驗回路主要由屏蔽泵、穩(wěn)壓器、預(yù)熱器、實驗段、直流加熱電源、冷卻器和相應(yīng)的數(shù)據(jù)測量系統(tǒng)組成。實驗段為豎直放置的窄矩形通道，寬度為50 mm，高度為1.5 mm和2.5 mm。實驗段的加熱方式為電加熱，有效加熱長度為0.8 m。工質(zhì)是去離子水，實驗壓力范圍0.5～5 MPa，入口過冷度范圍50～80 K，質(zhì)量流量13.3～266.7 kg/(m2·s)，熱流密度范圍10～300 kW/m2。實驗段外壁面上點焊熱電偶，測量外壁溫度計算得到內(nèi)壁溫度。焊點選擇在矩形壁面寬邊的中心線位置，布置間距為100 mm。實驗段參數(shù)列于表2[19]。本文選取截面尺寸為2.5 mm×5 mm通道內(nèi)4組不同壓力下的實驗工況進行模擬，具體運行參數(shù)列于表3。

2.2 模擬方法

1) 計算區(qū)域選取

由于窄縫通道具有較大的寬高比，除接近邊緣處，在寬度方向上參數(shù)分布基本均勻。為降低計算量，通過選取通道中心的縱截面進行二維模擬的方法簡化。為驗證該簡化的合理性，對工況Nuthel-659分別進行三維和二維計算，計算條件保持一致。

三維與二維計算結(jié)果的對比如圖4所示，可見，二維與三維計算得到的中心縱截面上的壁面溫度計算結(jié)果相差不大；除邊角非加熱區(qū)域外，三維計算結(jié)果的空泡份額分布在寬度方向上分布均勻。因此采用通道中心縱截面進行二維模擬的簡化方法可行。

表2 Nuthel窄縫通道沸騰實驗段參數(shù)Table 2 Parameter of Nuthel narrow rectangular channel boiling experiment section

表3 模擬工況運行參數(shù)Table 3 Operating parameter of simulation case

a——壁面中心線溫度分布對比；b——三維模型出口處空泡份額分布圖4 二維模型有效性驗證Fig.4 Validation of 2D model

2) 網(wǎng)格劃分

使用四邊形結(jié)構(gòu)化網(wǎng)格對矩形流道進行空間離散。為進行網(wǎng)格敏感性分析，劃分7套不同網(wǎng)格，具體劃分方式和計算結(jié)果列于表4，測試工況是Nuthel-627號。

對比網(wǎng)格#1、#2、#5和#7的結(jié)果可知，高度方向均勻劃分5個控制體時計算發(fā)散，將節(jié)點分布改為非均勻布置，使近壁面網(wǎng)格厚度為0.25 mm時則計算收斂；當劃分10或15個控制體時結(jié)果較為接近，但后者壁面網(wǎng)格y+值低于10，第1層節(jié)點落在速度過渡層內(nèi)，不適用標準壁面函數(shù)。因此，計算中高度方向劃分10個控制體。對比網(wǎng)格#3～#7的結(jié)果可知，長度方向上網(wǎng)格數(shù)量對壁面過熱度影響較小。綜上，選取#3網(wǎng)格對Nuthel實驗進行模擬。

表4 網(wǎng)格敏感性分析Table 4 Grid sensitivity analysis

3) 材料物性

計算中使用運行壓力下的水物性。由于實驗中入口過冷度較高，因此液相物性使用入口溫度和飽和溫度兩點線性插值方法確定。過冷沸騰中氣相溫度幾乎維持在飽和溫度，因此氣相物性使用飽和狀態(tài)下的水蒸氣物性。

本文模型中的淬滅熱流考慮了壁面材料的影響，而兩實驗中沸騰區(qū)壁面材料均為鋼，因此取不銹鋼的物性典型值，密度8 000 kg/m3，比熱500 kJ/(kg·K)。

4) 邊界條件

入口邊界為流量入口，給定液相質(zhì)量流量和溫度，氣相體積份額為0；出口為壓力出口，表壓為0。入口和出口非加熱段為絕熱邊界；認為加熱功率均勻分布，即加熱段為均勻熱流密度邊界。

2.3 計算結(jié)果分析

1) 壁面過熱度和近壁面空泡份額

分別使用本文模型和RPI模型進行模擬，得到壁面溫度和近壁面空泡份額的計算結(jié)果如圖5所示。

圖5 壁面溫度和近壁面空泡份額計算結(jié)果Fig.5 Calculation results of wall temperature and near-wall void fraction

由于過冷沸騰階段的高換熱效率，壁面溫度沿流動方向出現(xiàn)拐點，壁溫上升減緩，出現(xiàn)一段壁溫平臺，近壁面空泡份額相應(yīng)上升。通道壁面自沸騰起始點的平均過熱度計算結(jié)果與實驗結(jié)果的對比列于表5，可見RPI模型的計算結(jié)果偏低，而本文模型計算得到的沸騰區(qū)過熱度與實驗值符合良好，對壁面過熱度的預(yù)測精度較高。此外，本文模型近壁面空泡份額的計算值較低。

單相對流區(qū)的壁溫計算結(jié)果存在較大偏差，導(dǎo)致沸騰起始點的預(yù)測不準確。分析認為偏差主要由熱效率不穩(wěn)定造成：單相對流區(qū)的壁面溫度與壁面熱流密度約為線性關(guān)系，預(yù)先通過單相換熱實驗測量的熱效率可能在沸騰實驗過程中改變，因此模擬中的給定熱流與實驗條件存在偏差。而在過冷沸騰區(qū)域熱流偏差的影響較小，其原因為：過冷沸騰區(qū)壁面過熱度與壁面熱流約為0.25～0.3次冪函數(shù)關(guān)系，例如Jens-Lottes公式中ΔTsup∝q0.25，相同的熱流偏差在過冷沸騰區(qū)域引起的壁面過熱度誤差較小，僅為單相區(qū)域的2.4%。

2) 熱流密度分配

熱流密度分配的計算結(jié)果如圖6所示。本文模型的計算結(jié)果中，滑移熱流和單相對流熱流占據(jù)主導(dǎo)，蒸發(fā)熱流份額相對較小。而RPI模型的計算結(jié)果中，蒸發(fā)熱流在通道后段占據(jù)主導(dǎo)，單相對流熱流和淬滅熱流份額都相對較小。

表5 壁面平均過熱度計算結(jié)果與實驗結(jié)果的誤差Table 5 Error between calculated and experimental results of average superheat of wall

圖6 熱流密度分配計算結(jié)果Fig.6 Calculation result of wall heat flux distribution

本文模型的蒸發(fā)熱流份額較小，導(dǎo)致近壁面蒸發(fā)率較低，因此近壁面空泡份額的計算結(jié)果較小。此外，本文模型計算結(jié)果中淬滅熱流所占份額較小，幾乎可以忽略，說明在模擬工況中壁面顯熱的影響較低。

3 結(jié)論

本文針對矩形窄縫通道中的汽泡行為，考慮汽泡滑移現(xiàn)象對沸騰換熱的影響，構(gòu)建了包含滑移熱流的壁面熱流分配模型，并建立機理性的汽泡受力模型和滑移模型計算汽泡脫離直徑、浮升直徑和滑移距離等輔助參數(shù)，開發(fā)了一套適用于矩形窄縫通道內(nèi)向上流動沸騰的壁面沸騰模型，選用Nuthel窄縫通道沸騰實驗的1～4 MPa中低壓過冷沸騰工況進行驗證，主要結(jié)論如下：

1) 本文模型的壁面溫度計算結(jié)果相比RPI模型較高，模擬工況中沸騰區(qū)的壁面平均過熱度誤差由最大80%降至最大17%，說明滑移熱流的引入可以更精確地預(yù)測矩形窄縫通道內(nèi)向上流動沸騰的壁面溫度；

2) 本文模型的熱流分配計算結(jié)果中，滑移熱流和單相對流熱流份額占據(jù)主導(dǎo)，蒸發(fā)熱流份額較低，蒸發(fā)率較低，因此近壁面空泡份額計算結(jié)果相比RPI模型較低；淬滅熱流所占份額較小，說明在模擬工況中壁面顯熱的影響較低。