欠熱沸騰模型參數(shù)敏感性分析

郜 冶，李小暢，劉平安

（哈爾濱工程大學(xué) 航天與建筑工程學(xué)院，黑龍江 哈爾濱 150001）

近十年來，基于歐拉-歐拉兩流體六方程模型的RPI（Rensselaer Polytechnic Institute）欠熱沸騰數(shù)值模型［1］得到了較大的發(fā)展，越來越多地應(yīng)用于反應(yīng)堆工程熱工水力研究［2-3］。然而，由于沸騰機(jī)理尚不十分明確，現(xiàn)有的基于RPI模型的欠熱沸騰數(shù)值計(jì)算方法主要建立在眾多經(jīng)驗(yàn)或半經(jīng)驗(yàn)封閉子模型的基礎(chǔ)上，這些子模型適用范圍有限，且模型參數(shù)無通用確定方法，因此，對工況的依賴性較大。封閉子模型中的氣泡脫離壁面直徑、氣泡脫離頻率及成核面密度模型關(guān)系到氣泡份額和壁面溫度等重要熱工參數(shù)計(jì)算結(jié)果的準(zhǔn)確性。許多學(xué)者針對不同實(shí)驗(yàn)工況及流動工質(zhì)提出了不同的封閉子模型，但均缺乏普適性［4］。另外，氣泡脫離壁面后的運(yùn)動主要受氣液兩相間相互作用力影響，氣泡的運(yùn)動反過來又會影響近壁面氣泡份額和換熱效率，因此兩相間相互作用力模型同樣十分重要?？梢娀诂F(xiàn)有模型研究各封閉子模型及相間作用力模型的具體影響方式對正確使用RPI模型進(jìn)行數(shù)值計(jì)算意義重大。

本文基于Nylund等［5］的FT-6a實(shí)驗(yàn)，研究3個重要封閉子模型（氣泡脫離壁面直徑、氣泡脫離頻率及氣泡成核面密度）及兩個重要相間非曳力模型（升力和湍流耗散力）對氣泡軸向及徑向分布和壁面過熱度等計(jì)算結(jié)果的影響方式。

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

RPI欠熱沸騰數(shù)值計(jì)算方法基于歐拉-歐拉兩流體六方程模型、壁面熱流量分區(qū)模型及相間傳輸模型。本文主要給出壁面熱流量分區(qū)模型及相間傳輸模型，歐拉-歐拉兩流體六方程模型作為最基本的守恒方程本文不再給出。

1.1 壁面熱流量分區(qū)模型

Kurul等［1］提出的RPI模型的核心機(jī)理是基于Judd等［6］的實(shí)驗(yàn)和理論研究：通道內(nèi)發(fā)生壁面沸騰時可將壁面熱流量分為單相對流傳熱、蒸發(fā)傳熱及氣液置換傳熱3部分。壁面熱流量分區(qū)模型基本方程如下：式中：qW、qC、qQ和qE分別為壁面總熱流量、單相對流傳熱量、氣液置換傳熱量和蒸發(fā)傳熱量；hC為單相對流換熱系數(shù)，通過壁面函數(shù)求解；Tw、Tl分別為壁面溫度與近壁面液相溫度；kl、ρl、cpl分別為液相導(dǎo)熱系數(shù)、密度與比定壓熱容；tw、f 分 別 為 氣 泡 脫 離 周 期 與 頻 率；ρv 為 氣相密度；F 為氣泡對壁面的影響因子，取值為2；Nw為壁面成核面密度；Dw為氣泡脫離壁面直徑；Av、Al分別為壁面上氣相與液相的單位覆蓋面積，Av=min（1，πF2DwNw／4），Al=1-Av；hfv為汽化潛熱。由于該模型在計(jì)算中將氣相溫度固定為系統(tǒng)壓力下的飽和溫度，因此僅適用于氣泡體積分?jǐn)?shù)較低（αmax＜0.8）的情況。

式（1）中氣泡脫離壁面直徑、氣泡脫離頻率及壁面成核面密度未知，需要相應(yīng)子模型對其進(jìn)行封閉，本文選用如下使用相對較為廣泛的模型［7］：

式中：Tsub、Tsup分別為液相欠熱度和壁面過熱度，Tsup=Tw-Tsat，Tsub=Tsat-Tl；Tref為參考溫度；Dref為氣泡脫離壁面參考直徑；Nref為壁面成核參考面密度；Cf為氣泡脫離阻力因子；Dmax為氣泡脫離壁面時的最大直徑。式（2）中多處涉及近壁面液相溫度，為減小對近壁面網(wǎng)格的依賴，本文取y+=250時第1層網(wǎng)格節(jié)點(diǎn)上的液相溫度。

1.2 相間傳輸模型

1）質(zhì)量傳輸

氣液兩相間的質(zhì)量傳輸包括過熱壁面上的汽化和主流區(qū)的蒸發(fā)或冷凝。過熱壁面上液相的汽化速率可由式（1）中的蒸發(fā)傳熱量得到：

氣泡脫離壁面后在主流區(qū)的蒸發(fā)和冷凝速率如下：

式中：Гlv、Гvl分別為蒸發(fā)和冷凝速率；Ai為相界面積，根據(jù)氣泡直徑Db求解；hi為相界面換熱系數(shù)，根據(jù)Ranz-Marshall關(guān)系式獲得［8］。

2）能量傳輸

氣液兩相間的能量傳輸以相界面作為參考，兩相間的熱量傳輸為：

式中，qlv與qvl分別為液相向氣相傳熱量及氣相向液相傳熱量，二者大小相同，正負(fù)相反。

3）動量傳輸

兩相間動量傳輸用相互作用力表示，流動沸騰中總的相互作用力FT包括曳力FD和非曳力FND兩類，其中非曳力又包含升力FL、湍流耗散力FTP及壁面潤滑力FWL：

各相互作用力的求解模型為：

式中：CD為曳力系數(shù)，由Ishii-Zuber關(guān)系式確定［9］；v為速度矢 量；Δ為哈密頓算子；Db為氣泡直徑；CL為升力系數(shù)；yw為壁面法向距離；Cw1、Cw2為 無 量 綱 參 數(shù)，分 別 取-0.01 和0.05［10］；FTD，F(xiàn)avre、FTD，Lopze分別為湍流耗散力 的Favre平 均 模 型［11］和Lopze模 型［12］，其 中CD與曳力模型相同，CTD為湍流耗散系數(shù)；νt、σt分別為液相湍流黏性系數(shù)及湍流Schmidt數(shù)；kl為液相湍流動能。

4）氣泡直徑及相界面積

本文采用Kurul等［1］推薦的基于當(dāng)?shù)匾合嗲窡岫葰馀葜睆侥Ｐ停?/p>

基于該氣泡直徑模型，相界面積則可采用如下對稱模型獲得：

2 數(shù)值計(jì)算方法

數(shù)值計(jì)算基于ANSYS CFX 14.5軟件包，該程序提供了歐拉-歐拉兩流體六方程模型框架，在加熱壁面邊界條件的處理中嵌入了RPI模型，并提供了部分相間傳輸模型，未提供的模型可利用User Fortran以源項(xiàng)的形式加入到守恒方程中。

本文以Nylund等［5］的FT-6a實(shí)驗(yàn)為研究對象，實(shí)驗(yàn)?zāi)Ｐ图俺叽鐓?shù)如圖1a所示。實(shí)驗(yàn)以水為工質(zhì)，系統(tǒng)壓力為5.0MPa，入口欠熱度為4.5K，入口流量為1 163kg／（m2·K），均勻壁面熱流量為522kW／m2。取實(shí)驗(yàn)?zāi)Ｐ椭芟虻?／10和軸向加熱段前2m 作為數(shù)值計(jì)算對象，如圖1a中陰影部分所示。氣相物性參數(shù)取系統(tǒng)壓力下的飽和蒸汽參數(shù)；液相則取進(jìn)出口的平均值。液相湍流模擬采用k-ω SST 模型［13］，氣相則采用零方程模型。計(jì)算中建立了3套網(wǎng)格：286×50、391×100、605×200，各網(wǎng)格計(jì)算的平均氣泡份額差值小于3%，選取第2套網(wǎng)格作為計(jì)算對象，如圖1b所示。

圖1 FT-6a幾何模型及計(jì)算網(wǎng)格Fig.1 Sketch of FT-6ageometry and computational mesh

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

分析RPI模型中3 個子模型（Dw、Nw及f）及兩個相間非曳力模型（FL及FTD）對計(jì)算結(jié)果的影響。從式（1）及（2）可看出，RPI模型主要由Dw、Nw及f 確定，對于一具體的流動沸騰，這些封閉子模型主要由式（2）中的Dref、Cf及Nref決定。從式（7）則可看出，升力主要受升力系數(shù)CL影響，F(xiàn)avre湍流耗散力模型參數(shù)CD已由曳力模型確定，Lopze湍流耗散力模型則主要受耗散系數(shù)CTD影響。因此，針對上述所選定的RPI封閉子模型及相間非曳力模型的研究可通過分析如下5個參數(shù)對數(shù)值計(jì)算結(jié)果的 敏 感 性 來 實(shí) 現(xiàn)：Dref、Cf、Nref、CL及CTD。根據(jù)文獻(xiàn)［7］，5個參數(shù)推薦值列于表1。

表1 模型參數(shù)參考值Table 1 Reference values of model parameters

FT-6a實(shí)驗(yàn)測量了氣泡沿軸向和徑向的分布，未給出壁面過熱度，本文采用Jens 關(guān)聯(lián)式［15］估算充分發(fā)展的欠熱流動沸騰壁面過熱度：

式中：q 為 壁 面 熱 流 量，MW／m2；p 為 系 統(tǒng) 壓力，Pa。基于式（10）所計(jì)算的FT-6a實(shí)驗(yàn)壁面過熱度為9.49K。另外，Bartolomej等［16］基于圓管的類似工況（4.5 MPa，0.57 MW／m2）下的欠熱沸騰實(shí)驗(yàn)結(jié)果約為10K，因此可基本確定Jens關(guān)聯(lián)式計(jì)算結(jié)果的準(zhǔn)確性。

3.1 子模型參數(shù)敏感性分析

計(jì)算時升力及湍流耗散力分別采用Tomiyama和Favre模型，其他子模型參數(shù)取表1中的值。圖2為Dref對軸向平均氣泡份額及壁面過熱度的影響。從圖2可看出，隨參考直徑的增大，平均氣泡份額逐漸增大，壁面過熱度則隨之降低。這是由于隨壁面氣泡直徑的增大，蒸發(fā)傳熱量及氣液置換傳熱量隨之增大，由蒸發(fā)傳熱產(chǎn)生的氣泡份額增加，此時傳熱效率升高，壁面溫度降低。參考直徑達(dá)1.3mm 時，平均氣泡份額與實(shí)驗(yàn)結(jié)果符合較好，但壁面過熱度大幅低于Jens關(guān)聯(lián)式計(jì)算值，參考直徑取0.9mm 時壁面過熱度與關(guān)聯(lián)式符合較好。該圖表明Dw對平均氣泡份額和壁面溫度的計(jì)算結(jié)果影響均較大。

圖3示出Nref對平均氣泡份額及壁面過熱度的影響。圖3表明，隨Nref的增大，氣泡份額逐漸增大，壁面過熱度則隨之減小。顯然這與物理事實(shí)相符，壁面成核點(diǎn)越多，氣泡產(chǎn)生量越多，傳熱效率提高，壁面溫度降低。當(dāng)Nref＞8×105m-2時，計(jì)算結(jié)果不再變化，這是由于式（1）中的qE受max函數(shù)限制。在Nref的有效影響范圍內(nèi)，成核面密度對氣泡份額及壁面溫度影響均較明顯。

圖4為Cf對平均氣泡份額及壁面過熱度的影響。從圖4 可看出，當(dāng)減小阻力因子時（＜1），該參數(shù)對氣泡份額的計(jì)算結(jié)果幾乎無影響，對壁面過熱度的影響也較小。當(dāng)增大阻力因子時（＞1），氣泡份額有所下降，壁面過熱度則大幅升高。

圖2 Dref對軸向平均氣泡份額及壁面過熱度的影響Fig.2 Influence of bubble departure diameter on axial average void fraction and wall superheat

圖3 Nref對平均氣泡份額及壁面過熱度的影響Fig.3 Influence of wall nucleation site density on axial average void fraction and wall superheat

圖4 Cf 對平均氣泡份額及壁面過熱度的影響Fig.4 Influence of bubble detachment frequency on axial average void fraction and wall superheat

圖5 RPI子模型對線平均氣泡份額徑向分布的影響Fig.5 Influence of RPI sub-models on radial distribution of line-averaged void fraction

圖5為RPI各子模型計(jì)算的線平均氣泡份額在1 598mm 高處徑向分布與實(shí)驗(yàn)結(jié)果的對比。從圖5可看出，合適的模型參數(shù)可得到與實(shí)驗(yàn)數(shù)據(jù)基本一致的計(jì)算結(jié)果。圖中各曲線基本平行說明各封閉子模型對氣泡的徑向分布無任何影響，只影響氣泡份額大小。

綜合上述，脫離壁面氣泡直徑及成核密度對氣泡份額及壁面過熱度計(jì)算結(jié)果影響相對較大，氣泡脫離頻率對二者的影響相對較小。不能單獨(dú)對比某一參數(shù)的實(shí)驗(yàn)值來確定計(jì)算結(jié)果的可靠性。對于FT-6a 實(shí)驗(yàn)，推薦參考直徑0.9mm，參考成核面密度8×105m-2，氣泡脫離頻率阻力因子1.0。

3.2 非曳力模型參數(shù)敏感性分析

通過氣泡徑向分布計(jì)算結(jié)果來分析非曳力模型參數(shù)敏感性。圖6為采用不同升力系數(shù)計(jì)算得到的氣泡徑向分布與實(shí)驗(yàn)結(jié)果的對比。從圖6可看出，不同升力系數(shù)主要影響棒束間隙（0.01m）及絕熱壁面附近（0.035m）氣泡份額大小。升力系數(shù)越大，棒束間隙處氣泡份額越小，絕熱壁面附近氣泡份額越大。Tomiyama關(guān)系式與實(shí)驗(yàn)符合較好，CL取0.5 時與Tomiyama關(guān)系式的計(jì)算結(jié)果接近。

圖6 升力系數(shù)敏感性分析Fig.6 Sensitivity analysis of lift coefficient

圖7為不同湍流耗散力模型計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果的對比。從圖7可看出，與升力不同，湍流耗散力主要影響加熱壁面附近氣泡份額大小。從Lopez模型可看出，湍流耗散系數(shù)越大，加熱壁面附近氣泡份額越小。這意味著湍流耗散力有促進(jìn)氣泡遠(yuǎn)離加熱壁面的作用。Favre模型與CTD取0.5時的Lopez模型計(jì)算結(jié)果接近，且與實(shí)驗(yàn)數(shù)據(jù)符合較好。

圖7 湍流耗散力模型敏感性分析Fig.7 Sensitivity analysis of turbulent dispersion force

為進(jìn)一步觀察相間非曳力模型對氣泡徑向分布的影響，圖8示出了1 598 mm 截面處升力及湍流耗散力對氣相分布的影響。從圖8a可看出，CL較?。?.01）時，氣泡更多地向主流區(qū)運(yùn)動；CL較大（1.0）時，氣泡被束縛在壁面附近，Tomiyama模型計(jì)算的結(jié)果則處于二者之間。圖8b表明，CTD較?。?.1）時，氣泡聚集在壁面附近；CTD較大（0.5）時，氣泡則向主流區(qū)擴(kuò)散。分析結(jié)果與圖7基本一致。

圖8 升力（a）及湍流耗散力（b）對氣相分布的影響Fig.8 Influence of lift force（a）and turbulent dispersion force（b）on void fraction contours

4 結(jié)論

1）氣泡脫離壁面直徑模型對氣泡份額及壁面溫度的計(jì)算影響均較大，隨壁面氣泡直徑的增大，氣泡份額增加，壁面溫度降低。

2）成核面密度在其有效影響范圍內(nèi)對壁面過熱度的影響較大，對氣泡份額相對較小。成核面密度越大，氣泡份額越大，壁面過熱度越低。

3）對氣泡脫離頻率而言，當(dāng)氣泡脫離阻力因子在小于1的范圍內(nèi)變化時，對氣泡份額的計(jì)算結(jié)果幾乎無影響，對壁面過熱度的影響也較小；當(dāng)阻力因子大于1時，則對二者的影響均較大。

4）各RPI封閉子模型對氣泡徑向分布無任何影響，只影響氣泡份額大小。

5）不能通過對比單個參數(shù)的實(shí)驗(yàn)測量值來驗(yàn)證計(jì)算的可靠性，應(yīng)綜合對比多個實(shí)驗(yàn)值，以確定最佳模型參數(shù)。

6）相間非曳力模型對氣泡徑向分布的計(jì)算影響較大，升力有阻止氣泡離開加熱壁面的作用，湍流耗散力有促進(jìn)氣泡向主流區(qū)運(yùn)動的作用。

本文可為今后改善欠熱沸騰數(shù)值計(jì)算的準(zhǔn)確性在模型及參數(shù)選取方面提供參考。

［1］ KURUL N，PODOWSKI M Z.On the modeling of multidimensional effects in boiling channels［C］∥27th National Heat Transfer Conference.USA：ANS，1991：28-31.

［2］ LO S，OSMAN J.CFD modeling of boiling flow in PSBT 5×5bundle［J］.Science and Technology of Nuclear Installations，2012（1）：1-8.

［3］ 李小暢，郜冶.壓水堆子通道欠熱沸騰數(shù)值驗(yàn)證及交混翼研究［J］.原子能科學(xué)技術(shù)，2013，47（12）：2 208-2 215.LI Xiaochang，GAO Ye.Numerical validation of subcooled boiling in subchannel of PWR and Research on mixing vane［J］.At Energy Sci Technol，2013，47（12）：2 208-2 215（in Chinese）.

［4］ CHEUNG S C P，VAHAJI S，YEOH G H，et al.Modeling subcooled flow boiling in vertical channels at low pressures，Part 1：Assessment of empirical correlations［J］.International Journal of Heat and Mass Transfer，2014，75（1）：736-753.

［5］ NYLUND K M，BECKET R，EKLUND O，et al.Measurements of hydrodynamics characteris-tics，instability thresholds，and burnout limits for 6-rod clusters in natural and forced circulation［R］.Sweden：FRtGG-I，1967.

［6］ JUDD R L，HWANG K S.A comprehensive model for nucleate pool boiling heat transfer including microlayer evaporation［J］.ASME J Heat Transfer，1976，94（4）：623-629.

［7］ KREPPER E，KONCAR B，EGOROV Y.CFD modeling of subcooled boiling concept，validation and application to fuel assembly design［J］.Nucl Eng Des，2006，237（7）：716-731.

［8］ RANZ W E，MARSHALL W R.Evaporation from drops［J］.Chem Eng Prog，1952，48（4）：173-180.

［9］ ISHII M，ZUBER N.Drag coefficient and relative velocity in bubbly，droplet or particulate flows［J］.AIChE Journal，1979，25（5）：843-855.

［10］ANTAL S P，LAHEY R T，F(xiàn)LAHERTY J E.Analysis of phase distribution in fully developed laminar bubbly two-phase flow［J］.Multiphase Flow，1991，17（5）：635-652.

［11］BURNS A D B，F(xiàn)RANK，T，HAMILL I，et al.The Favre averaged drag model for turbulent dispersion in Eulerian［C］∥5th International Conference on Multiphase Flow.Yokohama，Japan：ICMF，2004：1-17.

［12］LOPEZ D B.Turbulent bubbly flow in a triangular duct［D］.New York：Rensselaer Polytechnic Institute，1991.

［13］KONCAR B，KREPPER E，EGOROV Y.CFD modeling of subcooled flow boiling for nuclear engineering applications［C］∥International Conference of Nuclear Energy for New Europe.Bled Slovenia：NSS，2005：1-14.

［14］TOMIYAMA A，TAMAI H，ZUN I，et al.Transverse migration of single bubbles in simple shear flows［J］.Chem Eng Sci，2002，57（11）：1 849-1 858.

［15］JENS W H，LOTTES P A.Analysis of heat transfer，burnout，pressure drop and density data for high pressure water，ANL-4627［R］.Michigan：ANL，1951.

［16］BARTOLOMEJ G G，CHANTURIYA V M.Experimental study of true void fraction when boiling subcooled water in vertical tubes［J］.Thermal Engineering，1967，14（2）：123-128.