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

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

(1.三峽大學(xué) 機(jī)械與動(dòng)力學(xué)院，湖北 宜昌 443002；2.湖北省水電機(jī)械設(shè)備設(shè)計(jì)與維護(hù)重點(diǎn)實(shí)驗(yàn)室，湖北 宜昌 443002)

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

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

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

1 數(shù)值模型

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

表1 MTR堆芯參數(shù)

1.1 幾何描述

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

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

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

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

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

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

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

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

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

湍流動(dòng)能方程：

Gk+Gb-ρε-YM-Sk

(1)

擴(kuò)散方程：

(2)

其中：

(3)

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

1.3 網(wǎng)格劃分

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

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

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

1.4 邊界條件

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

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

2 結(jié)果分析

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

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

圖4 正常運(yùn)行時(shí)冷卻劑在不同截面的流速分布

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié) 論

以IAEA 10 MW MTR堆為對(duì)象，采用COMSOL程序?qū)ζ錁?biāo)準(zhǔn)燃料組件部分結(jié)構(gòu)進(jìn)行三維建模并求解了堆芯發(fā)生流道堵塞時(shí)熱工水力參數(shù)變化情況，得出以下結(jié)論：

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

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

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