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

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

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

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

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

1 數(shù)學物理模型

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

1.1 子通道模型

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

連續(xù)方程：

(1)

式中：i、j為相鄰通道編號；Ai為通道i流通截面積，m2；ρi為通道i流體密度，kg/m3；mi為通道i軸向質量流量，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型臨界熱流密度機理模型

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

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

微液層蒸干所對應的熱流密度可表示為：

qDNB=ρfδhfgUB/LB

(5)

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

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

FB=FD

(6)

最終可得到汽塊速度的表達式為：

(7)

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

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

(8)

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

(9)

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

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

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

1.3 干涸型臨界熱流密度機理模型

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

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

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

(10)

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

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

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

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

(11)

(12)

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

2) 液滴的夾帶和沉積率

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

md=kdC

(13)

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

(14)

(15)

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

me=me1+me2

(16)

對于前者，Okawa提出了如下關系式：

(17)

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

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

(18)

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

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

mv=q/hfg

(19)

2 子通道模型與機理模型的耦合

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

2.1 子通道模型與DNB型機理模型的耦合

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

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

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

(20)

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

2.2 子通道模型與干涸型機理模型的耦合

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

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

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

(21)

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

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

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

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

3 驗證和分析

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

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

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

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

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

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

CL=a1Rea2

(22)

而Beyerlein等[18]建議CL是平均空泡份額和Re的函數(shù)，即它與紊流波動和當?shù)仄莘植记闆r有關?；诖?，本文中模型采用的CL表達式為：

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

(23)

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

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

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

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

4 結論

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

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