白宇飛， 楊星團， 張 震

(1.清華大學(xué) 核能與新能源技術(shù)研究院，北京 100084；2.清華大學(xué) 先進核能技術(shù)協(xié)同創(chuàng)新中心，北京 100084；3.先進反應(yīng)堆工程與安全教育部重點實驗室(清華大學(xué))，北京 100084)

沸騰通道內(nèi)的兩相流動不穩(wěn)定性可能導(dǎo)致有害機械振動，干擾控制系統(tǒng)，使部件熱疲勞并惡化傳熱性能，因而受到廣泛關(guān)注. Boure等[1]、Kakac等[2]、Ruspini等[3]基于大量的實驗與理論研究對該問題進行了系統(tǒng)的分類. 對于沸水反應(yīng)堆和蒸汽發(fā)生器，尤其要注意避免在高出口干度條件下發(fā)生由摩擦壓降主導(dǎo)的第二類密度波振蕩. 在臨界的反應(yīng)堆內(nèi)中子通量密度的空間分布是不均勻的，對于動力堆中常見的圓柱形堆芯，其軸向功率服從余弦函數(shù)分布，其徑向功率服從貝塞爾函數(shù)分布. 盡管在實際中采取各類功率分布展平措施，其熱中子通量密度分布不均勻系數(shù)仍在1.4左右[4]. 此外，由于各種原因，冷卻劑在堆芯內(nèi)數(shù)量龐大的平行通道中的流量分配并不是均勻的，通常采用增加堆芯入口流動阻力的方式展平進口流量[5].

國內(nèi)外學(xué)者對這些軸向和徑向上的不均勻分布進行了實驗與理論研究. Dijkman等[6]、Dutta等[7]通過數(shù)值模擬研究發(fā)現(xiàn)軸向熱流余弦分布具有強烈的穩(wěn)定效應(yīng). 然而，Narayanan等[8]、Lu等[9]在研究多種軸向熱流分布的基礎(chǔ)上指出軸向熱流余弦分布在低入口過冷度下使系統(tǒng)穩(wěn)定性增強，在高入口過冷度下使系統(tǒng)穩(wěn)定性減弱. 黃軍等[10]、Guo等[11]對低壓下平行通道系統(tǒng)中的流動不穩(wěn)定性進行了實驗研究，分析了兩個通道不均勻加熱和入口節(jié)流的影響；夏庚磊等[12-13]采用RELAP5程序進行數(shù)值模擬，分析了系統(tǒng)參數(shù)的影響. Xia等[14]研究了徑向不均勻加熱和入口節(jié)流對流量分配不穩(wěn)定性的影響，并分析了通道壓降與流量變化的關(guān)系. 上述對軸向熱流不均勻分布的研究所得結(jié)論尚未統(tǒng)一，因此有必要探究其影響流動不穩(wěn)定性的機理.

本文利用RELAP5/MOD4.0程序?qū)Σ⒙?lián)沸水堆子通道系統(tǒng)進行數(shù)值模擬，對比了軸向熱流線性分布及余弦分布工況下的流動不穩(wěn)定邊界，分析了系統(tǒng)和通道的臨界相變數(shù)隨徑向功率和入口節(jié)流系數(shù)不均勻性增大的變化情況，系統(tǒng)地研究不均勻加熱與入口節(jié)流對平行通道內(nèi)密度波不穩(wěn)定性的影響，以期為核動力設(shè)備的設(shè)計和運行提供參考和依據(jù).

1 研究方法

1.1 模型的建立

強迫循環(huán)回路中的兩個并聯(lián)通道作為一種典型的邊界條件被廣泛應(yīng)用于流動不穩(wěn)定性的實驗與理論研究中，其特點是平行通道的總流量保持恒定. 本研究采用基于兩流體、非平衡、非均相兩相流模型的大型熱工水力瞬態(tài)分析程序RELAP5/MOD4.0進行數(shù)值模擬. 程序中采用的基本守恒方程如下.

氣相和液相質(zhì)量守恒方程：

式中：t為時間，s；α為體積分數(shù)；ρ為密度，kg/m3；A為流道橫截面積，m2；x為軸向坐標(biāo)，m；v為流速，m/s；Γ為傳質(zhì)速率，kg/(m3·s)；下標(biāo)g、f分別為氣相液相.

氣相和液相動量守恒方程：

式中：p為壓力，Pa；B為質(zhì)量力，m/s2；FgW和FfW為氣和液相壁面曳力系數(shù)，s-1；vgI和vfI為相界面處氣和液相流速，m/s；FgI和FfI為氣和液相相界面曳力系數(shù)，s-1；C為虛擬質(zhì)量系數(shù).

氣相和液相能量守恒方程：

式中：U為比內(nèi)能，J/kg；q為體積釋熱率，W/m3；h為焓值，J/kg；Φ為能量耗散項，W/m3.

根據(jù)Guo等[11]研究中的實驗裝置建立的系統(tǒng)節(jié)點圖如圖1所示.

圖1 RELAP5程序節(jié)點圖

兩個幾何結(jié)構(gòu)相同的圓形截面豎直通道均由入口段、加熱段和上升段組成，過冷水從系統(tǒng)底部進入，在向上流動的過程中發(fā)生相變，并以一定干度的氣液兩相混合物的狀態(tài)流出. 采用3個由單一接管部件(Single Junction)連接的管型部件(Pipe)模擬每個通道，加熱功率通過內(nèi)部產(chǎn)生熱量的熱構(gòu)件(Heat Structure)控制，采用對流邊界條件使熱量從熱構(gòu)件傳至通道內(nèi)流體；分支部件(Branch)用于模擬上下兩個聯(lián)箱，與通道相連處設(shè)置了流動損失系數(shù)以模擬進出口節(jié)流閥的作用；時間相關(guān)控制體部件(Time-dependent Volume)用于給定入口流體的物性以及出口壓力；時間相關(guān)接管部件(Time-dependent Junction)用于給定系統(tǒng)的入口總質(zhì)量流量. 經(jīng)過節(jié)點敏感性分析，綜合考慮計算精度與耗時，將入口段、加熱段和上升段分別劃分為20、50和10個控制體.

1.2 流動不穩(wěn)定的判斷

本研究采用固定入口總質(zhì)量流量同時線性增大加熱功率的方式使流量發(fā)生振蕩，初始時保持較大的功率增長速度以迅速接近流動不穩(wěn)定邊界，之后采用較小的功率增長速度以利于對流動不穩(wěn)定起始的判定. 經(jīng)驗證，這種方式與階躍增大加熱功率的方式所得流動不穩(wěn)定邊界相同. 典型的加熱功率及通道進出口質(zhì)量流量變化如圖2所示.

系統(tǒng)壓力2 MPa，入口溫度168 ℃，入口總質(zhì)量流量0.053 kg/s，兩通道入口節(jié)流系數(shù)均為37. 由于兩個通道特性相同，功率較低時通道內(nèi)處于相同質(zhì)量流量的穩(wěn)定流動狀態(tài)，隨加熱功率增大，兩通道入口流量發(fā)生等振幅反相振蕩，振幅逐漸增大而周期保持不變，在振蕩充分發(fā)展后通道的入口發(fā)生倒流，不利于核動力裝置的安全可靠運行. 由圖2(b)可知平行通道系統(tǒng)中的單個通道的進出口流量反向振蕩且出口流量的振幅較小. 此外，發(fā)生不穩(wěn)定時通道出口處干度約為0.4. 結(jié)合以上特征可知,通道內(nèi)發(fā)生了第二類密度波振蕩. 參照文獻[11]，本研究選取單個通道入口質(zhì)量流量相對振幅達到10%作為判斷流動不穩(wěn)定發(fā)生的臨界值.

1.3 模擬與實驗數(shù)據(jù)的驗證

本研究采用Ishii等[15]提出的無量綱數(shù)表示流動不穩(wěn)定邊界，相變數(shù)Npch表征出口流體的相變程度，過冷度數(shù)Nsub表征入口流體的過冷度，且

出口定干度線有助于確定系統(tǒng)的運行特性

Nsub=Npch-xeνfg/νf.

式中：Q為加熱功率，kW；M為質(zhì)量流量，kg/s；hfg為汽化潛熱，kJ/kg；νfg為飽和氣液比容差，m3/kg；νf為飽和液相比容，m3/kg；hf為飽和液相焓，kJ/kg；hin為入口流體焓，kJ/kg；xe為出口熱平衡含氣率.

(a)兩通道入口流量變化

(b)平行通道系統(tǒng)中的單個通道進出口流量變化

從可靠性及數(shù)值穩(wěn)定性的角度出發(fā)[16]，本研究選取了RELAP5/MOD4.0程序中的非平衡非均相模型和半隱式數(shù)值格式對Guo等[11]研究中的16組實驗進行數(shù)值模擬，計算結(jié)果與實驗結(jié)果的對比如圖3所示.

圖3 計算結(jié)果與實驗結(jié)果的對比

Fig.3 Comparison of computational results and experimental data

由圖3可知，不同壓力下的流動不穩(wěn)定邊界分布在一個條形區(qū)域內(nèi)，左側(cè)為流動穩(wěn)定區(qū)間，右側(cè)為流動不穩(wěn)定區(qū)間. 大多數(shù)計算結(jié)果中的臨界加熱功率略小于實驗，相對誤差在-0.9%～-14.1%，驗證了模擬的準(zhǔn)確性，也證實了RELAP5/MOD4.0程序模擬平行通道內(nèi)密度波不穩(wěn)定性的有效性與適用性.

2 計算結(jié)果與分析

本文以典型的內(nèi)部流體豎直向上流動的沸水堆子通道[16]作為研究對象，用兩個并聯(lián)的通道代表沸水堆堆芯數(shù)量龐大的平行通道，對軸向、徑向不均勻加熱以及入口節(jié)流工況下的平行通道密度波不穩(wěn)定性進行研究. 其中單個通道的直徑為0.012 4 m，長度為3.657 6 m，壁面粗糙度為0.000 025 m，入口阻力系數(shù)變化范圍為0～46，出口阻力系數(shù)為5. 系統(tǒng)壓力為7 MPa，并設(shè)置系統(tǒng)的總質(zhì)量流量為0.264 kg·s-1. 此物理模型不含絕熱的入口段和上升段，經(jīng)過節(jié)點敏感性分析，將單個受熱通道劃分為48個控制體.

2.1 軸向不均勻加熱

首先對軸向加熱功率線性變化的工況進行研究，以確定軸向熱流不均勻分布影響流動不穩(wěn)定性的機理. 如圖4所示，設(shè)置了4組線性增大和減小的軸向熱流密度分布，其最大熱流與平均熱流的比值分別為1.4和1.2. 線性增大加熱功率使流量發(fā)生振蕩，流動不穩(wěn)定起始過程中兩通道入口流量隨時間的變化與均勻熱流工況類似，得到的流動不穩(wěn)定邊界如圖5所示，其中軸向熱流均勻分布的結(jié)果作為對比. 在不同入口過冷度下，軸向熱流線性增大使系統(tǒng)穩(wěn)定性增強而軸向熱流線性減小使系統(tǒng)穩(wěn)定性減弱. 隨軸向熱流分布斜率逐漸增大，流動不穩(wěn)定邊界對應(yīng)的出口熱平衡含氣率從不到0.3逐漸增大到超過0.9，系統(tǒng)穩(wěn)定性得到了顯著提升. 其根本原因是軸向熱流不均勻分布會影響沸騰邊界的位置，改變單相段及兩相段長度，進而改變單相段壓降占通道總壓降的份額. 圖6為入口溫度210 °C對應(yīng)過冷度數(shù)4.706，加熱功率115 kW時的空泡份額沿通道高度的分布，顯然單相段長度隨軸向熱流分布斜率的增大而增大. 由于總加熱功率一定，不同軸向熱流分布工況下的出口空泡份額相同. 此外，由圖5可得軸向熱流不均勻分布對流動不穩(wěn)定的影響隨入口過冷度的增大而增大. 這是因為計算中沸騰邊界均位于通道的前半部分，且隨入口過冷度增大不斷向通道出口移動，使不同軸向熱流分布下從通道入口到沸騰邊界的熱流密度積分差異增大.

圖4 軸向熱流線性分布

圖5 軸向熱流線性分布對應(yīng)流動不穩(wěn)定邊界

Fig.5 Flow instability boundaries obtained from linear distribution of axial heat flux

圖6 相同功率時空泡份額沿通道的分布

Fig.6 Void fraction along the channel at same heating power

未受擾動的均勻堆芯的軸向功率服從余弦函數(shù)分布，對于實際反應(yīng)堆中堆芯壽期的不同階段，由于控制棒插入的深度不同，功率的峰值可能增大并向堆芯的頂部或底部偏移. 這3種情況下的軸向熱流密度分布如圖7所示，其中峰值偏向入口和出口的熱流密度分布在Leung等[17]研究的基礎(chǔ)上修改得到. 3種工況所得流動不穩(wěn)定邊界與均勻熱流工況的對比如圖8所示. 在不同入口過冷度下，軸向熱流余弦分布對應(yīng)臨界出口熱平衡含氣率均大于均勻加熱，系統(tǒng)穩(wěn)定性增強. 峰值偏向通道入口和出口的熱流分布所對應(yīng)流動不穩(wěn)定邊界分別處于余弦分布的兩側(cè)，因此與軸向熱流余弦分布相比，峰值偏向入口會降低系統(tǒng)穩(wěn)定性而峰值偏向出口會提高系統(tǒng)穩(wěn)定性，這對指導(dǎo)反應(yīng)堆的安全穩(wěn)定運行具有重要意義. 此外，與熱流均勻分布相比，峰值偏向入口的分布在較低入口過冷度時增大系統(tǒng)的穩(wěn)定性，在較高入口過冷度時降低系統(tǒng)的穩(wěn)定性.

圖7 軸向熱流余弦及峰值偏向入口及出口的分布

Fig.7 Cosine, inlet-peaked and outlet-peaked axial heat flux distribution

圖8 軸向熱流余弦及峰值偏離中心的分布對應(yīng)流動不穩(wěn)定邊界

Fig.8 Flow instability boundaries obtained from cosine and peaked axial heat flux distribution

從入口到通道不同位置處的歸一化熱流積分如圖9所示. 在整個通道長度上，熱流分布的峰值越靠近通道出口，則積分熱流越小，功率相同時的單相段越長，所占通道總壓降份額越大，因此系統(tǒng)越穩(wěn)定. 本文不同入口過冷度下的沸騰邊界都位于通道的前半部分，此時余弦分布的積分熱流小于均勻分布，因此余弦加熱的穩(wěn)定性比均勻加熱更好. 此外，在圖8中以均勻與峰值偏向入口熱流分布所得流動不穩(wěn)定邊界的交點作為入口過冷度高低的分界，并在圖9中以這兩種分布對應(yīng)積分熱流的交點將前一半通道分成兩個區(qū)間. 入口過冷度較低時，沸騰邊界處于區(qū)間I，入口過冷度較高時，沸騰邊界處于區(qū)間II，積分熱流曲線可以很好地解釋不同流動不穩(wěn)定邊界的相對位置.

圖9 積分熱流沿通道長度的分布

在Narayanan等[8]、Lu等[9]的研究中，熱流余弦分布在高入口過冷度下使系統(tǒng)穩(wěn)定性減弱，這一結(jié)果與本文并不矛盾. 若沸騰邊界位于通道上半部分，此時余弦分布的積分熱流大于均勻分布，導(dǎo)致系統(tǒng)穩(wěn)定性減弱. 因此軸向熱流不均勻分布時系統(tǒng)的穩(wěn)定性取決于沸騰邊界的位置以及積分熱流沿通道的分布. 在密度波振蕩邊界處，單相段與兩相段的質(zhì)量流量和壓降存在相位差，單相段越長進而其壓降占通道總壓降的份額越大，使兩相段的壓降振幅不能平衡較大的單相段壓降振幅，導(dǎo)致流量振蕩的衰減，系統(tǒng)穩(wěn)定性增強.

2.2 徑向不均勻加熱

在通道軸向熱流均勻分布的基礎(chǔ)上，保持通道二的功率為通道一的1.4倍，兩通道入口流量隨加熱功率的變化如圖10所示(入口溫度255 ℃，對應(yīng)過冷度數(shù)2.015). 加熱功率較低時兩通道內(nèi)處于不同流量的穩(wěn)定流動狀態(tài). 因為加熱功率較大的通道內(nèi)沸騰更加劇烈，在并聯(lián)通道壓降相同的條件下，對應(yīng)流量較小，即兩個通道內(nèi)的流量不均勻分配. 隨加熱功率線性增大，兩個通道內(nèi)流體同時發(fā)生反相振蕩，形成密度波流動不穩(wěn)定.

通道二功率為通道一的1.4倍時的流動不穩(wěn)定邊界如圖11所示，其中系統(tǒng)的臨界相變數(shù)用發(fā)生不穩(wěn)定時的總功率和流量計算，兩個通道的臨界相變數(shù)分別用發(fā)生不穩(wěn)定時各自的功率和流量計算. 還給出了兩通道功率相同時的流動不穩(wěn)定邊界以便對比分析，此邊界同時表征了系統(tǒng)和單個通道的穩(wěn)定性. 由圖11可知，兩通道功率不同時系統(tǒng)的流動不穩(wěn)定邊界與功率相同時較為接近，說明徑向功率不均勻分布對系統(tǒng)穩(wěn)定性的影響并不顯著. 加熱功率較小的通道所對應(yīng)臨界相變數(shù)較小而加熱功率較大的通道所對應(yīng)臨界相變數(shù)較大，其原因是流量的不均勻分配. 顯然，系統(tǒng)的不穩(wěn)定由加熱功率較大的通道導(dǎo)致. 在一些入口過冷度下功率較大的通道出口產(chǎn)生過熱蒸汽，可能發(fā)生換熱惡化甚至導(dǎo)致通道燒毀，應(yīng)盡量避免此情況發(fā)生.

圖10 徑向不均勻加熱時兩通道入口流量變化

Fig.10 Variation of inlet mass flow rates under nonuniform radial heating

圖11 徑向不均勻加熱對應(yīng)流動不穩(wěn)定邊界

Fig.11 Flow instability boundaries under nonuniform radial heating

系統(tǒng)及通道的臨界相變數(shù)隨兩通道功率比的變化如圖12所示(入口溫度210 ℃，對應(yīng)過冷度數(shù)4.706). 由圖12可知，兩個通道內(nèi)的穩(wěn)定性差異隨不均勻性增大而增大，系統(tǒng)的臨界相變數(shù)隨功率比變化不明顯，說明徑向功率不均勻分布對系統(tǒng)的穩(wěn)定性影響較小. 其原因可能是并聯(lián)通道中和了兩個通道的差異，綜合了系統(tǒng)參數(shù). 徑向功率不均勻分布引起的流量不均勻分配及傳熱惡化問題比流動不穩(wěn)定性更值得關(guān)注.

圖12 徑向加熱不均勻性對系統(tǒng)和通道穩(wěn)定性的影響

Fig.12 Effect of nonuniformity of radial heating on stability of system and each channel

2.3 徑向不均勻入口節(jié)流

當(dāng)兩個通道入口節(jié)流系數(shù)不同時，也會發(fā)生類似圖10中的流量不均勻分配情況. 節(jié)流系數(shù)較大的通道內(nèi)流量較小而節(jié)流系數(shù)較小的通道內(nèi)流量較大以保證并聯(lián)通道兩端壓降相同. 將兩個通道的入口節(jié)流系數(shù)分別設(shè)置為23和46，所得系統(tǒng)及通道的流動不穩(wěn)定邊界與兩種入口節(jié)流系數(shù)相等工況的對比如圖13所示. 不均勻節(jié)流時系統(tǒng)穩(wěn)定性好于較小入口節(jié)流系數(shù)的均勻節(jié)流(kin,1=kin,2=23)，差于對應(yīng)較大入口節(jié)流系數(shù)的均勻節(jié)流(kin,1=kin,2=46). 入口節(jié)流系數(shù)較小的通道穩(wěn)定性相對均勻節(jié)流(kin,1=kin,2=23)增強，而入口節(jié)流系數(shù)較大的通道穩(wěn)定性相對均勻節(jié)流(kin,1=kin,2=46)減弱. 顯然，系統(tǒng)的不穩(wěn)定由入口節(jié)流系數(shù)較小的通道導(dǎo)致.

圖13 徑向不均勻入口節(jié)流對應(yīng)流動不穩(wěn)定邊界

Fig13. Flow instability boundaries under nonuniform radial inlet throttling

為研究入口節(jié)流不均勻性的影響，設(shè)置兩通道入口節(jié)流系數(shù)之和為46，減小通道一的節(jié)流系數(shù)同時增大通道二的節(jié)流系數(shù)，所得流動不穩(wěn)定的臨界相變數(shù)如圖14所示(入口溫度210 ℃，對應(yīng)過冷度數(shù)4.706). 兩通道的穩(wěn)定性差異隨入口節(jié)流不均勻性的增大而增大. 系統(tǒng)的穩(wěn)定性介于兩者之間，當(dāng)不均勻性較小時，系統(tǒng)穩(wěn)定性基本不發(fā)生變化，當(dāng)不均勻性較大時，系統(tǒng)穩(wěn)定性有減小的趨勢，此結(jié)果與文獻[18]中的結(jié)論相符. 這是因為不均勻性較大時，改變相同大小的入口節(jié)流系數(shù)對兩個通道的相對影響不同，入口節(jié)流系數(shù)較小的通道一引入的不穩(wěn)定效應(yīng)逐漸占據(jù)主導(dǎo).

圖14 徑向入口節(jié)流不均勻性對系統(tǒng)和通道穩(wěn)定性的影響

Fig.14 Effect of nonuniformity of radial inlet throttling on stability of system and each channel

為進一步探究提高系統(tǒng)穩(wěn)定性的方法，將通道一的入口節(jié)流系數(shù)固定為23，調(diào)節(jié)通道二的入口節(jié)流系數(shù)，所得流動不穩(wěn)定的臨界相變數(shù)如圖15所示(入口溫度210 ℃，對應(yīng)過冷度數(shù)4.706). 隨通道二入口節(jié)流系數(shù)增大，系統(tǒng)穩(wěn)定性逐漸增強. 其中通道二的穩(wěn)定性增長較為顯著，盡管通道一的入口節(jié)流系數(shù)不變，其穩(wěn)定性也有較小的增長. 說明增大任一通道的入口節(jié)流系數(shù)總是有利的，除自身外，其他通道與整個系統(tǒng)的穩(wěn)定性都會增強.

圖15 增大單個通道的入口節(jié)流系數(shù)對系統(tǒng)和通道穩(wěn)定性的影響

Fig15. Effect of increasing inlet throttling coefficient of single channel on stability of system and each channel

3 結(jié) 論

1)軸向熱流不均勻分布時系統(tǒng)的穩(wěn)定性取決于沸騰邊界的位置以及積分熱流沿通道的分布.

2)與軸向熱流均勻分布相比，軸向熱流線性增大使系統(tǒng)穩(wěn)定性增強而軸向熱流線性減小使系統(tǒng)穩(wěn)定性減弱；入口過冷度較低時，軸向熱流余弦分布提高系統(tǒng)穩(wěn)定性，入口過冷度較高時，系統(tǒng)穩(wěn)定性可能增強也可能減弱. 與軸向熱流余弦分布相比，峰值偏向入口會降低系統(tǒng)穩(wěn)定性;而峰值偏向出口會提高系統(tǒng)穩(wěn)定性.

3)徑向不均勻加熱和入口節(jié)流時通道內(nèi)流量發(fā)生不均勻分配，系統(tǒng)的流動不穩(wěn)定由相對不穩(wěn)定的通道即功率較大或入口節(jié)流系數(shù)較小的通道導(dǎo)致.

4)徑向不均勻加熱對系統(tǒng)穩(wěn)定性影響較小. 對于徑向不均勻入口節(jié)流，系統(tǒng)穩(wěn)定性隨不均勻性的增大有減小的趨勢，可采用增大任一通道入口節(jié)流系數(shù)的方法提高系統(tǒng)的穩(wěn)定性.