尚澤敏，楊立新,*，袁小菲，劉 偉，劉 余，蔣汀嵐

(1.北京交通大學(xué) 機(jī)械與電子控制工程學(xué)院，北京 100044； 2.中國(guó)核動(dòng)力研究設(shè)計(jì)院 核反應(yīng)堆系統(tǒng)設(shè)計(jì)技術(shù)重點(diǎn)實(shí)驗(yàn)室，四川 成都 610213)

采用流動(dòng)換熱模式是冷卻高功率密度反應(yīng)堆燃料元件的有效方法。與單相強(qiáng)制對(duì)流相比，過(guò)冷流動(dòng)沸騰能大幅提高換熱效率，在高熱流密度的場(chǎng)景中得到越來(lái)越多的應(yīng)用。研究過(guò)冷流動(dòng)沸騰的流動(dòng)特性和傳熱特性對(duì)反應(yīng)堆的運(yùn)行安全及經(jīng)濟(jì)性具有重要意義。

隨著計(jì)算機(jī)硬件的飛速發(fā)展和商業(yè)軟件的不斷進(jìn)步，利用計(jì)算流體力學(xué)(CFD)方法模擬流動(dòng)沸騰現(xiàn)象引起了越來(lái)越多的關(guān)注。流動(dòng)沸騰的數(shù)值模擬需要多種模型相互耦合，Judd等[1]在實(shí)驗(yàn)觀察和理論分析的基礎(chǔ)上，提出了一種反映池沸騰過(guò)程的機(jī)理模型，將壁面總熱流分為液相對(duì)流傳熱、蒸發(fā)和淬冷熱流3部分。隨后，Kurul等[2]將該模型編譯成CFD代碼來(lái)模擬垂直管道內(nèi)的過(guò)冷沸騰。目前流動(dòng)沸騰數(shù)值計(jì)算方法大多基于Kurul等給出的壁面沸騰模型，又稱(chēng)為PRI模型。秦浩等[3]基于OpenFOAM平臺(tái)對(duì)管內(nèi)流動(dòng)沸騰現(xiàn)象進(jìn)行了數(shù)值模擬，驗(yàn)證了模型的有效性。Filho等[4]基于FLUENT對(duì)歐拉雙流體模型中氣液界面換熱系數(shù)及氣泡脫離頻率進(jìn)行了研究。李權(quán)等[5]基于歐拉雙流體模型對(duì)非均勻加熱圓管的臨界熱流密度(CHF)及臨界位置進(jìn)行了預(yù)測(cè)。鄭樂(lè)樂(lè)等[6]使用STAR-CCM+基于歐拉雙流體模型對(duì)管內(nèi)核態(tài)沸騰進(jìn)行了數(shù)值模擬，并對(duì)DNB型臨界熱流密度進(jìn)行了預(yù)測(cè)。上述研究已成功模擬了多種條件下的流動(dòng)沸騰現(xiàn)象，CFD方法可靠性得到越來(lái)越多的驗(yàn)證，同時(shí)其中相關(guān)模型也得到了諸多改進(jìn)。但目前針對(duì)歐拉雙流體模型框架內(nèi)眾多的輔助模型進(jìn)行參數(shù)敏感性分析的研究較少，且對(duì)各種相間作用力的計(jì)算與否尚無(wú)統(tǒng)一認(rèn)識(shí)。

因此，本研究基于商用CFD分析軟件STAR-CCM+，采用歐拉雙流體模型結(jié)合 RPI 模型對(duì)管內(nèi)過(guò)冷流動(dòng)沸騰進(jìn)行數(shù)值模擬?；趯?shí)驗(yàn)結(jié)果對(duì)兩相計(jì)算中網(wǎng)格模型、湍流模型、沸騰子模型及相間作用力模型的參數(shù)進(jìn)行敏感性分析。

1 模型描述

1.1 控制方程

歐拉雙流體模型將兩相流場(chǎng)中各相均視作連續(xù)介質(zhì)且充滿整個(gè)流場(chǎng)，各相流動(dòng)變量在交界面上發(fā)生間斷并產(chǎn)生質(zhì)量、動(dòng)量和能量傳遞現(xiàn)象。由于流動(dòng)沸騰過(guò)程中氣泡與液相之間存在較強(qiáng)的相互作用力，如曳力和非曳力等，因此更適合采用兩流體模型描述。

體積守恒方程(流場(chǎng)中氣體與液相體積分?jǐn)?shù)(αl、αg)之和為1)：

αl+αg=1

(1)

連續(xù)性方程：

(2)

動(dòng)量方程：

(3)

能量方程：

(4)

1.2 相間質(zhì)量、能量及動(dòng)量傳遞

氣液相間的質(zhì)量傳遞速率取決于從液體側(cè)和蒸汽側(cè)到相間飽和界面的傳熱。假設(shè)兩相間的能量傳遞均通過(guò)交界面進(jìn)行，則有：

Qlg=hlga(Tsat-Tl)

(5)

Qgl=hgla(Tsat-Tg)

(6)

式中：Qlg和Qgl分別為氣液界面與液相和氣相間的傳熱；a為相界面面積密度；hlg和hgl分別為氣液界面與液相和氣相間的換熱系數(shù)，通常用下式計(jì)算：

hlg=Nulλl/db

(7)

hgl=Nugλg/db

(8)

其中：λl和λg分別為液相和氣相的熱導(dǎo)率；db為氣泡直徑。對(duì)于液相，Nul通常根據(jù)Ranz-Marshall關(guān)系式[7]計(jì)算。

(9)

而氣相的努塞爾數(shù)Nug一般給定為常數(shù)。

(10)

式中，hfg為汽化潛熱。

如上所述，流動(dòng)沸騰過(guò)程中氣泡與液相之間存在較強(qiáng)的相互作用力Mk，習(xí)慣上將其分為曳力FD和非曳力，其中非曳力包括升力FL、湍流耗散力FTD、壁面潤(rùn)滑力FWL、虛擬質(zhì)量力FVM，則有：

Mk=FD+FTD+FL+FWL+FVM

(11)

這5種相間作用力的計(jì)算公式及常用模型列于表1，并簡(jiǎn)單描述了各作用力對(duì)氣泡的影響。

1.3 壁面沸騰模型

流體流動(dòng)過(guò)程中受到加熱，壁面總熱量qwall可分為液相對(duì)流傳熱ql、蒸發(fā)傳熱qq、氣液置換(淬冷)傳熱qe和氣相對(duì)流傳熱qg4部分。壁面沸騰模型公式如下：

qwall=(ql+qq+qe)(1-Kdry)+Kdryqg

(12)

式中，Kdry為蒸汽壁面接觸面積分?jǐn)?shù)，采用下式計(jì)算：

(13)

(14)

其中：αδ為近壁面網(wǎng)格單元的平均空泡份額；αdry為臨界空泡份額，默認(rèn)值為0.9。4部分熱量計(jì)算公式可參考文獻(xiàn)[13]。其中，淬冷傳熱量qe計(jì)算公式中的Kquench為氣泡影響壁面面積分?jǐn)?shù)，可根據(jù)Kurul等[2]提出的模型計(jì)算：

(15)

式中，F(xiàn)A為在汽化核心密度和氣泡影響面積分?jǐn)?shù)間進(jìn)行比例縮放的系數(shù)。

蒸發(fā)熱流量計(jì)算公式中的Dd、Nw和f分別為氣泡脫離直徑、汽化核心密度和氣泡脫離頻率。這些參數(shù)需要額外的輔助模型才能完成對(duì)方程的封閉，輔助模型來(lái)源于實(shí)驗(yàn)或由理論推導(dǎo)總結(jié)出的經(jīng)驗(yàn)關(guān)系式。眾多汽化核心密度模型中，形式較為簡(jiǎn)單的Lemmert-Chawla模型[14]得到廣泛應(yīng)用。該模型認(rèn)為汽化核心密度是壁面過(guò)熱度ΔTsup的函數(shù)，其在兩相計(jì)算中的魯棒性較好，容易獲得初始解。然而，此模型基于低壓實(shí)驗(yàn)數(shù)據(jù)整理而來(lái)，實(shí)際上壓力對(duì)汽化核心密度有顯著影響，一般壓力越高，氣泡尺寸越小，更難脫離加熱面，從而使得汽化核心密度增大。為解決Lemmert-Chawla模型在高壓工況下適應(yīng)性不足這一問(wèn)題，在Hibiki-Ishii汽化核心密度模型的基礎(chǔ)上，Li等[15]基于公開(kāi)獲取的實(shí)驗(yàn)數(shù)據(jù)，開(kāi)發(fā)了新的汽化核心密度模型，此模型考慮了壓力、壁面過(guò)熱度及壁面接觸角的影響。本研究使用的輔助模型列于表2。

表1 相間相互作用力模型Table 1 Model of interfacial momentum transfer

2 計(jì)算工況及幾何建模

Bartolomei等[20]對(duì)垂直上升管內(nèi)的過(guò)冷流動(dòng)沸騰進(jìn)行了實(shí)驗(yàn)研究，對(duì)不同管徑、熱流密度、進(jìn)口質(zhì)量流速、進(jìn)口過(guò)冷度及壓力下的軸向空泡份額及溫度分布進(jìn)行了測(cè)量。本文選用壓力4.5 MPa、熱流密度570 kW/m2、進(jìn)口質(zhì)量流速900 kg/(m2·s)、進(jìn)口過(guò)冷度60 K的工況進(jìn)行數(shù)值模型研究。實(shí)驗(yàn)中加熱段長(zhǎng)2 m、管道直徑15.4 mm。

為減小進(jìn)口段效應(yīng)的影響，幾何建模時(shí)在進(jìn)口與加熱段之間留出200 mm的流動(dòng)發(fā)展空間(大于10倍當(dāng)量直徑)。在加熱段與出口之間設(shè)置100 mm長(zhǎng)的絕熱段以保持?jǐn)?shù)值計(jì)算的穩(wěn)定性。采用二維軸對(duì)稱(chēng)邊界條件，使用結(jié)構(gòu)化網(wǎng)格對(duì)流體域進(jìn)行劃分。建模的幾何結(jié)構(gòu)如圖1所示。

圖1 實(shí)驗(yàn)段幾何結(jié)構(gòu)示意圖Fig.1 Geometry of experimental section

3 模型敏感性分析

3.1 網(wǎng)格及湍流模型的敏感性分析

為研究網(wǎng)格對(duì)計(jì)算結(jié)果的影響，本文從徑向和軸向兩個(gè)維度進(jìn)行分析。設(shè)置不同的徑向分割數(shù)使加熱面處第1層網(wǎng)格高度從0.003 mm到0.4 mm變化。當(dāng)前工況下，對(duì)應(yīng)壁面y+范圍為1～150，網(wǎng)格參數(shù)設(shè)置列于表3。隨后，在徑向20層均勻分割的條件下，設(shè)置不同的軸向拉伸層數(shù)(85、115、145及175)生成另外4組網(wǎng)格。

選取標(biāo)準(zhǔn)k-ω模型、SSTk-ω模型、標(biāo)準(zhǔn)k-ε模型、Realizablek-ε模型及Realizablek-ε兩層模型共5種湍流模型進(jìn)行求解。不同湍流模型需配合相應(yīng)的壁面處理方法，STAR-CCM+中提供了低y+、高y+及全y+3種壁面處理方法。根據(jù)網(wǎng)格和湍流模型的不同組合并匹配相應(yīng)的壁面處理類(lèi)型得到表4所列的算例。其中，標(biāo)準(zhǔn)k-ε模型和Realizablek-ε模型均為高雷諾數(shù)模型，僅可使用高y+壁面處理方法。Realizablek-ε兩層模型僅可使用兩層全y+壁面處理方法。采用Realizablek-ε兩層模型和SSTk-ω兩層模型求解的不同徑向網(wǎng)格分布下沿程空泡份額示于圖2。由圖2可知，徑向網(wǎng)格層數(shù)對(duì)計(jì)算結(jié)果影響較大。

表3 徑向網(wǎng)格分布Table 3 Radial grid distribution

結(jié)合圖2a、b可發(fā)現(xiàn)，隨著徑向網(wǎng)格層數(shù)的增加，空泡份額計(jì)算值與實(shí)驗(yàn)值的偏差逐漸增大。Mesh1和Mesh2網(wǎng)格數(shù)量最多，但在采用Realizablek-ε兩層模型的計(jì)算中未能獲得收斂解。最精細(xì)的網(wǎng)格計(jì)算偏差反而最大，可見(jiàn)對(duì)于采用歐拉雙流體模型進(jìn)行的兩相計(jì)算，并非計(jì)算域中劃分的網(wǎng)格越多越準(zhǔn)確。對(duì)流動(dòng)沸騰現(xiàn)象的準(zhǔn)確模擬很大程度上依賴(lài)于近壁區(qū)域傳熱傳質(zhì)現(xiàn)象的建模，加熱面處第1層的網(wǎng)格高度過(guò)低(最小y+<40)會(huì)導(dǎo)致流動(dòng)沸騰計(jì)算中的局部參數(shù)，尤其是近壁區(qū)空泡份額發(fā)生劇烈變化，進(jìn)而影響計(jì)算的穩(wěn)定性，調(diào)整松弛因子依然難以獲得收斂解或計(jì)算結(jié)果出現(xiàn)較大的偏差。較高的第1層網(wǎng)格可保證氣泡的瞬態(tài)行為(如產(chǎn)生、生長(zhǎng)和脫離)被網(wǎng)格覆蓋，相關(guān)參數(shù)被平均化處理，進(jìn)而保證計(jì)算過(guò)程的穩(wěn)定性。

表4 湍流模型、網(wǎng)格及壁面處理算例Table 4 Computational cases of turbulence model, grid and wall function

圖2 不同湍流模型求解的各徑向網(wǎng)格空泡份額的變化Fig.2 Void fraction under each grid solved by different turbulence models

保持相同的徑向網(wǎng)格分布，不同軸向網(wǎng)格拉伸層數(shù)對(duì)空泡份額分布的影響示于圖3。由圖3可知，對(duì)當(dāng)前具有較高長(zhǎng)徑比(L/D>150)的流域，軸向網(wǎng)格層數(shù)對(duì)空泡份額分布的影響很小，因此后續(xù)計(jì)算設(shè)置軸向拉伸層數(shù)為115，既能保證計(jì)算精度，又可減少網(wǎng)格量。

圖3 軸向網(wǎng)格拉伸層數(shù)對(duì)空泡份額的影響Fig.3 Void fraction under different axial stretch layers

Mesh5網(wǎng)格尺寸下，采用不同湍流模型進(jìn)行求解所得的空泡份額及壁面溫度分布如圖4所示。由圖4可知，各湍流模型對(duì)壁面溫度分布的計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果的差異較小，且不同湍流模型間的計(jì)算結(jié)果差異不大。標(biāo)準(zhǔn)k-ε模型對(duì)軸向空泡份額分布的計(jì)算存在較大偏差，出口處的空泡份額高于其他湍流模型，且對(duì)氣相出現(xiàn)位置的計(jì)算與其他模型存在較大差異。Realizablek-ε模型及Realizablek-ε兩層模型能很好地反映軸向空泡份額分布。

3.2 沸騰模型敏感性分析

為進(jìn)行沸騰子模型的敏感性分析，將Lemmert-Chawla(LC)、Hibiki-Ishii(HI)和李權(quán)(Li)等開(kāi)發(fā)的汽化核心密度模型、Tolubinsky-Kostanchuk(TK)和Kocamustafaogullari(KI)氣泡脫離直徑模型以及3個(gè)氣泡影響面積縮放系數(shù)(1、2和4)進(jìn)行組合得到表5所列算例。

圖4 Mesh5網(wǎng)格下不同湍流模型的計(jì)算結(jié)果Fig.4 Calculation results of different turbulence models under Mesh5

首先對(duì)不同的氣泡影響面積縮放系數(shù)FA進(jìn)行敏感性分析，隨后分析汽化核心密度和氣泡脫離直徑模型的影響。

表5 壁面沸騰模型計(jì)算工況Table 5 Computational cases of different boiling models

不同模型計(jì)算結(jié)果的對(duì)比如圖5所示，其中圖5a為3組不同的FA對(duì)計(jì)算結(jié)果的影響，圖5b為不同沸騰子模型組合的計(jì)算結(jié)果對(duì)比。由圖5a可知，F(xiàn)A對(duì)參數(shù)分布的影響很小，計(jì)算中可取FA為2，該值由Bartolomei等[20]根據(jù)實(shí)驗(yàn)測(cè)得。由圖5b可知，除LC-KI模型組合外，其余5組模型均能對(duì)壁面溫度分布做出較好的預(yù)測(cè)。HI-KI模型組合和Li-KI模型組合對(duì)壁面溫度的計(jì)算結(jié)果較為接近，LC-TK模型計(jì)算的壁面溫度與實(shí)驗(yàn)結(jié)果最接近。此外可看出，采用基于局部過(guò)冷度得到氣泡脫離直徑的TK模型求解的空泡份額與實(shí)驗(yàn)值存在較大差異，過(guò)早地預(yù)測(cè)了氣泡產(chǎn)生的位置，但隨著空泡份額的增大，偏差逐漸減小。不同沸騰子模型組合對(duì)主流溫度的計(jì)算結(jié)果在加熱段的前半段差異較小，靠近出口段存在一定偏差。結(jié)合壁面溫度及空泡份額分布曲線可得，兩相計(jì)算中推薦采用HI-KI或Li-KI沸騰子模型組合。

圖5 不同沸騰模型計(jì)算結(jié)果對(duì)比Fig.5 Comparison of different boiling models

3.3 相間作用力模型敏感性分析

為評(píng)估相間作用力對(duì)計(jì)算結(jié)果的影響，將其分為曳力及非曳力兩部分進(jìn)行對(duì)比。在同一曳力模型(Tomiyama模型)下，通過(guò)不同非曳力模型組合得到表6所列算例。湍流耗散力選用Burns模型，壁面潤(rùn)滑力選用Antal模型，虛擬質(zhì)量力采用Auton模型，升力模型為T(mén)omiyama模型。隨后，再對(duì)比7種曳力模型對(duì)計(jì)算結(jié)果的影響。

不同非曳力模型組合下的軸向參數(shù)分布對(duì)比示于圖6。由圖6可知，非曳力中的湍流耗散力對(duì)空泡份額分布及主流溫度分布有較大影響，計(jì)算中若不考慮湍流耗散力，則會(huì)導(dǎo)致空泡份額升高，原因是湍流耗散力促使氣泡由加熱面向過(guò)冷的主流中擴(kuò)散，使得更多氣泡冷凝，從而降低空泡份額。而升力、虛擬質(zhì)量力、壁面潤(rùn)滑力對(duì)軸向參數(shù)分布的影響較小。不同模型組合均能對(duì)壁面溫度分布做出較好的預(yù)測(cè)。

不同非曳力模型對(duì)徑向空泡份額及軸向壁面溫度分布的影響示于圖7。由圖7可知，不考慮非曳力時(shí)，氣泡會(huì)集中分布于靠近加熱面處，從而引起壁面溫度飛升。增加湍流耗散力后，徑向空泡份額呈現(xiàn)展平趨勢(shì)。增加升力會(huì)使空泡份額的峰值位置由壁面向主流轉(zhuǎn)移。而虛擬質(zhì)量力、壁面潤(rùn)滑力對(duì)徑向空泡份額分布影響很小，可忽略不計(jì)。根據(jù)以上分析可知，采用歐拉雙流體模型進(jìn)行兩相計(jì)算時(shí)，加入湍流耗散力及升力的影響即可獲得較好的結(jié)果。

表6 相間相互作用力計(jì)算工況Table 6 Computational cases of different interaction forces

不同曳力模型的計(jì)算結(jié)果示于圖8。由圖8可知，對(duì)于當(dāng)前工況，各曳力模型均能較好地反映出軸向空泡份額的分布，其中TY曳力模型和對(duì)稱(chēng)曳力模型的計(jì)算結(jié)果與實(shí)驗(yàn)值吻合最好。

圖6 不同非曳力模型組合下的軸向參數(shù)分布Fig.6 Axial parameter distribution under different non-drag force model combinations

圖7 非曳力模型對(duì)徑向空泡份額及軸向壁面溫度分布的影響Fig.7 Influence of non-drag force model on radial void fraction and axial wall temperature distribution

4 結(jié)論

本研究基于歐拉雙流體模型，結(jié)合壁面沸騰模型對(duì)管內(nèi)過(guò)冷流動(dòng)沸騰進(jìn)行了數(shù)值模擬?；趯?shí)驗(yàn)結(jié)果對(duì)模型參數(shù)進(jìn)行了敏感性分析。通過(guò)對(duì)比分析識(shí)別出了對(duì)結(jié)果影響較大的參數(shù)并給出了推薦設(shè)置。

圖8 曳力模型對(duì)計(jì)算結(jié)果的影響Fig.8 Influence of drag force model on calculation result

采用歐拉雙流體模型進(jìn)行過(guò)冷流動(dòng)沸騰計(jì)算時(shí)，需注意加熱面y+不能太低，過(guò)低的y+值會(huì)導(dǎo)致收斂困難或產(chǎn)生較大偏差，本研究中加熱面y+處于40～150范圍內(nèi)可獲得較好的計(jì)算結(jié)果。對(duì)其他工況或幾何結(jié)構(gòu)開(kāi)展兩相CFD計(jì)算時(shí)可采用文中網(wǎng)格敏感性分析的方法來(lái)調(diào)整加熱面處網(wǎng)格的高度。湍流模型推薦選用Realizablek-ε模型或Realizablek-ε兩層模型。壁面沸騰模型中的汽化核心密度和氣泡脫離直徑模型推薦選用HI-KI或Li-KI模型組合。各相間作用力中對(duì)計(jì)算影響較大的有湍流耗散力及升力，而虛擬質(zhì)量力及壁面潤(rùn)滑力的影響可忽略。