基于OpenFOAM的過冷流動沸騰數(shù)值模擬

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

(西安交通大學(xué) 陜西省先進(jìn)核能技術(shù)重點實驗室，動力工程多相流國家重點實驗室，陜西 西安 710049)

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

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

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

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

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

1.1 相間動量交換模型

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

Flg=-Fgl

(1)

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

M=Md+Ml+Mwl+Mtd+Mvm

(2)

1.2 壁面熱流密度分配模型

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

Qtot=Qc+Qq+Qe

(3)

但上述表達(dá)式未考慮氣泡覆蓋壁面時壁面與氣相之間的單相對流換熱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)驗關(guān)系式，本文采用Lavieville等[10]提出的模型：

(5)

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

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

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

(6)

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

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

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

Qq=Awhq(Tw-Tl)

(7)

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

(8)

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

twait=0.8/f

(9)

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

Qe=mwhlg

(10)

式中：hlg為汽化潛熱；mw為汽化的傳質(zhì)速率，其經(jīng)驗關(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計算關(guān)系式為：

K=4.8e-Ja/80

(13)

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

(14)

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

1.3 輔助模型

1) 氣泡脫離直徑模型

氣泡從加熱壁面上脫離時的直徑稱為氣泡脫離直徑，其大小取決于氣泡長大過程中的受力情況。本文采用Tolubinsky等[13]根據(jù)實驗觀測得出的經(jīng)驗?zāi)Ｐ陀嬎銡馀菝撾x直徑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]計算，該模型認(rèn)為核化密度N與壁面過熱度ΔTsup有關(guān)：

N=(CNΔTsup)n

(16)

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

3) 氣泡脫離頻率模型

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

(17)

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

2.1 算例描述

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

圖2 實驗段示意圖Fig.2 Schematic diagram of experimental section

2.2 計算方法

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

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

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

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

2.4 計算結(jié)果分析

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

圖4 計算值與實驗值的對比Fig.4 Comparison between calculation results and experimental data

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

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

3 結(jié)論

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