趙大衛(wèi) 熊萬玉 劉文興 昝元峰

摘 要：基于微液層蒸干DNB型沸騰臨界機(jī)理，本文構(gòu)建了非均勻加熱下的臨界熱流密度（CHF）預(yù)測理論模型。對上游非均勻熱流進(jìn)行沿程積分確定臨界點(diǎn)的來流工況；采用臨界點(diǎn)當(dāng)?shù)責(zé)崃髅芏扰卸ǚ序v危機(jī)的觸發(fā)；從而實(shí)現(xiàn)上游非均勻加熱累積效應(yīng)及臨界觸發(fā)點(diǎn)當(dāng)?shù)鼐植刻匦缘慕Y(jié)合。采用軸向截尾余弦功率分布下矩形通道內(nèi)的CHF實(shí)驗(yàn)數(shù)據(jù)進(jìn)行驗(yàn)證分析：模型預(yù)測的臨界觸發(fā)位置及當(dāng)?shù)谻HF值與實(shí)驗(yàn)結(jié)果符合較好，90%工況下的臨界熱流密度偏差控制在±30%以內(nèi)。對比功率因子修正法的預(yù)測結(jié)果發(fā)現(xiàn)：臨界位置預(yù)測的延后及功率分布影響因子的高估使得功率因子修正法得到的當(dāng)?shù)谻HF值明顯低于非均勻臨界實(shí)驗(yàn)值。

關(guān)鍵詞：臨界熱流密度 截尾余弦功率曲線 微液層蒸干模型

中圖分類號：TK123 文獻(xiàn)標(biāo)識碼：A 文章編號：1674-098X（2015）08（a）-0114-03

臨界熱流密度（CHF）是核燃料元件表面發(fā)生傳熱惡化的熱通量，是冷卻劑流動沸騰機(jī)理發(fā)生突然轉(zhuǎn)變的結(jié)果。因此，臨界熱流密度（CHF）是保證反應(yīng)堆安全運(yùn)行的重要熱工水力限制參數(shù)。國內(nèi)外已開展了大量臨界熱流密度研究，之前所展開的CHF實(shí)驗(yàn)及理論研究絕大部分集中于均勻加熱下的CHF特性。但在反應(yīng)堆運(yùn)行過程中，堆芯內(nèi)的實(shí)際軸向功率分布近似呈截尾余弦的非均勻曲線形式。為了給反應(yīng)堆設(shè)計及安全分析提供更為精確的輸入條件，針對非均勻加熱下的CHF研究一直為核反應(yīng)堆熱工水力研究領(lǐng)域的重要關(guān)注點(diǎn)。

在之前所開展的非均勻CHF研究中，功率因子修正法[1]基于質(zhì)量和能量守衡關(guān)系式和不同的沸騰危機(jī)觸發(fā)機(jī)理模型進(jìn)行推導(dǎo)得到的上游功率非均勻分布形式對當(dāng)?shù)嘏R界熱流密度的影響系數(shù)FC，通過修正均勻加熱下的CHF預(yù)測值qcr，EU得到非均勻CHF預(yù)測值qcr，non。文中將基于微液層蒸干模型假設(shè)構(gòu)建非均勻加熱下的CHF預(yù)測理論模型。通過考慮臨界觸發(fā)點(diǎn)上游非均勻加熱對來流工況的累積效應(yīng)，并結(jié)合當(dāng)?shù)鼐植繜崃魈匦缘膶?shí)現(xiàn)對非均勻加熱下CHF預(yù)測。采用非均勻加熱下矩形窄縫通道內(nèi)的CHF實(shí)驗(yàn)數(shù)據(jù)對上述兩種預(yù)測方法進(jìn)行對比分析。

1 非均勻加熱下的微液層蒸干模型

1.1 模型假設(shè)

Lee和Mudawar[2] 構(gòu)建微液層蒸干模型進(jìn)行DNB型沸騰臨界的預(yù)測。本文對該模型的基本假設(shè)進(jìn)行補(bǔ)充，將其擴(kuò)展至非均勻加熱下的臨界預(yù)測中。

（1）在流道凈蒸汽產(chǎn)生點(diǎn)（NVG）之后，汽泡在脫離壁面后相互融合形成汽塊，汽塊在壁面非常薄的液膜上滑移，汽塊的當(dāng)量直徑大致等于壁面上汽泡脫離點(diǎn)處的汽泡直徑DB。

（2）汽塊滑移速度UB為主流速度和汽塊相對速度的疊加，汽塊相對速度由加載在汽塊上的浮力和拖拽力相平衡得到。

（3）汽塊的長度LB受限于Helmholtz臨界波長。

（4）臨界點(diǎn)上游的非均勻熱流對主流體的貢獻(xiàn)不斷向下流傳遞累積確定了臨界點(diǎn)的來流工況，而沸騰臨界的觸發(fā)由當(dāng)?shù)責(zé)崃魉鶝Q定。

沸騰臨界所對應(yīng)的臨界熱流密度可以表示為：

（1）

1.2 模型描述

1.2.1 汽塊軸向受力和汽塊移動速度

汽塊的移動速度UB通過軸向方向施加在汽塊上的浮力FB和拖拽力FD間的平衡確定[2]：

（2）

以上各式中：采用Chan和Prince推薦的適用于高壓條件下拖拽系數(shù)CD表達(dá)式[3]：

（3）

壁面邊界層內(nèi)的速度分布使用Karman速度分布方程。近壁面區(qū)的汽塊一般處在緩沖區(qū)的范圍，為簡化計算，本文模型直接采用緩沖區(qū)速度公式計算當(dāng)?shù)亓魉俜植?，可得距離壁面y= δ+ DB/2處的當(dāng)?shù)亓魉倏梢员磉_(dá)為：

（4）

1.2.2 汽塊直徑和汽塊長度

汽塊的當(dāng)量直徑DB由Levy模型計算如下[4]：

（5）

汽塊的長度LB為Helmholtz臨界波長：

（6）

式中Usb是微液層中液體流速，由于加熱壁面附近的微液層非常薄且微液層中液體流速相對汽塊移動速度非常小，可假設(shè)Usb等于零。

1.2.3 汽塊徑向受力和微液層厚度

Lee和Mudawar提出了兩個施加在汽塊上的相反方向的力，分別為蒸發(fā)力和側(cè)面提升力[2]。近壁面滑移汽塊還需考慮壁面潤滑力FWL和Marangoni力FM[5、6]。汽塊下微液層厚度δ通過施加在汽塊上的徑向力的平衡來計算：

（7）

求出δ以后，可以通過式（4）計算得到UBL，然后分別用式（2）和（6）計算得到UB和LB，再將新的UB和LB代入式（7）計算得到新的δ，通過迭代計算直到以上幾個參數(shù)值達(dá)到收斂。最后，將迭代計算得到的δ，UB和LB代入式（1）來計算臨界熱流密度值。

2 0程序流程

在軸向熱流非均勻分布條件下，當(dāng)?shù)責(zé)崃髅芏萹loc（z）可以表示為：

（8）

式中：

qave為整個流道上的平均熱流密度/kW·m-2；

Fp（z）為沿程功率分布因子。

對于非均勻加熱工況，臨界觸發(fā)點(diǎn)并不一定只出現(xiàn)在流道出口處，而是隨入口工況的改變，在功率峰至出口間不斷遷移。因此非均勻加熱下的臨界位置確定，需要將假定臨界觸發(fā)點(diǎn)zpcr從功率峰至出口方向不斷遷移，通過比較臨界觸發(fā)時所對應(yīng)的最小平均熱流密度，來搜尋臨界觸發(fā)點(diǎn)。

在假定臨界觸發(fā)點(diǎn)zpcr處，對應(yīng)臨界熱流密度qloc，cr（zpcr）計算流程如下。

（1）基于初始所假設(shè)的整個實(shí)驗(yàn)段上平均熱流密度qave0計算初始功率分布qloc（z）；將假定臨界觸發(fā)點(diǎn)zpcr上游的局部功率進(jìn)行積分，求得臨界觸發(fā)點(diǎn)的來流工況。

（2）在臨界觸發(fā)點(diǎn)處的汽塊及液膜參數(shù)δ， UB， LB的計算中采用臨界觸發(fā)點(diǎn)的當(dāng)?shù)責(zé)崃髅芏萹loc（zpcr）。

（3）基于沸騰危機(jī)觸發(fā)判別公式（1）求得在zpcr處觸發(fā)臨界所需當(dāng)?shù)責(zé)崃髅芏萹，loc，cr（zpcr）。

（4）判斷臨界觸發(fā)熱流q，loc，cr（zpcr）與初始假設(shè)當(dāng)?shù)責(zé)崃鱭loc（zpcr）間的偏差，若當(dāng)?shù)責(zé)崃髌畲笥谠O(shè)定值，返回步驟（1），重新設(shè)定平均熱流密度，循環(huán)迭代（1）-（4）直至當(dāng)?shù)責(zé)崃髌盥淙朐O(shè)定偏差。

（5）當(dāng)?shù)責(zé)崃髅芏萹，loc，cr（zpcr）即為在zpcr處觸發(fā)臨界所對應(yīng)熱流密度，同時求出實(shí)驗(yàn)段此時所對應(yīng)平均熱流密度qave（zpcr）。依次求解功率峰至實(shí)驗(yàn)段出口間各點(diǎn)觸發(fā)臨界所對應(yīng)當(dāng)?shù)責(zé)崃髅芏萹，loc，cr（zpcr）及平均熱流密度qave（zpcr）。其中最小平均熱流密度qave（zpcr）所對應(yīng)位置，即為非均勻加熱下的臨界觸發(fā)位置，該點(diǎn)所對應(yīng)當(dāng)?shù)責(zé)崃髅芏燃礊榕R界熱流密度qloc，cr。

3 結(jié)果驗(yàn)證及分析

為了驗(yàn)證本文所構(gòu)建的非均勻加熱CHF理論模型，采用矩形通道內(nèi)的96組非均勻加熱CHF實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對比分析。實(shí)驗(yàn)中，沿豎直流道方向上的軸向功率曲線為典型截尾余弦曲線，其最大功率因子為1.55。非均勻加熱CHF實(shí)驗(yàn)數(shù)據(jù)和非均勻加熱液膜蒸干模型及修正因子兩種方法的預(yù)測結(jié)果如圖1所示。其中修正因子法采用Tong功率因子關(guān)系式。

Tong[2]基于沸騰危機(jī)觸發(fā)時近壁面汽泡層內(nèi)焓升恒定的假設(shè)推導(dǎo)了非均勻加熱下CHF修正因子關(guān)系式：

（9）

非均勻加熱時，實(shí)驗(yàn)及預(yù)測的CHF觸發(fā)位置均分布在功率峰值至實(shí)驗(yàn)段出口之間的區(qū)域；非均勻加熱液膜蒸干模型預(yù)測的臨界觸發(fā)位置及當(dāng)?shù)谻HF值均與實(shí)驗(yàn)結(jié)果相吻合；修正因子法所得到的臨界觸發(fā)位置相對滯后，沸騰臨界觸發(fā)位置基本位于實(shí)驗(yàn)流道的出口，所預(yù)測當(dāng)?shù)谻HF值也偏小。

非均勻加熱時，實(shí)驗(yàn)及預(yù)測的CHF觸發(fā)位置均分布在功率峰值至實(shí)驗(yàn)段出口之間的區(qū)域；非均勻加熱液膜蒸干模型預(yù)測的臨界觸發(fā)位置及當(dāng)?shù)谻HF值均與實(shí)驗(yàn)結(jié)果相吻合；修正因子法所得到的臨界觸發(fā)位置相對滯后，沸騰臨界觸發(fā)位置基本位于實(shí)驗(yàn)流道的出口，所預(yù)測當(dāng)?shù)谻HF值也偏小。

在圖2中非均勻加熱下沸騰臨界點(diǎn)當(dāng)?shù)谻HF值的預(yù)測結(jié)果對比中：非均勻加熱液膜蒸干模型對90%實(shí)驗(yàn)工況的預(yù)測結(jié)果能夠落入±30%的誤差范圍內(nèi)。而修正因子法的大部分預(yù)測值明顯低于實(shí)驗(yàn)中的當(dāng)?shù)谻HF，僅有部分高質(zhì)量流速下的實(shí)驗(yàn)工況預(yù)測誤差在±30%之內(nèi)。對于截尾余弦功率分布的加熱方式，軸向功率因子在流道中段的最高功率峰之后是不斷下降的。圖1中，修正因子法對臨界位置預(yù)測值明顯偏后，進(jìn)而帶來所對應(yīng)當(dāng)?shù)毓β室蜃拥南陆怠Ｔ谙嗤骄鶡崃髅芏认?，使得修正因子法預(yù)測臨界觸發(fā)位置所對應(yīng)當(dāng)?shù)責(zé)崃髅芏绕汀榱伺卸ㄐ拚蜃臃A(yù)測結(jié)果偏低是否僅是臨界位置預(yù)測滯后所引起的，圖3給出了非均勻加熱下臨界觸發(fā)對應(yīng)平均熱流密度的預(yù)測結(jié)果。

從圖3臨界觸發(fā)對應(yīng)平均熱流密度的預(yù)測對比結(jié)果可以看出，非均勻加熱液膜蒸干模型中所預(yù)測平均熱流密度的偏差在±10%的范圍內(nèi)。而修正因子法得到的平均熱流密度均低于實(shí)驗(yàn)值，臨界位置預(yù)測滯疊加平均熱流密度預(yù)測偏低，使得圖2中修正因子法預(yù)測當(dāng)?shù)谻HF值明顯低于實(shí)驗(yàn)值。

在功率因子修正法中，是采用Tong非均勻功率分布影響因子Fc來修正均勻加熱下的CHF預(yù)測值qcr，EU得到非均勻CHF預(yù)測值qcr，non：

（10）

其中，均勻加熱CHF由Debortoli關(guān)系式[7]進(jìn)行計算。圖4中給出了與非均勻加熱CHF實(shí)驗(yàn)工況范圍基本一致下的均勻加熱CHF實(shí)驗(yàn)結(jié)果與Debortoli關(guān)系式預(yù)測結(jié)果的對比。雖然Debortoli關(guān)系式預(yù)測CHF在±20%的偏差范圍內(nèi)分布較為發(fā)散，但預(yù)測結(jié)果并沒有出現(xiàn)系統(tǒng)的負(fù)偏差?；谛拚蜃臃ㄓ嬎阒械墓剑?0），可以推測采用公式（9）針對本實(shí)驗(yàn)中的非均勻功率分布所得到功率分布影響因子Fc被高估，進(jìn)而使得非均勻CHF預(yù)測值偏低。

4 結(jié)語

在Lee和Mudawar所建立的均勻加熱微液層蒸干模型基礎(chǔ)上，本文采用臨界點(diǎn)上游非均勻熱流積分求得來流工況，由臨界點(diǎn)當(dāng)?shù)責(zé)崃髅芏葋頉Q定壁面處的沸騰危機(jī)觸發(fā)特性，構(gòu)建了預(yù)測非均勻加熱下CHF的理論模型。從而實(shí)現(xiàn)非均勻加熱臨界觸發(fā)中上游非均勻加熱累積效應(yīng)及當(dāng)?shù)鼐植刻匦缘慕Y(jié)合。采用該模型對軸向截尾余弦功率分布下矩形通道內(nèi)的CHF結(jié)果進(jìn)行了預(yù)測：模型預(yù)測的臨界觸發(fā)位置與實(shí)驗(yàn)臨界觸發(fā)位置基本一致，偏差控制在±20%之內(nèi)；在800～3200 kg/m2s的質(zhì)量流速范圍內(nèi)，臨界觸發(fā)點(diǎn)處CHF預(yù)測值與非均勻CHF實(shí)驗(yàn)值的偏差控制在±30%左右，且無系統(tǒng)偏差。對比功率因子修正法的預(yù)測結(jié)果發(fā)現(xiàn)：臨界位置預(yù)測的延后及功率分布影響因子Fc的高估使得功率因子修正法得到的當(dāng)?shù)谻HF值遠(yuǎn)低于非均勻臨界實(shí)驗(yàn)值。針對本實(shí)驗(yàn)中軸向截尾余弦加熱CHF實(shí)驗(yàn)，該文所構(gòu)建非均勻加熱微液層蒸干模型的預(yù)測結(jié)果明顯優(yōu)于功率因子修正法的預(yù)測結(jié)果。

參考文獻(xiàn)

[1]Tong L S. Prediction of departure from nucleate boiling for an axially non-uniform heat flux distribution[J].of Nuclear Energy，1967（21）：51-56.

[2]Lee C.H.， Mudawar I.，A mechanistic critical heat flux model for subcooled flow boiling based on local bulk flow conditions[J].Int of Multiphase Flow，1988（14）：711-728.

[3]Chan B.K.C， Prince R.G.H. Distillation studies-viscous drag on a gas bubble rising in a liquid[J].Aiche Journal，1965（11）：176.

[4]Levy S.Forced convection subcooled boiling - prediction of vapor volumetric fraction， Int[J].Heat Mass Transf，1967（0）：951.

[5]Antal S.P.， Lahey R.T.， Jr.， et al.Analysis of phase distribution in fully-developed laminar bubbly 2-phase flow[J].Int of Multiphase Flow，1991（17）：635-652.

[6]Lahey R.T， Jr， Drew D.A. The three-dimensional time and volume averaged conservation equations of two-phase flow[J]. Advances in nuclear science and technology， Plenum press， New York， NY， USA，1988（20）：1-69.

[7]Debortoli R A， Green S J.Forced convection heat transfer burnout studies for water in rectangular channels and round tubes at pressures above 500 Psia[J].Burnout，1958（9）：188.