基于COMSOL的板狀燃料組件流道堵塞事故研究

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

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

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

對于窄縫結(jié)構(gòu)內(nèi)部溫度場和流場的研究，由于測量技術(shù)限制，實驗實施困難，相關(guān)研究主要采取數(shù)值模擬的方法[2].計算流體動力學(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)象進行模擬，研究堵塞程度導致的堆型流場變化和溫度變化.

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

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個通道的半部分進行建模.

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

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

1.2 數(shù)學模型

假設(shè)冷卻劑為牛頓流體，在流道中的流動為無相變湍流流動.根據(jù)Daxin Gong[9]對不同湍流模型在板狀燃料組件計算適用性的研究，選取realizablek-ε兩方程模型作為求解湍流流場的數(shù)學模型.該模型是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)移和導出效率降低，熱量滯留在Ch1流道內(nèi)導致冷卻劑溫升明顯.左側(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ā)生相鄰兩流道均被堵塞的情況，整塊燃料板釋放的熱量將無法順利導出，溫度升高程度更嚴重，可能對組件結(jié)構(gòu)完整性產(chǎn)生較大威脅.