秦 浩，王明軍，李林峰，田文喜，蘇光輝，秋穗正

(西安交通大學(xué) 陜西省先進(jìn)核能技術(shù)重點(diǎn)實(shí)驗(yàn)室，動(dòng)力工程多相流國(guó)家重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710049)

過(guò)冷流動(dòng)沸騰指主流中液相溫度低于飽和溫度，而近壁面處發(fā)生局部沸騰的兩相流動(dòng)換熱現(xiàn)象。相比于單相對(duì)流換熱，過(guò)冷沸騰具有較高的換熱效率，因此在能源與動(dòng)力工程等行業(yè)中被廣泛應(yīng)用。國(guó)內(nèi)外學(xué)者均在探索基于計(jì)算流體動(dòng)力學(xué)(CFD)方法來(lái)研究?jī)上嗔鲃?dòng)換熱問(wèn)題[1-2]。Krepper等[3]采用Bartolomei等[4]開(kāi)展的壓力范圍在1.5～15 MPa間的實(shí)驗(yàn)數(shù)據(jù)對(duì)商用軟件CFX中的過(guò)冷沸騰模型進(jìn)行了驗(yàn)證。Cheung等[5]研究了核化密度、氣泡脫離直徑和氣泡脫離頻率等不同模型組合對(duì)低壓過(guò)冷流動(dòng)沸騰實(shí)驗(yàn)的適用性，發(fā)現(xiàn)任一經(jīng)驗(yàn)關(guān)系式組合的計(jì)算結(jié)果都不能與所有低壓實(shí)驗(yàn)數(shù)據(jù)良好符合。Zhang等[6]研究了不同湍流模型對(duì)數(shù)值模擬結(jié)果的影響，認(rèn)為k-ε模型較k-ω模型效果好。陳二峰等[7]和李松蔚等[8]修正了CFX中的氣泡脫離直徑模型，模擬了低壓工況下的過(guò)冷沸騰。

目前主流的商用CFD軟件中均包含兩相流的求解模塊，可基于歐拉兩流體模型或VOF(體積分?jǐn)?shù))模型求解兩相流場(chǎng)。但商用軟件封裝嚴(yán)格，開(kāi)放程度弱，對(duì)用戶的自主開(kāi)發(fā)有明顯的限制，在求解復(fù)雜的兩相流問(wèn)題時(shí)，很多模型難以修改或添加，從而不一定能滿足用戶的科研需求。而采用C++語(yǔ)言編寫的開(kāi)源CFD類庫(kù)OpenFOAM以其面向?qū)ο缶幊?、大?guī)模并行能力及開(kāi)放的編程環(huán)境等優(yōu)勢(shì)逐步在學(xué)術(shù)及工程領(lǐng)域得到廣泛應(yīng)用。目前，OpenFOAM在海洋與船舶工程、化學(xué)工程等方向的數(shù)值模擬已取得了良好的效果?；谠撻_(kāi)源平臺(tái)添加模型或開(kāi)發(fā)新的求解器以模擬過(guò)冷流動(dòng)沸騰現(xiàn)象將是兩相流研究的一個(gè)重要方向。

本文以4.5 MPa下豎直圓管內(nèi)過(guò)冷流動(dòng)沸騰現(xiàn)象為研究對(duì)象，闡述數(shù)值模擬所采用的數(shù)學(xué)物理模型，包括歐拉兩流體模型及相關(guān)輔助模型，基于OpenFOAM平臺(tái)進(jìn)行模擬。

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

本文基于歐拉兩流體六方程模型，考慮了氣液兩相間質(zhì)量、動(dòng)量和能量的傳遞，并引入壁面熱流密度分配模型以描述壁面沸騰現(xiàn)象。

1.1 相間動(dòng)量交換模型

由牛頓第三定律可知，氣相對(duì)液相的作用力Fgl與液相對(duì)氣相的作用力Flg大小相等、方向相反，即：

Flg=-Fgl

(1)

因此本節(jié)僅對(duì)液相對(duì)氣相的作用力進(jìn)行描述。彌散在液相中的氣泡受到來(lái)自液相的作用力可分為曳力Md和非曳力，非曳力又包括升力Ml、壁面潤(rùn)滑力Mwl、湍流耗散力Mtd和虛擬質(zhì)量力Mvm，則氣泡總的受力情況可表示為：

M=Md+Ml+Mwl+Mtd+Mvm

(2)

1.2 壁面熱流密度分配模型

廣泛應(yīng)用于流動(dòng)沸騰模擬的壁面熱流密度分區(qū)模型是RPI模型[9]，它由提出者Kurul和Podowski就職的倫斯勒理工學(xué)院(RPI)而得名。該模型將從壁面?zhèn)鬟f到流體中的熱流密度分為：1) 單相過(guò)冷流體強(qiáng)迫對(duì)流換熱Qc；2) 氣泡脫離壁面時(shí)單相過(guò)冷液體重新覆蓋壁面而引入的淬滅熱流Qq；3) 液相蒸發(fā)引入的熱流Qe。因此，壁面總熱流Qtot表示為：

Qtot=Qc+Qq+Qe

(3)

但上述表達(dá)式未考慮氣泡覆蓋壁面時(shí)壁面與氣相之間的單相對(duì)流換熱Qv。為完善物理模型，提高模擬的準(zhǔn)確性，改進(jìn)的熱流密度分配模型逐漸得到了廣泛應(yīng)用。

Qtot=f(αl)(Qc+Qq+Qe)+(1-f(αl))Qv

(4)

式中，f(αl)為關(guān)于液相份額αl的經(jīng)驗(yàn)關(guān)系式，本文采用Lavieville等[10]提出的模型：

(5)

式中，αl,crit為液相份額的臨界值，本文取為0.2。f(αl)函數(shù)的圖像示于圖1，當(dāng)液相份額低于0.2即空泡份額大于0.8時(shí)，f(αl)的值迅速減小，壁面熱流密度主要表現(xiàn)為氣相單相對(duì)流換熱。

單相液體的對(duì)流換熱Qc表達(dá)式為：

Qc=(1-Aw)hc(Tw-Tl)

(6)

式中：Aw為氣泡影響的區(qū)域面積；Tw為壁面溫度；Tl為主流液體溫度；hc為單相液體對(duì)流換熱系數(shù)。

圖1 f(αl)函數(shù)圖像Fig.1 f(αl) function graph

壁面驟冷過(guò)程中的換熱Qq可表示為：

Qq=Awhq(Tw-Tl)

(7)

式中，hq為驟冷換熱系數(shù)。

(8)

式中：f為氣泡脫離頻率；twait為氣泡等待時(shí)間；kl為液相導(dǎo)熱系數(shù)；ρl為液相密度；cpl為液相比熱容。氣泡等待時(shí)間twait可近似表示為：

twait=0.8/f

(9)

蒸汽產(chǎn)生帶走的熱量Qe通過(guò)汽化過(guò)程中的傳質(zhì)速率表示：

Qe=mwhlg

(10)

式中：hlg為汽化潛熱；mw為汽化的傳質(zhì)速率，其經(jīng)驗(yàn)關(guān)系式為：

(11)

式中：ρg為氣相密度；dw為氣泡脫離直徑；N為核化密度。

由于換熱機(jī)理不同，單位加熱面積可劃分為核態(tài)沸騰氣泡影響區(qū)Ab和非影響區(qū)1-Ab[11]。Valle等[12]研究認(rèn)為，氣泡影響面積份額可寫成單位加熱壁面上核化數(shù)目和氣泡投影面積的表達(dá)式：

(12)

式中，K為沸騰模型常系數(shù)。Valle等[12]推薦的K計(jì)算關(guān)系式為：

K=4.8e-Ja/80

(13)

Ja=ρlcpl(Tsat-Tl)/ρghlg

(14)

式中，Tsat為液相飽和溫度。

1.3 輔助模型

1) 氣泡脫離直徑模型

氣泡從加熱壁面上脫離時(shí)的直徑稱為氣泡脫離直徑，其大小取決于氣泡長(zhǎng)大過(guò)程中的受力情況。本文采用Tolubinsky等[13]根據(jù)實(shí)驗(yàn)觀測(cè)得出的經(jīng)驗(yàn)?zāi)Ｐ陀?jì)算氣泡脫離直徑dw。

dw=max(drefexp(-(Tsat-Tl)/Tref)，dmax)

(15)

式中：參考直徑dref和氣泡最大直徑dmax分別取0.6 mm和1.4 mm；參考溫度Tref通常取45 K。

2) 壁面核化密度模型

單位加熱壁面面積上氣化核心的數(shù)目稱為核化密度。本文采用Lemmert-Chawla模型[14]計(jì)算，該模型認(rèn)為核化密度N與壁面過(guò)熱度ΔTsup有關(guān)：

N=(CNΔTsup)n

(16)

式中：參數(shù)CN一般取210；n=1.805。

3) 氣泡脫離頻率模型

氣泡脫離頻率是加熱壁面上氣泡脫離速度的量度。Cole等[15]通過(guò)高速攝像技術(shù)研究了臨界熱流密度附近的沸騰現(xiàn)象，并推薦氣泡脫離頻率f的關(guān)系式為：

(17)

2 算例與結(jié)果分析

2.1 算例描述

Bartolomei等[16]于1967年設(shè)計(jì)搭建了豎直加熱圓管內(nèi)流動(dòng)沸騰實(shí)驗(yàn)回路，進(jìn)行了不同壓力、熱流密度、質(zhì)量流量以及入口過(guò)冷度下的過(guò)冷沸騰實(shí)驗(yàn)，測(cè)量獲得了壁面溫度、流體截面平均溫度、截面平均空泡份額等沿圓管高度的分布。文獻(xiàn)[16]中公開(kāi)的實(shí)驗(yàn)段及實(shí)驗(yàn)條件示于圖2。本文將以此實(shí)驗(yàn)為研究對(duì)象，進(jìn)行數(shù)值模擬。

圖2 實(shí)驗(yàn)段示意圖Fig.2 Schematic diagram of experimental section

2.2 計(jì)算方法

為減少計(jì)算量，本文取計(jì)算區(qū)域?yàn)槎S模型，即取圓管通道的1個(gè)豎直矩形截面，上下面分別為速度入口和壓力出口邊界條件，左側(cè)為對(duì)稱邊界條件，右側(cè)為加熱面。液相考慮其物性隨溫度的變化，氣相采用飽和溫度下的物性值。曳力采用Schiller-Naumann模型，壁面潤(rùn)滑力采用Antal模型，湍流耗散力采用Burns模型。由于升力和虛擬質(zhì)量力對(duì)該問(wèn)題的計(jì)算結(jié)果影響不大，且添加升力后會(huì)使程序計(jì)算的收斂性變差，因此本文在計(jì)算中不添加升力模型。液相湍流模型采用標(biāo)準(zhǔn)k-ε模型。梯度離散采用Gauss linear格式，其他采用一階迎風(fēng)差分格式。計(jì)算殘差要求小于10-4。

2.3 網(wǎng)格無(wú)關(guān)性分析

過(guò)冷沸騰由于其局部加熱的特性，近壁面參數(shù)與網(wǎng)格相關(guān)，因此不能采用近壁面參數(shù)作為網(wǎng)格無(wú)關(guān)性判斷的依據(jù)。本文選取沿軸向長(zhǎng)度變化的截面平均空泡份額作為網(wǎng)格無(wú)關(guān)性的考核依據(jù)，計(jì)算結(jié)果如圖3所示。由圖3可知，軸向控制體網(wǎng)格數(shù)對(duì)計(jì)算結(jié)果的影響不明顯，而當(dāng)徑向控制體網(wǎng)格數(shù)為30與35時(shí)，計(jì)算結(jié)果相差不大。所以，經(jīng)過(guò)網(wǎng)格無(wú)關(guān)性計(jì)算，在徑向及軸向分別取30和2 000個(gè)節(jié)點(diǎn)，總網(wǎng)格數(shù)為60 000。此時(shí)，網(wǎng)格尺度y+在40～90之間，滿足兩相流計(jì)算要求(20

圖3 網(wǎng)格無(wú)關(guān)性分析Fig.3 Grid independence analysis

2.4 計(jì)算結(jié)果分析

截面平均空泡份額、液相平均溫度和壁面溫度的數(shù)值模擬結(jié)果與實(shí)驗(yàn)值的對(duì)比示于圖4?？张莘蓊~與反應(yīng)堆的穩(wěn)定性及堆芯中子動(dòng)力學(xué)密切相關(guān)，因此空泡份額的準(zhǔn)確預(yù)測(cè)在核反應(yīng)堆的設(shè)計(jì)及運(yùn)行中非常重要。由圖4可知，數(shù)值模擬結(jié)果與實(shí)驗(yàn)值符合良好，說(shuō)明OpenFOAM中嵌入的歐拉兩流體模型及輔助模型可很好地預(yù)測(cè)該實(shí)驗(yàn)條件下的流動(dòng)沸騰現(xiàn)象。

圖4 計(jì)算值與實(shí)驗(yàn)值的對(duì)比Fig.4 Comparison between calculation results and experimental data

圖5 計(jì)算結(jié)果云圖Fig.5 Cloud map of simulation result

計(jì)算的液相流速、液相溫度及空泡份額示于圖5。可看出，在流道出口處氣相占比較大，擠占了大部分流通面積，致使液相流速提高；加熱管道中，在0.9 m處開(kāi)始產(chǎn)生氣泡，1.3～2.0 m處橫截面上空泡份額的峰值并不是出現(xiàn)在壁面上，而是在距離加熱壁面約R/5位置處，這是氣液兩相間作用力導(dǎo)致的，其中，壁面潤(rùn)滑力推動(dòng)氣泡向流道中心運(yùn)動(dòng)，湍流耗散力使氣泡沿空泡份額梯度方向運(yùn)動(dòng)，故二者均有展平空泡份額在徑向上分布的效果。

3 結(jié)論

本文基于開(kāi)源CFD平臺(tái)OpenFOAM，采用歐拉兩流體模型對(duì)4.5 MPa下豎直圓管內(nèi)的過(guò)冷流動(dòng)沸騰現(xiàn)象進(jìn)行了數(shù)值模擬。計(jì)算中考慮了氣相與液相間的曳力、壁面潤(rùn)滑力和湍流耗散力，熱流密度分配模型采用了考慮氣相傳熱的四分模型，氣泡脫離直徑、核化密度和氣泡脫離頻率模型分別采用Tolubinsky模型、Lemmert-Chawla模型和Cole模型。計(jì)算得到了空泡份額、液相平均溫度和壁面溫度沿圓管軸向的分布，計(jì)算值與實(shí)驗(yàn)值符合良好，證明上述模型對(duì)該工況下過(guò)冷流動(dòng)沸騰現(xiàn)象有較好的模擬能力。