袁顯寶 劉 超 譚 偉 張永紅 張彬航 毛璋亮 周建軍 唐海波 劉芙蓉

(1.三峽大學(xué) 機械與動力學(xué)院，湖北 宜昌 443002；2.湖北省水電機械設(shè)備設(shè)計與維護重點實驗室，湖北 宜昌 443002)

堆芯流道堵塞是嚴重的堆芯事故之一.由于板狀燃料組件的設(shè)計緊湊，當(dāng)發(fā)生燃料腫脹、堆內(nèi)材料碎片或者異物進入堆芯循環(huán)等異常事件時，就可能引發(fā)堆芯堵流事故.此時，反應(yīng)堆內(nèi)冷卻劑的流動受到阻礙導(dǎo)致傳熱惡化，燃料釋放熱量的滯留將引起組件溫度升高并威脅燃料包殼的完整性，造成安全事故，導(dǎo)致反應(yīng)堆停止運行[1].

對于窄縫結(jié)構(gòu)內(nèi)部溫度場和流場的研究，由于測量技術(shù)限制，實驗實施困難，相關(guān)研究主要采取數(shù)值模擬的方法[2].計算流體動力學(xué)(computational fluicl dynamics, CFD)方法可構(gòu)建精細的三維幾何結(jié)構(gòu)，直接求解研究對象的能量守恒方程；并對堵塞區(qū)域進行合適的模擬，得到清晰的可視化結(jié)果，在反應(yīng)堆熱工水力研究中得到了廣泛應(yīng)用.宋磊[3]等使用Fluent軟件對板狀燃料組件進行了三維CFD計算，求解了單流道95%堵塞和全部堵塞造成的燃料組件內(nèi)的熱工水力變化.董化平、樊文遠[4-5]等利用Fluent研究了板型燃料組件和多層環(huán)形板狀燃料組件在發(fā)生入口堵塞后物理場的變化情況，探索了使用動網(wǎng)格技術(shù)模擬堵塞面的方法在研究該類事故中的應(yīng)用效果.Amgad Salama[6]對典型材料測試堆板狀燃料組件發(fā)生彎曲堵流現(xiàn)象進行模擬，研究堵塞程度導(dǎo)致的堆型流場變化和溫度變化.

目前國內(nèi)對板狀燃料組件發(fā)生堵塞的研究相對較少，本文基于COMSOL Multiphysics有限元軟件的熱流耦合技術(shù)，開展了板狀燃料組件單流道堵塞工況下的流體流動和流固傳熱數(shù)值研究，對板狀燃料組件流道堵塞事故的預(yù)防和事故嚴重性評估具有一定的參考價值.

1 數(shù)值模型

計算模型基于國際原子能機構(gòu)IAEA(International Atomic Energy Agency)10 MW 輕水冷卻和慢化的理想化池式材料測試堆(MTR，Material Test Reactor)，堆芯參數(shù)[3、6]見表1.

表1 MTR堆芯參數(shù)

1.1 幾何描述

MTR堆芯共包含21盒組件，其中1個標準燃料組件包含23塊燃料板，燃料板以平行等間距的方式排列.組件俯視圖如圖1所示.燃料板由燃料芯體和

圖1 燃料元件結(jié)構(gòu)(單位：mm)

外部覆蓋的包殼組成，冷卻劑沿著Z軸負方向通過板間縫隙并帶走燃料芯體所釋放的熱量.模型的上方和下方各設(shè)置高度為300 mm的空腔，使流體流動更符合堆芯運行時燃料組件周圍的冷卻環(huán)境[3-5].表2給出了模型幾何參數(shù).

表2 燃料元件的幾何參數(shù)

文獻調(diào)研結(jié)果表明，板狀燃料堵塞事故主要影響堵塞流道及相鄰的兩個流道[7]，故相關(guān)研究常選取組件邊緣起2到3塊燃料板及相鄰冷卻劑流道作為研究對象[1].考慮到板狀燃料布置和幾何的對稱性，選取了單組件邊緣2塊板和3個通道的半部分進行建模.

當(dāng)組件內(nèi)流道發(fā)生堵塞時，由于阻塞物體的形狀和厚度是不可預(yù)測的，根據(jù)典型堵流事故實例[8]，堵塞原因以堆外異物掉入為主，極有可能將整個流道入口完全覆蓋.依照樊文遠[4]提出的4種假定的入口堵塞形式，采取在流道入口設(shè)置剛性無厚度薄面的方式來模擬阻塞工況，屬于流道堵塞工況中的中心堵塞.圖2中左側(cè)圖形為模型的原本幾何形狀，考慮到燃料組件縱橫比較大，視覺上縮放為圖2右側(cè)容易辨識的三維幾何圖形.以下腔室一頂點為原點構(gòu)建三維模型，其俯視圖的a、b、c、d4個頂點與圖1平面圖中所選的研究區(qū)域相對應(yīng).

圖2 組件流道發(fā)生堵塞的幾何模型

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

假設(shè)冷卻劑為牛頓流體，在流道中的流動為無相變湍流流動.根據(jù)Daxin Gong[9]對不同湍流模型在板狀燃料組件計算適用性的研究，選取realizablek-ε兩方程模型作為求解湍流流場的數(shù)學(xué)模型.該模型是k-ε兩方程模型的修正形式，能對范圍較廣的湍流流動現(xiàn)象進行可信模擬，同時具有良好的魯棒性和計算經(jīng)濟性，通過引用兩個新的變量湍流動能k與湍流耗散率ε來求解湍流粘度μt.使用壁函數(shù)來近似替代近壁區(qū)域流體的法向速度變化.對于湍流區(qū)域，控制湍流動能k的輸運方程和湍流動能耗散率ε的輸運方程組為：

湍流動能方程：

Gk+Gb-ρε-YM-Sk

(1)

擴散方程：

(2)

其中：

(3)

式中：ρ為流體密度；μ為流體動力粘度；ν為流體的運動粘性系數(shù)；uj為流體在xj方向的流動速度；Gk為平均速度梯度面產(chǎn)生的湍動能；Gb為受浮升力而產(chǎn)生的湍動能；YM為可壓縮湍流中脈動膨脹對整體湍流擴散率的貢獻，在不可壓縮流體的計算中通常忽略掉；C1ε為常量，計算中取1.44；C2為常量，計算中取1.9；C3ε為浮力對ε的影響系數(shù)；σk為k的普朗特數(shù)，計算中取1.0；σε為ε的湍流普朗特數(shù)，計算中取 1.2；Sk為源項；Sε為源項.

1.3 網(wǎng)格劃分

由于板狀燃料間隙狹窄，僅有2.23 mm，燃料板厚度僅有1.27 mm，而模型在高度方向整體長1 200 mm，縱橫比較大，對網(wǎng)格精度要求較高.冷卻劑在由上腔室進入流道和從流道中流出時，湍流強度高，速度變化劇烈，因此需要對相應(yīng)區(qū)域的網(wǎng)格加強細化.

參考Amgad Salama[6]對燃料板組件堵流事故進行研究時二維計算模型的網(wǎng)格劃分方式，采取結(jié)構(gòu)化網(wǎng)格和非結(jié)構(gòu)化網(wǎng)格相結(jié)合的方式對模型進行網(wǎng)格劃分，如圖3所示.對燃料板及流道內(nèi)冷卻劑區(qū)域進行劃分時采取結(jié)構(gòu)化網(wǎng)格，上下腔的流體區(qū)域采取非結(jié)構(gòu)化網(wǎng)格，邊界層設(shè)為3層.進行無關(guān)性驗證后，網(wǎng)格包含400萬個計算單元，整體平均質(zhì)量為0.896 6.

圖3 板狀燃料元件模型的網(wǎng)格劃分

1.4 邊界條件

設(shè)置湍流模塊和流體與固體傳熱模塊，湍流模塊中依次設(shè)置入口邊界條件為速度入口、出口邊界條件為0、壓力出口以及壁條件為壁函數(shù).為了簡化模型設(shè)置了兩個對稱邊界.采用內(nèi)壁條件模擬堵塞面，該邊界位于Ch1通道頂端，改變該邊界的面積占比來控制堵塞程度.對稱面、堵塞面位置詳見圖2.流體與固體傳熱模塊中分別選定燃料板為固體域、冷卻劑所在區(qū)域為流體域，流入溫度為311.15 K，出口邊界和對稱邊界與湍流模塊相一致.設(shè)置燃料芯體為熱源，釋熱功率按體積均勻分布.建模所需具體參數(shù)和初始值設(shè)定見表3.

表3 邊界條件和初始值[3]

2 結(jié)果分析

2.1 正常運行狀態(tài)

正常運行時，冷卻劑自頂向下流動，在燃料板頂面受到阻礙形成3條分流，分別流經(jīng)3個狹窄通道(Ch1、Ch2、Ch3)，在下腔室匯聚，最終從下腔室底面流出.圖4中云圖為模型中分別位于上腔室、燃料區(qū)、下腔室3個區(qū)域XY截面上冷卻劑的流速分布.

圖4 正常運行時冷卻劑在不同截面的流速分布

圖5為對應(yīng)圖4中截線Y=16.625 mm上冷卻劑流速在X軸的分布.冷卻劑以4.1 m/s速度進入上腔室，即將進入流道時，流速分布較平坦，右側(cè)的近壁區(qū)流速出現(xiàn)衰減；冷卻劑受燃料板阻隔后分流，在流道內(nèi)流速分布趨于一致，流速峰值在流道中心，達到7 m/s；經(jīng)流道流出后，右側(cè)近壁區(qū)的流速在流動方向上衰減緩慢，相對流速較高，該結(jié)果與參考文獻[7]一致.

圖5 組件不同高度截線上冷卻劑流速隨X軸變化

圖6所示云圖為XZ面上的溫度分布，位置為Y=16.625 mm.由于各流道冷卻劑流量大致相同，在各燃料板功率分布一致的條件下，組件內(nèi)部的溫度場分布基本對稱.燃料板主要以對流換熱方式向冷卻劑轉(zhuǎn)移自身產(chǎn)熱.在流動方向上，隨著冷卻劑被加熱，與燃料板溫差減小，換熱效率逐漸降低，冷卻作用逐漸減小，燃料板的溫度逐漸上升，接近底部時達到342 K；在X軸方向上，燃料板與冷卻劑溫差較大，包殼附近溫度梯度較大，但流道中心大部分區(qū)域溫度較低，冷卻劑對燃料板釋熱仍有一定裕量.

圖6 Y=16.625 mm面的溫度分布

2.2 堵塞事故狀態(tài)

堵塞會對冷卻劑的正常流動造成較大影響，引起流量重新分配.如圖7(a)所示，在流道入口附近，堵塞面的存在迫使原本將流入Ch1的冷卻劑向旁邊流道轉(zhuǎn)移.流道出口附近，堵塞流道Ch1到下腔室的過渡區(qū)域發(fā)生較為明顯的漩渦和回流現(xiàn)象，如圖7(b)所示.經(jīng)過對比各流道入口橫截面上的平均流速，得到表4關(guān)于不同運行工況下各冷卻劑流量占比.

圖7 堵塞引起流量轉(zhuǎn)移和回流

表4 3種工況下各流道流量占比 (單位：%)

堵塞發(fā)生后，流量再分配和流道內(nèi)漩渦的存在使得相應(yīng)區(qū)域的溫度發(fā)生明顯變化，如圖8所示.由圖可知，未發(fā)生堵塞時，兩燃料板溫差較小，Ch1內(nèi)冷卻劑平均溫度與入口溫度相比無明顯升高；發(fā)生95%堵塞后，Ch1流道的冷卻劑流動速度減慢，燃料板熱量轉(zhuǎn)移和導(dǎo)出效率降低，熱量滯留在Ch1流道內(nèi)導(dǎo)致冷卻劑溫升明顯.左側(cè)燃料板出現(xiàn)不連續(xù)的高溫區(qū)域，Ch2、Ch3流道及右側(cè)燃料板受影響較??；全堵塞時，左側(cè)燃料板溫度出現(xiàn)了連續(xù)的高溫區(qū)域，Ch1流道內(nèi)冷卻劑溫度大幅升高且與左側(cè)燃料板溫度十分接近，但右側(cè)燃料和Ch2、Ch3流道的溫度變化仍不明顯.

圖8 組件中軸線的溫度分布

圖9為不同工況下三維截線(面Y=16.625 mm與面Z=600 mm的交線)上的溫度分布.在正常工況下，流經(jīng)3個流道的冷卻劑流量基本相同，燃料板向兩側(cè)傳遞的熱量基本相同，因此溫度基本呈對稱分布，Ch3流道由于右側(cè)無熱源，整體溫度稍低于Ch1.95%堵塞工況下，Ch1流道的流量急劇下降，整個流道內(nèi)的溫度大幅升高，部分區(qū)域達到了正常運行時燃料板處溫度340 K.堵塞引起的流量降低和回流現(xiàn)象使溫度峰向左側(cè)小幅度偏移，但影響范圍僅限Ch1流道和左側(cè)燃料板，此時Ch2流道的溫度上升現(xiàn)象并不明顯.由于堵塞使流經(jīng)Ch1的流量向Ch2和Ch3轉(zhuǎn)移，使得右側(cè)燃料板與冷卻劑之間的換熱更加充分，溫度小幅降低.

圖9 組件中軸線的溫度分布

全部堵塞工況下，Ch1流道的冷卻劑幾乎無法提供冷卻作用，燃料板溫度進一步升高，達到了360 K.但同95%堵塞狀況相似，全堵塞工況對右側(cè)燃料板傳熱無明顯影響.中軸線上的溫度分布狀況與參考文獻[3]中Fluent計算所得結(jié)果(E)偏差較小.

表5給出了3種工況下，組件不同區(qū)域的最高溫度.全堵塞工況相對于95%堵塞工況的各部分區(qū)域最高溫度均有上升，包殼的最高溫度為362.74 K，冷卻劑最高溫度360.45 K，均未超過環(huán)境壓力0.17 MPa時的飽和溫度388 K[7]，維持在安全范圍內(nèi)，因此不會在燃料板包殼表面產(chǎn)生過熱沸騰現(xiàn)象.

表5 堵塞前后不同區(qū)域中的最高溫度 (單位：K)

3 結(jié) 論

以IAEA 10 MW MTR堆為對象，采用COMSOL程序?qū)ζ錁藴嗜剂辖M件部分結(jié)構(gòu)進行三維建模并求解了堆芯發(fā)生流道堵塞時熱工水力參數(shù)變化情況，得出以下結(jié)論：

1)穩(wěn)態(tài)運行時，組件在X方向溫度呈周期性對稱分布，沿冷卻劑的流動方向溫度逐漸升高，接近底端時到達溫度峰值，約為342.73 K.

2)發(fā)生95%堵塞時，受阻塞影響冷卻劑流量將重新分配，堵塞流道Ch1內(nèi)的冷卻劑將形成漩渦并影響傳熱.由于冷卻劑仍保持流動，傳熱惡化的程度有限.完全堵塞時，左側(cè)料板單側(cè)喪失冷卻，溫度迅速升高，燃料包殼的最高溫度達到了362.74 K，但未超過冷卻劑飽和溫度，不會引起沸騰.單流道堵塞影響區(qū)域僅限于堵塞流道及相鄰的燃料板，相鄰流道溫度升高不明顯.

3)完全堵塞時，堵塞流道內(nèi)的冷卻劑溫度可達360.45 K，考慮到每塊燃料板僅有兩側(cè)流道，流道間隙較狹窄，存在相鄰兩流道同時受到阻塞的可能性.若發(fā)生相鄰兩流道均被堵塞的情況，整塊燃料板釋放的熱量將無法順利導(dǎo)出，溫度升高程度更嚴重，可能對組件結(jié)構(gòu)完整性產(chǎn)生較大威脅.