棒束子通道CHF機(jī)理模型開發(fā)及初步驗證

桂民洋，田文喜，吳 迪，陳榮華，張 魁，蘇光輝，秋穗正

(西安交通大學(xué) 核科學(xué)與技術(shù)學(xué)院，陜西 西安 710049)

在核反應(yīng)堆的設(shè)計與安全分析中，臨界熱流密度是重要的限制性熱工參數(shù)，它直接影響核反應(yīng)堆的安全性和經(jīng)濟(jì)性。臨界熱流密度是指沸騰傳熱機(jī)理發(fā)生變化而使核燃料元件表面發(fā)生傳熱惡化時發(fā)熱元件表面單位面積產(chǎn)生的熱量，如果反應(yīng)堆燃料元件表面發(fā)生沸騰臨界，將會導(dǎo)致燃料元件表面溫度過高從而造成加熱壁面燒毀，放射性物質(zhì)泄漏，進(jìn)而造成嚴(yán)重的反應(yīng)堆運行事故[1]。因此，準(zhǔn)確預(yù)測臨界熱流密度對反應(yīng)堆熱工水力設(shè)計和安全分析有著十分重要的意義。臨界熱流密度模型一般分為兩類：偏離泡核沸騰(DNB)型和干涸型。常用的預(yù)測方法主要有經(jīng)驗關(guān)系式和臨界熱流密度查詢表(CHF look-up table)[2]。

但受限于特定流體種類、熱工邊界、幾何尺寸等實驗條件，大多數(shù)實驗關(guān)系式的適用范圍不能輕易外推，或外推后預(yù)測精度下降。因此，一些學(xué)者轉(zhuǎn)而從臨界發(fā)生時的微觀機(jī)理入手，力圖獲得適用范圍更廣的臨界熱流密度機(jī)理模型。如Weisman等[3]、Ying等[4]提出的近壁面氣泡壅塞模型，他們認(rèn)為，發(fā)生臨界時壁面附近大量汽泡聚集形成一個汽泡層，當(dāng)汽泡層外沿處的液相流動減弱到不能穿透汽泡層到達(dá)加熱壁面時，汽泡層的汽泡會結(jié)合形成汽泡膜阻隔傳熱從而發(fā)生DNB；而Lee等[5]、Liu等[6]針對DNB提出了微液層蒸干模型，認(rèn)為壁面附近會有由多個小汽泡結(jié)合而形成一個拉長的大汽塊，汽塊下存在一層很薄的微液層，當(dāng)微液層蒸干時即發(fā)生臨界。對于干涸，其主要發(fā)生在環(huán)狀流末期，此時液滴的作用顯著，因而大多數(shù)學(xué)者主要關(guān)注環(huán)狀流流型區(qū)液滴的夾帶沉積現(xiàn)象[7-12]。

作為反應(yīng)堆堆芯熱工水力分析的重要工具，子通道程序中臨界熱流密度基本通過實驗關(guān)系式求得，而本文則基于臨界熱流密度機(jī)理模型，通過與子通道分析方法耦合，實現(xiàn)對棒束通道內(nèi)臨界熱流密度的預(yù)測，其中子通道程序用于計算通道中的焓分布和流量分布，為臨界熱流密度機(jī)理模型提供邊界條件，進(jìn)而判斷臨界的發(fā)生。

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

本研究中考慮3類模型，即子通道模型、DNB型和干涸型臨界熱流密度機(jī)理模型。其中，子通道模型基于SACOS子通道程序[13]。

1.1 子通道模型

子通道內(nèi)流體的流動包含軸向和橫向兩個方向，其控制方程組(連續(xù)方程、能量方程、軸向動量方程和橫向動量方程)如下。

連續(xù)方程：

(1)

式中：i、j為相鄰?fù)ǖ谰幪枺籄i為通道i流通截面積，m2；ρi為通道i流體密度，kg/m3；mi為通道i軸向質(zhì)量流量，kg/s。

能量方程：

(2)

軸向動量方程：

(3)

式中：pi為通道i的壓力，Pa；g為重力加速度，m/s2；θ為通道傾斜角度，(°)；fT為橫向動量因子；ui為通道i的軸向流速，m/s；uj為通道j的軸向速度，m/s；u*為通道間橫向攪混流體流速，m/s；fi為通道i的摩擦系數(shù)；Di為通道i的等效水力直徑，m；Ks為局部阻力系數(shù)；Δz為軸向控制體高度，m。

橫向動量方程：

(4)

式中：sij為通道i與j間隙長度，m；lij為通道i與j間攪混長度，m；KG為通道i與j間橫流阻力系數(shù)；ρ*為通道間橫流冷卻劑密度，kg/m3，為來流通道的冷卻劑參數(shù)。

1.2 DNB型臨界熱流密度機(jī)理模型

DNB型臨界熱流密度機(jī)理模型主要基于微液層蒸干模型，假設(shè)加熱壁面附近產(chǎn)生的小汽泡結(jié)合形成大汽塊，在汽塊下存在非常薄的液相層，稱為微液層，當(dāng)微液層蒸干即發(fā)生臨界，如圖1所示。

圖1 微液層蒸干模型示意圖

微液層蒸干所對應(yīng)的熱流密度可表示為：

qDNB=ρfδhfgUB/LB

(5)

式中：δ為汽塊下微液層厚度，m；hfg為汽化潛熱，kJ/kg；UB為汽塊移動速度，m/s；LB為汽塊長度，m。

在微液層蒸干模型中，汽塊下微液層厚度δ、汽塊移動速度UB與汽塊長度LB是模型求解的3個關(guān)鍵參數(shù)。其中，汽塊移動速度可通過軸向方向施加在汽塊上的浮力FB和拖拽力FD的平衡得到，即：

FB=FD

(6)

最終可得到汽塊速度的表達(dá)式為：

(7)

式中：UBL為汽塊徑向位置處的主流速度，m/s；CD為拖拽系數(shù)。

汽塊的長度LB假設(shè)為Helmholtz臨界波長，可得：

(8)

微液層的厚度δ由徑向方向施加在汽塊上的力的平衡來確定，主要考慮蒸發(fā)力FE、側(cè)面提升力FL[5]、壁面潤滑力FWL和Marangoni力FM[6]。最終可得其表達(dá)式為：

(9)

式中：DB為汽塊當(dāng)量直徑，m；q為加熱壁面熱流密度，W/m2；CWL為壁面潤滑系數(shù)；CL為提升力系數(shù)；G為通道質(zhì)量流密度，kg/(m2·s)。

圖2為微液層蒸干模型計算流程圖。

圖2 微液層蒸干模型計算流程圖

1.3 干涸型臨界熱流密度機(jī)理模型

干涸型臨界熱流密度機(jī)理模型主要基于環(huán)狀流液膜蒸干模型。加熱通道內(nèi)飽和沸騰區(qū)兩相流流型如圖3所示，主要經(jīng)歷了泡狀流、攪混流，最后形成環(huán)狀流直到液膜逐漸燒干。通道內(nèi)環(huán)狀流的特點是液膜在沿著通道壁面流動，同時夾帶有液滴的汽芯在通道中央流動。

圖3 環(huán)狀流液膜蒸干模型示意圖

在環(huán)狀流區(qū)域存在液膜蒸發(fā)、液滴夾帶和液滴沉積的復(fù)雜現(xiàn)象，干涸的發(fā)生可認(rèn)為是三者之間相互競爭的結(jié)果?？紤]沿通道軸向液膜質(zhì)量流量的變化，可得到液膜質(zhì)量流量的控制方程：

(10)

式中：Wf為液膜質(zhì)量流量，kg/s；md為液滴沉積率，kg/(m2·s)；me為液滴夾帶率，kg/(m2·s)；mv為液相蒸發(fā)率，kg/(m2·s)；P為通道濕周(加熱周長)，m。

一旦確定了環(huán)狀流起始點，通過引入液滴夾帶沉積的本構(gòu)關(guān)系式，由式(10)可得到沿軸向壁面液膜質(zhì)量流量的變化趨勢，當(dāng)液膜質(zhì)量流量為0時，即認(rèn)為發(fā)生干涸。

1) 環(huán)狀流起始點

環(huán)狀流起始點采用Wallis[14]的理論，其推薦如下由攪混流到環(huán)狀流轉(zhuǎn)變的關(guān)系式：

(11)

(12)

即可通過通道內(nèi)含汽率確定環(huán)狀流起始點。

2) 液滴的夾帶和沉積率

液滴的夾帶率和沉積率關(guān)系式選用Okawa模型[9]，它是基于大量的液滴實驗數(shù)據(jù)擬合而來，對于沉積率md，其與氣相中液滴濃度C有關(guān)：

md=kdC

(13)

式中，kd為液滴沉積系數(shù)，表達(dá)式為：

(14)

(15)

針對液滴夾帶率，從夾帶產(chǎn)生的機(jī)理出發(fā)，主要考慮兩類效應(yīng)：中心氣相對液相的剪切夾帶(液相雷諾數(shù)Ref>320)和氣液界面處氣泡的破裂造成的夾帶，即：

me=me1+me2

(16)

對于前者，Okawa提出了如下關(guān)系式：

(17)

式中：fi為相間剪切應(yīng)力系數(shù)；Jg為氣相表觀速度，m/s；δ為液膜厚度，m；σ為表面張力，N/m。

對于后者，采用了Ueda等[15]建立的夾帶關(guān)系式：

(18)

3) 液相的蒸發(fā)率

模型中假設(shè)環(huán)狀流區(qū)域液膜處于飽和狀態(tài)，則加熱壁面的熱量全部用于液相的蒸發(fā)，則液相的蒸發(fā)率為：

mv=q/hfg

(19)

2 子通道模型與機(jī)理模型的耦合

子通道程序已被證明可較好地計算棒束通道中的流量和焓，得到通道中的流量和溫度分布；而臨界熱流密度機(jī)理模型從微觀機(jī)理現(xiàn)象入手，可克服經(jīng)驗關(guān)系式對幾何條件、熱工邊界條件的限制要求，通過兩者的耦合計算，可充分發(fā)揮各自的作用。但一般的機(jī)理模型是基于單通道提出的，對于子通道還需考慮通道間的質(zhì)量、能量攪混，因而對于兩者耦合還需做特殊的處理。

2.1 子通道模型與DNB型機(jī)理模型的耦合

一般情況下，DNB主要發(fā)生在過冷沸騰區(qū)，若通道軸向均勻加熱，則DNB最先發(fā)生在出口；但通道軸向功率非均勻分布時，DNB在通道出口之前某處就有可能發(fā)生，因此，為充分考慮非均勻加熱，軸向通道將被分為若干個控制體，由子通道程序計算得到通道的流量、焓參數(shù)傳遞給機(jī)理模型，每個控制體作為單獨通道均調(diào)用機(jī)理模型進(jìn)行計算?？紤]到徑向功率的不均勻性(如1個子通道為4根加熱棒中間圍成的區(qū)域)，保守來看應(yīng)取徑向熱流密度最大值與式(5)所求的qDNB進(jìn)行比較，從而判斷該位置處是否發(fā)生DNB，耦合后的程序計算流程如圖4所示。

圖4 DNB模型與子通道分析的耦合計算

特別地，為考慮通道間能量攪混，模型中熱流密度需由進(jìn)出口焓進(jìn)行修正，即：

(20)

式中：qj為軸向控制體j中熱流密度，W/m2；hout為軸向控制體j出口焓，J/kg；hin為軸向控制體j入口焓，J/kg；G為軸向控制體j質(zhì)量流量，kg/(m2·s)；A為軸向控制體j流通面積，m2；Det為軸向控制體j加熱直徑，m；l為軸向控制體j軸向長度，m。

2.2 子通道模型與干涸型機(jī)理模型的耦合

液膜蒸干模型是基于單通道(如圓管、矩形通道)提出的，其與子通道模型的耦合需考慮子通道特殊的結(jié)構(gòu)形式，如圖5所示。典型的子通道包含角通道、邊通道和中心通道，一方面，機(jī)理模型需考慮通道間的質(zhì)量、能量攪混；另一方面，需考慮徑向功率的非均勻性，如圖3中的角通道和邊通道，由于存在非加熱壁面，因此壁面液膜軸向質(zhì)量方程中不存在蒸發(fā)項mv，對于中心通道，由于其周圍加熱棒功率的不一致性，導(dǎo)致液膜厚度不同，也需分開考慮。

圖5 棒束中液膜蒸干模型示意圖

因此，考慮干涸型機(jī)理模型在子通道中的應(yīng)用，式(10)液膜質(zhì)量流量的控制方程需改為以下形式：

(21)

式中：k為液膜標(biāo)識(k=1～N)；wl為子通道間液相橫向流質(zhì)量流量，kg/(m·s)；ζf為橫向流中液膜所占份額。

1個子通道中液膜可分為N份(對于角通道N=2，對于邊通道N=3，對于中心通道N=4)，它們共用同一區(qū)域的氣相和液相，對于液相橫向流，可由子通道程序計算得到，而液相又可分為壁面流動液膜和中心夾帶的液滴，程序中兩者的比例ζf可認(rèn)為和通道中液膜與液滴的濃度比相同。

總地來說，由子通道程序計算得到各通道環(huán)狀流起始點，可認(rèn)為起始點處液膜均勻覆蓋在各壁面上，同時由子通道程序可確定通道間流體的橫向流動，即可通過液膜質(zhì)量守恒方程及相關(guān)本構(gòu)方程確定各壁面軸向液膜質(zhì)量流量變化趨勢，當(dāng)某一壁面液膜質(zhì)量流量或空泡份額趨于0時即認(rèn)為發(fā)生干涸，耦合后的程序計算流程如圖6所示。

圖6 蒸干模型與子通道分析的耦合計算

3 驗證和分析

子通道程序與臨界熱流密度機(jī)理模型耦合的驗證基于EPRI的CHF實驗數(shù)據(jù)庫[16-17]，實驗工況包含壓水堆(PWR)和沸水堆(BWR)，可用于DNB型和干涸型臨界熱流密度機(jī)理模型的驗證，通道類型主要為4×4棒束和5×5棒束，具體的通道幾何與熱工邊界參數(shù)列于表1。

表1 臨界熱流密度機(jī)理模型驗證數(shù)據(jù)

本文中，首先由SACOS子通道程序建立棒束通道的幾何模型，給定邊界條件，包括出口壓力、入口流量、入口溫度和平均功率密度，通過子通道程序計算通道各處局部參數(shù)，進(jìn)而調(diào)用DNB或干涸型臨界熱流密度機(jī)理模型計算，逐漸增大平均功率密度，直至發(fā)生臨界。圖7～10示出了模型計算值與實驗值的對比結(jié)果，其中圖7、8為DNB型臨界熱流密度的對比結(jié)果，圖9、10為干涸型臨界熱流密度的對比結(jié)果。

從圖7可看出，驗證數(shù)據(jù)點為110個，94.4%的數(shù)據(jù)點誤差在20%之內(nèi)，機(jī)理模型對CHF的預(yù)測趨勢與實驗結(jié)果一致。從圖8可看出，機(jī)理模型的預(yù)測結(jié)果滿足一般規(guī)律，且與實驗結(jié)果相比有較好的一致性，臨界熱流密度(q)隨系統(tǒng)壓力的升高、入口過冷度的增大和入口流量的增大而增大。本文的DNB型臨界熱流密度機(jī)理模型是基于Lee和Mudawar[5]最早提出的微液層蒸干模型改進(jìn)而來，其中微液層的存在也被一些可視化實驗所證實。但機(jī)理模型中仍存在經(jīng)驗系數(shù)需人為確定，這也是模型計算誤差的主要原因。其主要表現(xiàn)在式(9)中的提升力系數(shù)CL，Lee和Mudawar提出了僅與雷諾數(shù)Re相關(guān)的CL表達(dá)式：

圖7 DNB型臨界熱流密度總體驗證結(jié)果

圖8 DNB型臨界熱流密度的變化

CL=a1Rea2

(22)

而Beyerlein等[18]建議CL是平均空泡份額和Re的函數(shù)，即它與紊流波動和當(dāng)?shù)仄莘植记闆r有關(guān)?；诖耍疚闹心Ｐ筒捎玫腃L表達(dá)式為：

CL=230Re-0.35-0.23exp(1.8α)

(23)

式中，CL的經(jīng)驗處理使得模型計算值存在偏差，且具有一定的不確定性。可預(yù)料，當(dāng)不斷修正上式中的常數(shù)時，可獲得與實驗值更為接近的預(yù)測結(jié)果。

從圖9可看出，驗證數(shù)據(jù)點為696個，95.8%數(shù)據(jù)點誤差在20%之內(nèi)，機(jī)理模型對CHF的預(yù)測趨勢與實驗結(jié)果一致，顯示了模型較好的預(yù)測能力。從圖10可看出，臨界熱流密度隨入口過冷度和流量的增大而增大，隨系統(tǒng)壓力的增大而減小，與實驗趨勢相同?？傮w來看，干涸型機(jī)理模型可較為真實地反映棒束通道中壁面附著的液膜逐漸減薄至發(fā)生干涸的物理過程，兩者存在的偏差有兩個方面的原因：一方面，機(jī)理模型中液滴夾帶沉積模型存在計算偏差，目前關(guān)于液滴的夾帶沉積模型多基于有限的實驗數(shù)據(jù)，且受限于特定的適用范圍，關(guān)系式的擴(kuò)展性較差，因而引入了一定的不確定性，另一方面，目前干涸型臨界熱流密度機(jī)理模型中并未考慮格架的影響，通常格架的存在會增強液滴的橫向運動，進(jìn)而改變液滴的沉積率，因此會影響臨界熱流密度。

圖9 干涸型臨界熱流密度總體驗證結(jié)果

圖10 干涸型臨界熱流密度的變化

4 結(jié)論

本研究以子通道模型為基礎(chǔ)，通過耦合DNB型和干涸型臨界熱流密度機(jī)理模型，實現(xiàn)了對棒束通道內(nèi)臨界熱流密度的預(yù)測，其中借助于子通道程序計算通道中的焓分布和流量分布，為臨界熱流密度機(jī)理模型提供邊界條件，進(jìn)而判斷臨界的發(fā)生。通過與臨界熱流密度實驗數(shù)據(jù)對比發(fā)現(xiàn)，耦合程序?qū)Π羰ǖ乐信R界熱流密度具有較好的預(yù)測精度，臨界熱流密度預(yù)測值隨相關(guān)熱工參數(shù)的變化規(guī)律與實驗相比也有較好的一致性。

本研究中，基于棒束子通道的臨界熱流密度機(jī)理模型可有效彌補目前實驗關(guān)系式適用范圍受限的不足，也可為后續(xù)臨界熱流密度實驗工況設(shè)計提供參考，但需注意的是，目前的研究中臨界熱流密度模型的驗證范圍仍有限，因此在后續(xù)的研究中有必要進(jìn)行更為充分的驗證。同時，對于干涸現(xiàn)象，中心氣相、液膜和液滴的理論被普遍接受，除液滴的夾帶和沉積關(guān)系式外，模型中的不確定因素較少，而對于DNB現(xiàn)象，除微液層蒸干模型，還有近壁面汽泡壅塞、邊界層分離和界面提升等理論，模型中的不確定性參數(shù)也較多，因此后續(xù)的研究也需進(jìn)一步探討臨界熱流密度機(jī)理模型的可靠性。