快中子通量實驗堆失流事故三維數(shù)值模擬

周 磊，張大林，劉雅鵬，王式保，王心安，王成龍，田文喜，秋穗正，蘇光輝

(西安交通大學 動力工程及多相流國家重點實驗室，核科學與技術學院，陜西 西安 710049)

為準確評估鈉冷快堆在各類事故工況下的瞬態(tài)響應以及各類參數(shù)能否滿足事故工況下的驗收準則，通常采用系統(tǒng)程序針對鈉冷快堆進行建模計算以獲取事故工況下鈉冷快堆的各種關鍵參數(shù)，如西安交通大學開發(fā)的THACS[1-2]、中國原子能科學研究院開發(fā)的FR-Sdaso[3]和FRSYS[4]、華北電力大學開發(fā)的SAC-CFR[5-6]，以及美國阿貢國家實驗室開發(fā)的SAS4A/SASSYS-1[7]、俄羅斯國家科學中心物理動力研究院開發(fā)的GRIF[8]、法國原子能委員會開發(fā)的CATHARE[9]等，這類系統(tǒng)程序通常采用集總參數(shù)法或一維分析方法，以及有限的二維和偽三維分析方法以分析各類工況下的各種現(xiàn)象，為保證計算結果的可靠性，通常需要大量的實驗數(shù)據(jù)用以提供各種經(jīng)驗關系式或物理模型。該類系統(tǒng)程序雖能滿足鈉冷快堆的設計和安全分析，但開發(fā)驗證成本高，適用范圍相對窄，同時在一定程度上忽略了反應堆系統(tǒng)局部現(xiàn)象，犧牲了核動力系統(tǒng)的經(jīng)濟性。伴隨著高性能計算機的發(fā)展，對于核動力系統(tǒng)中部分系統(tǒng)與部件采用的全三維模擬逐漸成為可能。相較于系統(tǒng)程序，三維CFD方法能有效捕捉反應堆復雜的三維流動換熱現(xiàn)象，識別潛在的安全影響因素，充分挖掘反應堆的經(jīng)濟性[10]。但不同于系統(tǒng)程序針對特定反應堆長期開展的大量驗證與確認工作以及采用偏保守的物理模型，三維CFD方法從物理模型上、數(shù)值算法上具有更加廣泛的適用范圍與適用對象，與真實的物理現(xiàn)象更為接近，但針對不同流態(tài)(如層流、湍流、過渡流)、不同工質(zhì)(如水、液態(tài)鈉、液態(tài)鉛、熔鹽等)、不同工況(如強迫循環(huán)、自然循環(huán)工況)、不同流動現(xiàn)象(如回流、射流、繞流等)、不同流型(如單相流、氣液兩相流、氣固兩相流等)，CFD方法仍缺乏普適的湍流模型和數(shù)值算法，特別是在針對反應堆全堆級別的事故工況下的瞬態(tài)模擬驗證結果仍有限。目前已有的全堆數(shù)值模擬[11-12]采用商業(yè)軟件針對實驗快堆的失給水工況和主泵卡死等事故工況進行全三維模擬，但建模時對堆芯組件結構進行了較大的簡化，無法評估鈉冷快堆重要的盒間流現(xiàn)象[13]在瞬態(tài)過程中對反應堆安全性能的影響，此外模擬結果缺乏足夠的實驗數(shù)據(jù)的驗證，另一種較通用的方法是采用跨維度耦合[14-15]的方式進行模擬，針對核心關注部分采用三維CFD方法進行模擬，而針對其余部分采用系統(tǒng)程序進行模擬，該種方法在一定程度上結合了系統(tǒng)程序與CFD方法的優(yōu)勢，保證了精度，但對于兩種不同程序的耦合又會帶來收斂問題，同時對于瞬態(tài)變化劇烈的工況，難以獲取較好的模擬結果。

因此，為確保CFD方法在針對鈉冷快堆的低Pr流體液態(tài)金屬鈉[16]，以及鈉冷快堆自然循環(huán)工況下[17]計算的準確性，本文利用快中子通量實驗堆(FFTF)進行失流事故實驗，并針對其建立相應的三維模型，其中冷池和熱池按照原型尺寸建模，堆芯組件簡化為多孔介質(zhì)，同時保留堆芯盒間幾何特性，針對無保護失流事故進行900 s的瞬態(tài)模擬，并將模擬結果與實驗數(shù)據(jù)進行對比，初步驗證CFD方法在針對鈉冷快堆自然循環(huán)工況模擬中的可行性與可靠性。

1 FFTF系統(tǒng)及失流事故簡介

FFTF[18]是由美國西屋公司設計建造的一座400 MW熱功率，采用氧化物燃料、液態(tài)金屬鈉冷卻的實驗堆。FFTF位于華盛頓漢福德研究中心，于1980年首次臨界，一直運行到1992年，并于1993年關閉，期間進行了大量實驗與測試，為鈉冷快堆的燃料、材料、部件積累了大量經(jīng)驗。1986年，為測試FFTF非能動安全特性，在該堆上進行了大量無保護瞬態(tài)實驗，其中包含了13次無緊急停堆的失流事故實驗(LOFWOS)，為液態(tài)金屬堆的程序驗證、鈉冷快堆的固有安全性提供了重要支撐。2018年，在IAEA的組織下，選取LOFWOS#13作為國際基準題，多個國家和組織參與了此次基準題計算，西安交通大學作為參與單位之一，開展了系統(tǒng)程序和三維瞬態(tài)計算。

圖1 冷卻系統(tǒng)示意圖Fig.1 Coolant system overview

FFTF是一座環(huán)路式鈉冷快堆，其典型環(huán)路及主要部件示于圖1。FFTF包含2個回路、3個環(huán)路，液態(tài)鈉從反應堆容器離開，分別進入3個熱腿，在經(jīng)過每條環(huán)路上的中間熱交換器時將熱量傳遞給二回路中無放射性的鈉，最后通過二回路的空氣熱交換器將熱量導入環(huán)境中，因此該堆并不直接用來發(fā)電。FFTF還采用了獨特的安全設計以保證反應堆的安全，如較為特殊的氣體膨脹組件，能減小無保護失流事故的事故后果，該組件插入到堆芯徑向反射區(qū)的靠近內(nèi)堆芯區(qū)，如圖2所示，組件內(nèi)部頂端通過不銹鋼插件限制了約0.028 3 m3的氬氣，造成頂端氬氣無法泄漏，底端與柵板聯(lián)箱連通，液鈉可進入組件。因此，組件的液鈉高度與入口鈉的壓力及氬氣溫度有關，額定工況下，由于主泵提供較大的揚程，使得組件內(nèi)液鈉的高度高于堆芯活性區(qū)，而當堆芯溫度升高或失去強迫循環(huán)后，液鈉的高度下降，堆芯徑向中子泄漏增加，為反應堆提供負反饋，從而保障反應堆的安全性。堆芯組件分布示于圖3，圖中DF為燃料組件，CR、SR為控制棒組件，GEM為氣體膨脹組件。圖3中具有特殊標記的組件PIOTA為裝有溫度測點的實驗組件，LOFWOS#13實驗開始前保持50%額定功率、100%額定流量、約7.5 h，使得整個系統(tǒng)達到熱平衡狀態(tài)。為保證實驗完整運行，反應堆安全系統(tǒng)在實驗前進行了調(diào)整，避免控制棒過早插入，對堆內(nèi)流量和通量反應堆關閉信號也進行了關閉或旁通。實驗開始的初始時刻，主回路的3個環(huán)路主泵同時關閉，轉入惰轉，約90 s后主泵惰轉結束，反應堆轉入自然循環(huán)，整個實驗過程中二回路泵仍保持額定運行狀態(tài)。在實驗進行過程中未進行其他動作與干預。實驗開始前反應堆穩(wěn)態(tài)參數(shù)列于表1。

圖2 氣體膨脹組件Fig.2 Gas expansion module

圖3 堆芯組件類型及分布Fig.3 Assembly type and distribution in core

表1 反應堆穩(wěn)態(tài)參數(shù)Table 1 Steady state parameter of reactor

1.1 網(wǎng)格模型和計算邊界

FFTF堆本體如圖4所示，堆內(nèi)構件眾多，結構復雜，尺寸跨度大，完全按照原始幾何建模劃分網(wǎng)格，將導致計算量巨大，現(xiàn)階段難以完成瞬態(tài)計算。因此，針對真實反應堆結構進行了一定程度的簡化。本文建模中忽略了一系列與主要流道不相關的幾何腔室與結構，如圖5所示，僅保留了與流動密切相關的幾何結構，如冷池、熱池、堆芯等關鍵結構，此外將各分隔板簡化為無厚度的面，對于位于柵板聯(lián)箱中的堆芯支持結構以及入口管腳，該結構復雜且數(shù)目眾多，主要起流量分配的作用，因此在幾何建模過程中忽略了該部分結構，并將該部分結構的壓降及流量分配的作用歸到燃料組件的入口處。堆芯的各類組件采用繞絲棒束組件，流道狹長，針對每盒組件進行了多孔介質(zhì)簡化，最終的網(wǎng)格模型如圖6所示，其中堆芯組件采用六面體網(wǎng)格，其余部分采用四面體網(wǎng)格，整體網(wǎng)格呈混合網(wǎng)格，各類型網(wǎng)格面連續(xù)一致。

圖4 FFTF堆本體Fig.4 Overview of reactor vessel

圖5 FFTF簡化幾何模型Fig.5 Simplified geometry model of FFTF

1.2 數(shù)學模型與邊界條件

1) 堆芯組件模型

由于對反應堆堆內(nèi)結構進行了一定程度的簡化，特別是堆芯中忽略了組件插入柵板聯(lián)箱中的組件分流結構，可能直接導致堆芯各類組件的流量分配與設計值產(chǎn)生較大差距；另外堆芯各類幾何結構復雜，形狀阻力難以確定，為保證簡化后流量分配與設計值保持一致，在穩(wěn)態(tài)下采用如下方法進行組件進口區(qū)域的多孔介質(zhì)參數(shù)調(diào)整。

圖6 堆本體網(wǎng)格模型Fig.6 Mesh model for reactor vessel

將堆芯內(nèi)各組件按照各組件的實際結構沿軸向標記為進口區(qū)、出口區(qū)、繞絲棒束區(qū)3個分區(qū)類型，在初始計算前進行如下參數(shù)設置：進口區(qū)黏性阻力系數(shù)初始指定為0，慣性阻力系數(shù)初始指定為0；出口區(qū)黏性阻力系數(shù)指定為0，慣性阻力系數(shù)指定為0，并在計算過程中保持不變；繞絲棒束區(qū)黏性阻力系數(shù)和慣性阻力系數(shù)通過經(jīng)驗關系式獲得，并在計算過程中保持不變。

在穩(wěn)態(tài)調(diào)試中，每隔n個迭代步獲取堆芯每盒組件的出口流量，并根據(jù)式(1)計算堆芯內(nèi)各組件出口流量W的設計值與計算值之間的相對誤差εi。

(1)

式中，下標i為堆芯內(nèi)組件序號，cal為計算值，des為設計值。

根據(jù)式(1)計算堆芯內(nèi)各組件出口流量的設計值與計算值之間的相對誤差εi，通過式(2)調(diào)整堆芯各組件進口區(qū)多孔介質(zhì)的黏性阻力系數(shù)和慣性阻力系數(shù)，直至各組件的流量與設計流量的誤差在允許范圍，認為穩(wěn)態(tài)調(diào)試完成，可開展瞬態(tài)計算。

(2)

式中：1/α為黏性阻力系數(shù)；C2為慣性阻力系數(shù)；下標n為當前迭代步，m為上一迭代步，u為用戶指定修正的參數(shù)。

2) 物理模型

本文湍流模型采用Realizablek-ε模型、標準壁面函數(shù)、SIMPLE算法，梯度離散格式為Green-Gauss Node Based格式，壓力離散格式為二階格式，動量、湍動能、湍流耗散率和能量方程采用二階迎風格式，時間離散格式采用一階顯式格式，由于密度變化不大，浮升力作用采用Boussinesq假設。

3) 邊界條件和初始條件

計算時，進口采用流量進口邊界和進口溫度邊界，出口采用壓力邊界。其中進口溫度邊界采用圖7所示的瞬態(tài)溫度變化，該瞬態(tài)曲線根據(jù)實際反應堆實驗數(shù)據(jù)獲得。進口流量變化(歸一化值)如圖8所示，該曲線由實驗所得，該值乘以穩(wěn)態(tài)流量即為瞬態(tài)流量入口邊界。堆芯各類組件的穩(wěn)態(tài)功率分布如圖9所示，而瞬態(tài)下堆芯功率變化曲線(歸一化值)如圖8所示，該值乘以穩(wěn)態(tài)功率即為瞬態(tài)功率變化值。需要指出的是，在瞬態(tài)計算過程假設堆芯各組件功率份額與穩(wěn)態(tài)下保持一致。

圖7 冷管道進口瞬態(tài)溫度變化Fig.7 Transient temperature change at inlet of cold leg

圖8 歸一化堆芯功率和入口流量Fig.8 Normalized of core power and cold leg inlet mass flow rate

圖9 穩(wěn)態(tài)堆芯功率分布Fig.9 Core power distribution under steady condition

2 計算結果與分析

2.1 穩(wěn)態(tài)計算結果

每盒組件經(jīng)過調(diào)整多孔介質(zhì)參數(shù)后流量的計算值與設計值的相對誤差示于圖10，每盒組件的流量的相對誤差均小于±2%，且絕大多數(shù)組件流量的相對誤差小于±1%，因此可認為經(jīng)過調(diào)整后的多孔介質(zhì)較為合理，能滿足堆芯的流量分配需求，堆芯穩(wěn)態(tài)計算結果是合理的。穩(wěn)態(tài)下堆芯各類組件的出口溫度統(tǒng)計示于圖11，可看出，堆芯組件的最高出口溫度約為699 K，最低出口溫度約為602 K。

2.2 失流事故計算結果

1) 熱池流動換熱特性分析

反應堆堆芯出口溫度即熱管段溫度的模擬值與實驗值的對比示于圖12，可看出，鈉池巨大的鈉裝量使得在整個瞬態(tài)變化過程中溫度實驗值變化很小，前200 s內(nèi)基本保持不變，而后700 s緩慢下降，且在整個瞬態(tài)過程僅有4 K的溫度下降。CFD模擬值較好地捕捉了反應堆出口溫度的變化趨勢。

圖10 組件流量計算值與設計值的相對誤差Fig.10 Relative error between calculated value and design value of assembly mass flow rate

圖11 穩(wěn)態(tài)下堆芯組件出口溫度Fig.11 Outlet temperature of assembly under steady condition

圖12 熱管段溫度對比Fig.12 Comparison for hot leg temperatures

堆本體xz截面的溫度分布云圖示于圖13。從圖13可看出，初始階段，由于堆芯組件流速較高，致使在熱池中心區(qū)域形成1個高溫區(qū)域，而熱池外圍溫度相對較低，當堆芯流速降低后，中心高溫區(qū)域逐漸萎縮，從熱池頂部至熱池底部逐漸形成明顯的熱分層現(xiàn)象，特別是在熱池底部，由于流速低，從堆芯出來的冷流體可能積聚于此，直接導致此處的熱分層最為明顯。而對于冷池與其他腔室，整體溫度變化不大。

熱池xz截面不同時刻的速度矢量圖和速度云圖示于圖14。在初始階段，堆芯出口的流體向上流到液面處，在熱池中心區(qū)域形成1個高速區(qū)，而在中心高速區(qū)外，形成了2個大的漩渦，而后從熱管流出，同時伴隨著流速的下降，中心高速區(qū)域逐漸減小，特別是在后期，流速較低時，熱池上部的流體幾乎保持靜止，流體在幾乎與出口相同高度位置向外流出反應堆。

圖13 xz截面溫度云圖Fig.13 Temperature contour of xz cross section

圖14 熱池xz截面的速度矢量圖和速度云圖Fig.14 Vector and contour of velocity for xz cross section of hot pool

2) 堆芯流動換熱特性分析

具有溫度測點的組件的出口溫度模擬值與實驗值的對比示于圖15，該組件在堆芯中的位置如圖3所示。對于Row 6組件的溫度預測，實驗值除出現(xiàn)2個峰值溫度外，在100～150 s左右，溫度出現(xiàn)1個較為平緩的趨勢，而后在150 s后又開始下降，模擬值沒有體現(xiàn)該溫度變化趨勢，但考慮到Row 6組件在堆芯中的位置較特殊，周圍為氣體膨脹組件和非燃料組件，是堆內(nèi)熱工和物理條件最為復雜的組件，從熱工角度來說，外側非燃料組件溫度較低，而內(nèi)側為燃料組件，溫度較高，此處溫度梯度較大，盒間的導熱或?qū)α鲹Q熱作用較強烈；從物理角度來說，瞬態(tài)下，由于氣體膨脹組件失去強迫循環(huán)，氣體膨脹組件液位下降，堆芯中子泄漏增加，可能導致該組件的功率份額分布與穩(wěn)態(tài)分布差距較大，屬于物理熱工強耦合作用，單獨考慮熱工影響，可能不能準確模擬該過程。

圖15 Row 6 PIOTA組件出口溫度Fig.15 Outlet temperature of Row 6 PIOTA

堆芯活性出口截面的溫度分布示于圖16?？煽闯觯研咀顭峤M件的位置在堆芯中不斷發(fā)生變化，剛開始位于堆芯中心區(qū)域，而后靠近燃料組件的外圍區(qū)域，最終呈現(xiàn)出中心低、中間環(huán)型區(qū)域高、外側低的馬鞍形溫度分布，而從堆芯y=0截面的溫度分布曲線可看出，剛開始時堆芯功率較高，y=0截面溫度分布差距較大，而伴隨瞬態(tài)進程，由于盒間的存在，溫度分布逐漸被展平。

不同時刻下盒間出口的流動特性示于圖17。在初始階段，從外側非燃料組件區(qū)域向下流入盒間，而從中心較高溫區(qū)的燃料區(qū)域流出盒間，而在后期，由于堆芯溫度分布產(chǎn)生了較大變化，使得盒間從外部和中心區(qū)域流入堆芯，而從外側燃料區(qū)域流出堆芯。通過對比盒間的矢量圖和堆芯xy截面的溫度分布云圖，可看出矢量圖流動方向基本與溫度保持一致，都是從低溫區(qū)流入盒間，從高溫區(qū)流出盒間，這是由于低溫時密度較大，而高溫時密度較小，在浮升力的作用下，呈現(xiàn)這種流動趨勢。

圖17 盒間出口截面流動特性Fig.17 Flow characteristic of outlet for inter wrapper

3 結論

本文采用CFD方法針對FFTF的LOSWOS#13無保護失流事故進行了900 s的瞬態(tài)模擬，并將模擬值與實驗值進行了對比，得出如下結論：

1) CFD方法能較好地捕捉鈉冷快堆全堆的熱池內(nèi)熱分層等復雜的三維流動換熱特性，特別是對熱池巨大的熱容及其在瞬態(tài)下的熱惰性響應；

2) 堆芯最熱組件在瞬態(tài)過程中的位置不斷發(fā)生變化，堆芯盒間流的流動方向受盒內(nèi)溫度影響較大，且盒間在堆芯流速低、功率小的情形下，能有效展平堆芯溫度。