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

(西安交通大學(xué) 動(dòng)力工程及多相流國(guó)家重點(diǎn)實(shí)驗(yàn)室，核科學(xué)與技術(shù)學(xué)院，陜西 西安 710049)

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

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

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

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

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

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

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

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

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

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

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

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

圖5 FFTF簡(jiǎn)化幾何模型Fig.5 Simplified geometry model of FFTF

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

1) 堆芯組件模型

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

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

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

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

(1)

式中，下標(biāo)i為堆芯內(nèi)組件序號(hào)，cal為計(jì)算值，des為設(shè)計(jì)值。

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

(2)

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

2) 物理模型

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

3) 邊界條件和初始條件

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

圖7 冷管道進(jìn)口瞬態(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 計(jì)算結(jié)果與分析

2.1 穩(wěn)態(tài)計(jì)算結(jié)果

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

2.2 失流事故計(jì)算結(jié)果

1) 熱池流動(dòng)換熱特性分析

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

圖10 組件流量計(jì)算值與設(shè)計(jì)值的相對(duì)誤差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 熱管段溫度對(duì)比Fig.12 Comparison for hot leg temperatures

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

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

圖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) 堆芯流動(dòng)換熱特性分析

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

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

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

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

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

3 結(jié)論

本文采用CFD方法針對(duì)FFTF的LOSWOS#13無(wú)保護(hù)失流事故進(jìn)行了900 s的瞬態(tài)模擬，并將模擬值與實(shí)驗(yàn)值進(jìn)行了對(duì)比，得出如下結(jié)論：

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

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