鉛鉍冷卻燃料棒束堵流事故CFD模擬與分析

樊亦江，余大利，劉書(shū)勇，郁杰

（1. 中國(guó)科學(xué)院合肥物質(zhì)科學(xué)研究院核能安全技術(shù)研究所 合肥 230031；2. 中國(guó)科學(xué)技術(shù)大學(xué)合肥 230027；3. 先進(jìn)反應(yīng)堆工程與安全教育部重點(diǎn)實(shí)驗(yàn)室，清華大學(xué) 北京 100084）

鉛鉍冷卻快堆（Lead-bismuth cooled Fast Reactor，LFR）具備固有安全性、較高的能量密度和較長(zhǎng)的燃料循環(huán)壽期等特點(diǎn)，是第四代核反應(yīng)堆研究的重點(diǎn)堆型，這得益于液態(tài)鉛鉍合金（Lead-Bismuth Eutectic，LBE）優(yōu)異的中子學(xué)特性、物理特性、傳熱特性和高熱容量［1］。由于鉛與鉛秘合金在高溫和高速流動(dòng)下會(huì)對(duì)結(jié)構(gòu)材料產(chǎn)生較強(qiáng)的腐蝕作用［2］，堆芯沉積的腐蝕產(chǎn)物容易引起堵流事故，另外繞絲斷裂脫落、燃料元件的熱膨脹和輻射腫脹也會(huì)導(dǎo)致堵流事故的發(fā)生［3］。發(fā)生堵流事故時(shí)，堵塊附近冷卻劑的流通面積減少，堵塊后方出現(xiàn)回流區(qū)，導(dǎo)致流動(dòng)傳熱惡化，局部溫度顯著升高，威脅包殼完整性。

由于液態(tài)金屬的分子普朗特?cái)?shù)很低，相對(duì)于一般流體具有更厚的溫度邊界層［4］，不能通過(guò)相似原理采用常規(guī)介質(zhì)進(jìn)行?；瘜?shí)驗(yàn)，而液態(tài)鉛鉍具有的高溫、不透明、腐蝕性強(qiáng)等特點(diǎn)也對(duì)開(kāi)展實(shí)驗(yàn)造成了困難［5］。國(guó)內(nèi)外針對(duì)液態(tài)鉛鉍在棒束通道內(nèi)流動(dòng)進(jìn)行的實(shí)驗(yàn)較少，德國(guó)的Pacio 等人［6，7］通過(guò)實(shí)驗(yàn)研究了19 棒束含繞絲燃料組件的壓降和傳熱特性，并在不同子通道設(shè)置實(shí)心堵塊探究了組件局部堵塞產(chǎn)生的熱工水力效應(yīng)。意大利的Martelli 等人［8］開(kāi)展了基于格架固定的37 棒束燃料組件的熱工水力實(shí)驗(yàn)。許多學(xué)者也開(kāi)展了相關(guān)的模擬研究，陳寶文等人［9］使用Fluent 對(duì)不同堵塊參數(shù)下的軸向包殼溫度分布進(jìn)行了分析。Chai 等人［10］使用STAR-CCM+針對(duì)19 棒束燃料組件中的不同堵流工況進(jìn)行了模擬，分析了堵塞面積對(duì)包殼最大溫升和子通道間橫流速度的影響。孫暢［11］使用STAR-CCM+對(duì)19 棒束燃料組件多組堵流工況進(jìn)行了模擬，分析了堵塊參數(shù)對(duì)回流區(qū)特性的影響。Wu 等人［12］開(kāi)發(fā)了適用于含繞絲燃料棒束組件的子通道代碼SACOSPB，結(jié)果與19 棒束和61 棒束燃料組件實(shí)驗(yàn)符合良好。趙鵬程等人［13］針對(duì)100 MW 小型模塊化自然循環(huán)鉛冷快堆SNCLFR-100，基于ATHLETMOD3.0 A 系統(tǒng)程序開(kāi)展了堆芯功率最大組件局部堵流事故瞬態(tài)分析。

本文基于Pacio 等開(kāi)展的含繞絲19 棒束組件實(shí)驗(yàn)，使用Fluent對(duì)其一組正常工況進(jìn)行模擬，通過(guò)與實(shí)驗(yàn)值對(duì)比作為模型驗(yàn)證，進(jìn)而在子通道內(nèi)設(shè)置不同參數(shù)的堵塊，基于堵塞位置的速度場(chǎng)、渦結(jié)構(gòu)以及包殼表面溫度的發(fā)展，分析堵塊的厚度和面積對(duì)組件堵塞位置流動(dòng)傳熱的影響。

1 正常工況模擬及驗(yàn)證

1.1 幾何模型

本文模擬基于Pacio 等人設(shè)計(jì)的19 棒束實(shí)驗(yàn)組件，組件由正六邊形的外套管和19 根繞絲纏繞的加熱棒組成，棒束和子通道布置及編號(hào)如圖1 所示，組件幾何參數(shù)列于表1。由于繞絲與加熱棒的接觸方式為線接觸，不利于劃分高質(zhì)量的網(wǎng)格，建立幾何模型時(shí)通常將繞絲嵌入包殼內(nèi)0.1 mm 使二者由線接觸轉(zhuǎn)變?yōu)槊娼佑|［14-16］。此種做法提升了接觸位置的網(wǎng)格質(zhì)量，同時(shí)由于繞絲嵌入深度很小，對(duì)工質(zhì)流動(dòng)傳熱的影響可以忽略。

圖1 棒束和子通道編號(hào)示意圖Fig.1 Diagram of number of rod and sub-channel

表1 燃料組件幾何參數(shù)falseTable 1 Main geometrical parameters of fuel assembly

1.2 材料物性

模擬采用的工質(zhì)是液態(tài)鉛鉍合金，其物性隨溫度變化的經(jīng)驗(yàn)式參考Jaeger 等人［17］的報(bào)告，公式列于表2。

表2 液態(tài)鉛鉍合金物性參數(shù)falseTable 2 Liquid lead-bismuth eutectic material properties

1.3 網(wǎng)格劃分及網(wǎng)格無(wú)關(guān)性檢驗(yàn)

網(wǎng)格采用Fluent Meshing 軟件對(duì)流體域進(jìn)行多面體網(wǎng)格劃分，該軟件對(duì)復(fù)雜幾何區(qū)域的包容性較好，流體域全長(zhǎng)1200 mm，其中設(shè)置了300 mm 的入口段區(qū)域，網(wǎng)格如圖2 所示。基于最小單元尺寸為0.45 m、0.5 m、0.55 mm 的標(biāo)準(zhǔn)劃分出Mesh1、Mesh2、Mesh3 三種不同尺寸網(wǎng)格，網(wǎng)格數(shù)量分別為1537 w、1364 w、1255 w。由于1 號(hào)加熱棒受到的加熱功率最高，其表面包殼溫升普遍高于其他位置棒束，故提取不同網(wǎng)格數(shù)量下1 號(hào)加熱棒表面溫度在軸向上的發(fā)展進(jìn)行網(wǎng)格敏感性分析，如圖3 所示。三種網(wǎng)格模擬出的參考位置溫度值隨高度變化比較接近，綜合考慮計(jì)算精度和計(jì)算資源的消耗，最終選擇Mesh2 的網(wǎng)格劃分標(biāo)準(zhǔn)，并使邊界層網(wǎng)格Y+～1。

圖2 流體區(qū)域網(wǎng)格Fig.2 Mesh around fluid region

圖3 網(wǎng)格敏感性分析Fig.3 Mesh sensitivity analysis

1.4 湍流數(shù)值模型

在常用的湍流模型中，k-ε模型可以較好地模擬湍流充分發(fā)展區(qū)域的流動(dòng)，k-ω 模型在模擬附著邊界層湍流和分離湍流時(shí)則更為精確，SSTk-ω 模型在主流和近壁面處分別應(yīng)用 了k-ε和k-ω 模 式，Chai Xiang 等 人［10］也基于含繞絲棒束組件驗(yàn)證了SSTk-ω 模型的準(zhǔn)確性，因此本文選擇SSTk-ω 湍流模型進(jìn)行模擬。

湍流普朗特?cái)?shù)是湍流動(dòng)量擴(kuò)散率和湍流能量擴(kuò)散率之比，是為了?；牧鳠嵬慷氲臒o(wú)量綱數(shù)，對(duì)于水、空氣等流體常取0.85，對(duì)于液態(tài)鉛鉍合金等低普朗特?cái)?shù)流體則需使用經(jīng)驗(yàn)式進(jìn)行修正，Kays 模型基于當(dāng)?shù)氐耐牧鲄?shù)計(jì)算，并且與LES 結(jié)果符合較好［4］，本文即選用Kays 模型進(jìn)行計(jì)算，見(jiàn)式（1），其中α為熱擴(kuò)散率，vt為湍流黏度，單位均為m2·s-1。

采用SIMPLEC 算法進(jìn)行求解，離散格式選用二階迎風(fēng)格式，殘差降為10-4作為收斂標(biāo)準(zhǔn)。

1.5 邊界條件

在計(jì)算過(guò)程中，在繞絲和棒束表面設(shè)置均勻熱流密度邊界條件，使總功率與實(shí)驗(yàn)值197 kW 相等，其余壁面均為絕熱條件。入口質(zhì)量流量設(shè)為19.18 kg/s，入口溫度設(shè)為473.15 K。出口條件設(shè)為壓力出口。

1.6 結(jié)果分析

Pacio 實(shí)驗(yàn)在加熱段不同高度截面上設(shè)置了若干熱電偶，用以測(cè)量加熱棒表面溫度和典型子通道中心溫度，模擬值與實(shí)驗(yàn)值對(duì)比如圖4 所示，不同高度下的平均相對(duì)誤差分別為1.18%、2.57%、3.05%，模擬結(jié)果較為可信。

圖4 不同高度下測(cè)點(diǎn)溫度與模擬結(jié)果對(duì)比圖Fig.4 Comparison of temperature at different heights between experimental data and simulation results

2 堵流工況模擬

2.1 堵流工況設(shè)置

堵塊在燃料組件中的存在形態(tài)和材料特性目前尚無(wú)定論，通常的做法是設(shè)定一定厚度的堵塊填充關(guān)注的子通道位置，進(jìn)而分析組件局部堵塞對(duì)流動(dòng)傳熱的影響。本文選用與實(shí)驗(yàn)包殼材料相同的不銹鋼作為堵塊介質(zhì)，進(jìn)行不同厚度、面積堵塊的工況模擬，工況設(shè)置見(jiàn)表3。

表3 堵流工況設(shè)置Table 3 Blockage accidents parameters

2.2 堵塊厚度對(duì)流動(dòng)傳熱的影響

選取Case 0、Case 1、Case 2 分析堵塊厚度對(duì)流動(dòng)傳熱的影響，在加熱段每間隔1 mm 提取出1 號(hào)燃料棒包殼表面最高溫度，如圖5 所示，堵塊所在位置以虛線表示。可見(jiàn)包殼最大溫升出現(xiàn)在堵塊中心附近，隨堵塊厚度增加有小幅度上升。在單一子通道堵塞時(shí)，堵塊厚度增加并不會(huì)顯著增加包殼表面峰值溫度。

圖5 堵塊厚度對(duì)包殼最高溫度軸向分布的影響Fig.5 Impacts of the thickness of the block on the axial peak temperature distribution of cladding

由于繞絲作用，包殼表面最大溫度在軸向方向上出現(xiàn)周期性的漲落變化。單一子通道堵塞對(duì)流動(dòng)的擾動(dòng)，會(huì)影響到堵塊下游約0.25 倍繞絲螺距后中心加熱棒的傳熱。

參照?qǐng)D1 在法向330 度方向過(guò)Case 1 堵塊中心作縱向截面，如圖6（a）所示。Case 1 中堵塊占據(jù)一個(gè)子通道的位置，來(lái)流流經(jīng)堵塊位置時(shí)，流量向周圍三個(gè)子通道重新分配，在堵塊邊緣的影響下發(fā)生強(qiáng)烈的邊界層分離，由于存在逆壓力梯度，側(cè)面形成剪切層和渦旋結(jié)構(gòu)，在堵塊后方形成回流區(qū)并逐漸過(guò)渡至正常流動(dòng)。

圖6（b）為基于Q 準(zhǔn)則［18］得到的瞬時(shí)渦場(chǎng)結(jié)構(gòu)，可較為清晰地識(shí)別出強(qiáng)度差距明顯的渦結(jié)構(gòu)，Q 的計(jì)算公式見(jiàn)式（2）。堵塊側(cè)面較大的渦旋結(jié)構(gòu)實(shí)際上加快了不同速度之間流體的輸運(yùn)，在堵塊側(cè)表面形成局部的高速流動(dòng)加強(qiáng)了換熱。周圍三個(gè)子通道的繞流流場(chǎng)在堵塊后方分別形成渦旋結(jié)構(gòu)并匯聚、貼附，在下游壁面進(jìn)行再附著，渦旋會(huì)在堵塊后方發(fā)展出軸向速度與主流方向相反的區(qū)域，隨著流動(dòng)發(fā)展，回流區(qū)逐漸向主流區(qū)過(guò)渡。

圖6 Case 1 堵塊附近速度云圖及瞬時(shí)渦場(chǎng)結(jié)構(gòu)Fig.6 Velocity cloud image and instantaneous vortex field structure around the block of Case 1

在相同位置對(duì)Case 2 繪制速度云圖及瞬時(shí)渦場(chǎng)結(jié)構(gòu)圖，如圖7 所示。與圖6 對(duì)比可知，堵塊厚度較小時(shí)，流動(dòng)未在堵塊側(cè)面完成再附著，在堵塊后方會(huì)形成高強(qiáng)度的渦旋可以增強(qiáng)回流區(qū)換熱。堵塊厚度增加后，部分流體在堵塊側(cè)表面提前完成再附著，堵塊后方渦旋強(qiáng)度大大降低，回流區(qū)整體流速降低，對(duì)下游影響范圍增加。圖7 （b）中回流區(qū)貼近壁面有一較明顯渦旋，是由于繞絲的交混效應(yīng)所產(chǎn)生。

圖7 Case 2 堵塊附近速度云圖及瞬時(shí)渦場(chǎng)結(jié)構(gòu)Fig.7 Velocity cloud image and instantaneous vortex field structure around block of Case 2

2.3 堵塊面積對(duì)流動(dòng)傳熱的影響

Case 1、Case 3、Case 4 中，堵塊分別占據(jù)1、3、6 個(gè)中心子通道，同樣提取出1 號(hào)燃料棒包殼表面最高溫度繪圖，如圖8 所示。可見(jiàn)單一子通道堵塞對(duì)下游流動(dòng)傳熱的影響并不顯著，隨著堵塊面積增加，堵塊位置的包殼熱點(diǎn)溫度顯著增加，堵塊上游出現(xiàn)了小幅溫升，這是由于堵塊對(duì)上游流動(dòng)的阻礙作用增加，靠近堵塊位置的工質(zhì)流速明顯降低。

圖8 堵塊面積對(duì)包殼最高溫度軸向分布的影響Fig.8 Impacts of the area of the block on the axial peak temperature distribution of rod cladding

堵塊下游包殼熱點(diǎn)溫度整體升高，由于Case 3 中子通道堵塞的不對(duì)稱性，堵塞下游一段距離內(nèi)，溫度波動(dòng)劇烈，但局限在堵塊后方約0.5倍繞絲螺距范圍內(nèi)，在加熱段末端逐漸穩(wěn)定。Case 4 中子通道的大面積堵塞對(duì)堵塊后方區(qū)域的影響范圍更大，峰值溫度的波動(dòng)周期更長(zhǎng)。

在相同位置過(guò)Case 3、Case 4 堵塊中心作縱向截面，如圖9、圖10 所示?？梢?jiàn)相鄰子通道堵塞面積增加后，堵塊后方渦旋結(jié)構(gòu)位置升高，回流區(qū)高度、面積增加，堵塊所影響的下游低速區(qū)范圍增加，從而引起更大范圍的傳熱惡化。同時(shí)由于繞絲只對(duì)流動(dòng)起導(dǎo)引作用，所以相鄰子通道流速降低，也會(huì)使繞絲的交混作用不能充分發(fā)揮。

圖9 Case 3 堵塊附近速度云圖及瞬時(shí)渦場(chǎng)結(jié)構(gòu)Fig.9 Velocity cloud image and instantaneous vortex field structure around the block of Case 3

圖10 Case 4 堵塊附近速度云圖及瞬時(shí)渦場(chǎng)結(jié)構(gòu)Fig.10 Velocity cloud image and instantaneous vortex field structure around the block of Case 4

2.4 堵塞對(duì)包殼表面周向溫度的影響

1 號(hào)燃料棒包殼表面溫度峰值位置處沿周向的溫度分布如圖11 所示。燃料棒周圍局部子通道的不均勻堵塞會(huì)使包殼表面出現(xiàn)100～200 K的溫差，產(chǎn)生的熱應(yīng)力可能導(dǎo)致包殼變形、功能性受損。

圖11 1 號(hào)加熱棒包殼表面溫度峰值位置處的周向溫度Fig.11 The circumferential temperature at the peak temperature position of no.1 fuel rod cladding surface

3 結(jié)論

本文對(duì)液態(tài)鉛鉍合金在含繞絲19 棒束燃料組件中的正常工況和堵流工況進(jìn)行了模擬與分析，將正常工況的模擬值與實(shí)驗(yàn)值對(duì)比，驗(yàn)證了模型的精確性。在此基礎(chǔ)上，進(jìn)一步分析了堵塊厚度和面積對(duì)燃料棒包殼溫度和堵塊附近流動(dòng)傳熱的影響規(guī)律。得出主要結(jié)論如下。

（1）堵流事故中，包殼溫度在堵塊中部位置出現(xiàn)局部熱點(diǎn)，威脅包殼完整性，堵塊后方出現(xiàn)的回流區(qū)使包殼表面的工質(zhì)流速降低，造成傳熱惡化。大面積堵塞會(huì)導(dǎo)致回流區(qū)高度和影響范圍顯著增加，對(duì)包殼冷卻的影響可延續(xù)至下游較遠(yuǎn)處。

（2）小面積堵塞對(duì)流動(dòng)的影響較小，即便堵塊厚度增加，包殼熱點(diǎn)溫度上升也不顯著，但局部的不均勻堵塞會(huì)加劇堵塞下游包殼溫度的波動(dòng)，并導(dǎo)致包殼表面溫度不均勻性增加，產(chǎn)生的熱應(yīng)力可能導(dǎo)致包殼變形。