孫 旭，周 琦,*，余慧鶯，朱慶福，夏兆東，寧 通，馬驍?shù)?/p>

(1.中國原子能科學(xué)研究院 反應(yīng)堆工程技術(shù)研究所，北京 102413；2.國家電投集團(tuán)科學(xué)技術(shù)研究院有限公司，北京 102209)

在核燃料循環(huán)領(lǐng)域，易裂變核材料在某些環(huán)節(jié)上是以溶液的形式存在，溶液中的水是中子的優(yōu)良慢化材料，使得溶液中具有較少的易裂變核材料便可達(dá)到臨界。同時(shí)溶液具有易流動(dòng)、易變形的特點(diǎn)，以及多相鈾钚混合體系、多體相互作用等使得核材料的臨界安全控制更加復(fù)雜[1]。據(jù)文獻(xiàn)報(bào)道，各國在核燃料循環(huán)領(lǐng)域共發(fā)生過22起臨界事故，其中有21起發(fā)生在溶液或漿液中[2]。臨界事故不僅會(huì)引起國際社會(huì)的廣泛關(guān)注，嚴(yán)重影響核能發(fā)展，還會(huì)造成大量放射性物質(zhì)的釋放、人員傷亡和環(huán)境污染。因此，有必要對(duì)溶液系統(tǒng)臨界事故進(jìn)程進(jìn)行研究，從而對(duì)臨界事故的預(yù)防、屏蔽設(shè)計(jì)及后果評(píng)價(jià)提供指導(dǎo)。

自2011年以來，中國原子能科學(xué)研究院一直承擔(dān)核燃料系統(tǒng)瞬發(fā)臨界事故分析過程與方法研究的任務(wù)，研制了包括GETAC-M(固體)、GETAC-S(溶液)、GETAC-WP(濕粉末)及GETAC-DP(干粉末)等一系列核燃料系統(tǒng)瞬態(tài)分析程序[3-4]。包括GETAC-S在內(nèi)的大多溶液系統(tǒng)瞬態(tài)分析程序使用點(diǎn)堆模型對(duì)事故進(jìn)程進(jìn)行估計(jì)，然而核燃料循環(huán)體系中溶液系統(tǒng)存在多種幾何形態(tài)，且存在多相溶液的不均勻混合等特征，點(diǎn)堆模型難以適應(yīng)復(fù)雜的計(jì)算條件。本文基于蒙特卡羅均勻化理論與有限體積方法，提出適用于瞬態(tài)分析問題的三維擴(kuò)散時(shí)空動(dòng)力學(xué)模型，將該模型與非穩(wěn)態(tài)傳熱模型、輻照裂解氣泡模型相耦合，對(duì)GETAC-S進(jìn)行升級(jí)。

1 基礎(chǔ)理論

核溶液系統(tǒng)發(fā)生臨界事故后，系統(tǒng)功率將在短時(shí)間內(nèi)迅速上升得到一個(gè)功率峰值，而后溫度反饋效應(yīng)和輻照裂解氣泡反饋效應(yīng)會(huì)使系統(tǒng)功率回到較低水平。隨溫度下降和氣泡消失，這些反饋效應(yīng)的減小和消失，系統(tǒng)功率會(huì)重新達(dá)到一個(gè)比第一個(gè)功率峰值低的值，如此反復(fù)形成功率震蕩，直至部分水被蒸發(fā)或易裂變材料濺出而使系統(tǒng)變?yōu)榇闻R界狀態(tài)停止功率震蕩[5]。為了模擬臨界事故的發(fā)展過程，需要從中子物理和熱工水力反饋兩方面進(jìn)行分析。物理上采用三維擴(kuò)散時(shí)空動(dòng)力學(xué)模型可為熱工計(jì)算提供更加準(zhǔn)確的功率分布，從而對(duì)事故進(jìn)程的發(fā)展給出更為精確的模擬結(jié)果。

三維擴(kuò)散時(shí)空動(dòng)力學(xué)的求解通常需要將系統(tǒng)均勻化近似后得到多群常數(shù)作為輸入條件，本文采用蒙特卡羅方法產(chǎn)生均勻化多群常數(shù)，并采用有限體積方法進(jìn)行三維擴(kuò)散時(shí)空動(dòng)力學(xué)的求解。

1.1 蒙特卡羅均勻化理論

蒙特卡羅方法一般使用連續(xù)能量截面，對(duì)三維實(shí)際問題建立精確的幾何模型，并考慮各向異性散射，相較于確定論輸運(yùn)計(jì)算方法，在能量、空間和角度等方面的近似更少，能夠準(zhǔn)確描述非均勻性幾何、邊界條件、共振等物理特性，是較為精確的輸運(yùn)方法。

蒙特卡羅均勻化的核心是將求解輸運(yùn)問題得到的各種事件統(tǒng)計(jì)結(jié)果轉(zhuǎn)化為均勻化常數(shù)，這些事件主要包括散射、裂變、俘獲吸收、(n，xn)等反應(yīng)，均勻化常數(shù)包括各種反應(yīng)的群截面、輸運(yùn)修正、散射矩陣、各階勒讓德項(xiàng)，次級(jí)中子產(chǎn)生矩陣、裂變產(chǎn)額、裂變能譜，以及動(dòng)力學(xué)計(jì)算需要的中子速度倒數(shù)、緩發(fā)中子先驅(qū)核產(chǎn)額與衰變常量等。

為簡(jiǎn)化表達(dá)式，引入內(nèi)積算符〈·，·〉來代表能量、空間或角度的積分形式，這樣第x類核反應(yīng)的反應(yīng)率計(jì)數(shù)估計(jì)可表示為：

(1)

式中：Σx為第x類核反應(yīng)的截面；φ為中子通量密度；V為體積；S為角度；E為中子能量；r為空間位置；Ω為中子運(yùn)動(dòng)的方向。

內(nèi)積函數(shù)的統(tǒng)一性，如空間均勻化與能量積分通量可表示為：

〈φ〉k,g≡〈φ,1〉k,g=

(2)

式中：k為空間網(wǎng)格編號(hào)；g為能群編號(hào)；Vk為第k網(wǎng)格體積；Eg為第g能群的能量上限。

(3)

1.2 多群時(shí)空中子動(dòng)力學(xué)擴(kuò)散方法

有限體積方法是一種處理偏微分方程的數(shù)值方法，在熱工、流體、力學(xué)等領(lǐng)域有著廣泛應(yīng)用，將有限體積方法用于中子擴(kuò)散方程的求解，能夠更加有效地完成中子物理參數(shù)與其他參數(shù)的交互，實(shí)現(xiàn)耦合與反饋計(jì)算，提高反應(yīng)堆物理分析的準(zhǔn)確性。

多群時(shí)空中子動(dòng)力學(xué)擴(kuò)散方程如式(4)所示，本文使用有限體積方法對(duì)其進(jìn)行離散求解[6]。

g=1,…,G

(4)

i=1,…,ND

(5)

式中：vg為第g群中子運(yùn)動(dòng)速度；φg為第g群中子通量密度；t為時(shí)間；Dg為第g群中子的擴(kuò)散系數(shù)；Σr，g為第g群中子的移出截面；Σs，g′→g為第g′群中子到第g群的散射截面；χg為第g群中子的裂變能譜；β為緩發(fā)中子份額；(νΣf)g′為第g′群中子的平均裂變中子數(shù)與裂變截面的乘積；i為緩發(fā)中子先驅(qū)核組編號(hào)，組數(shù)為ND；λi為第i組先驅(qū)核衰變常量；Ci為第i組先驅(qū)核濃度；G為能群數(shù)。

2 軟件實(shí)現(xiàn)

使用蒙特卡羅均勻化方法進(jìn)行多群常數(shù)的生成，然后采用有限體積方法進(jìn)行三維擴(kuò)散中子時(shí)空動(dòng)力學(xué)行為研究，采用非穩(wěn)態(tài)傳熱模型與輻照裂解氣泡模型進(jìn)行熱工反饋分析。

首先需要進(jìn)行穩(wěn)態(tài)中子擴(kuò)散方程的求解(臨界源問題)，得到初始狀態(tài)下求解區(qū)域內(nèi)的中子通量密度分布；然后進(jìn)行非穩(wěn)態(tài)擴(kuò)散方程的求解(固定源問題)，以得到隨時(shí)間變化的中子通量密度與緩發(fā)中子先驅(qū)核濃度，再通過將熱工反饋計(jì)算得到的溫度、空泡系數(shù)反饋到非穩(wěn)態(tài)中子擴(kuò)散方程的求解中，反復(fù)迭代，直到模擬時(shí)間結(jié)束。

2.1 基于蒙特卡羅方法的多群常數(shù)生成

蒙特卡羅程序可通過計(jì)數(shù)功能進(jìn)行多群常數(shù)的生成。這些多群截面可提供給各種能群結(jié)構(gòu)或空間網(wǎng)格劃分的確定論程序。從能量上看，蒙特卡羅程序可計(jì)算單獨(dú)的核素或元素的宏觀或微觀核反應(yīng)率，可得到一個(gè)或多個(gè)任意能群結(jié)構(gòu)的多群截面，便于研究能群結(jié)構(gòu)對(duì)多群截面精度的影響；從空間上看，蒙特卡羅程序還能對(duì)各種空間劃分方式進(jìn)行多群截面的計(jì)算，包括材料、柵元以及區(qū)域相關(guān)的笛卡爾坐標(biāo)系下的空間劃分，還支持重復(fù)結(jié)構(gòu)系統(tǒng)，如組件或堆芯的計(jì)算。利用蒙特卡羅方法生成多群常數(shù)，滿足軟件適用于復(fù)雜幾何、材料以及能譜結(jié)構(gòu)的目標(biāo)要求。采用OpenMC程序計(jì)算了多群截面、瞬發(fā)中子常數(shù)以及緩發(fā)中子常數(shù)等參數(shù)[7]。

2.2 基于OpenFOAM的三維擴(kuò)散時(shí)空動(dòng)力學(xué)程序

三維擴(kuò)散時(shí)空動(dòng)力學(xué)程序是在使用有限體積方法求解偏微分方程的OpenFOAM的基礎(chǔ)上進(jìn)行開發(fā)[8]。OpenFOAM是一個(gè)完全面向?qū)ο蟮拈_源CFD類庫，這些庫全部是由C++語言編寫，OpenFOAM針對(duì)偏微分方程具有多種離散求解方式，這些求解方法同樣適用于三維擴(kuò)散時(shí)空動(dòng)力學(xué)方程的求解。

基于OpenFOAM求解三維擴(kuò)散時(shí)空動(dòng)力學(xué)方程的基本流程如圖1所示。主要分為3部分，分別是軟件輸入、網(wǎng)格生成以及擴(kuò)散方程的求解。其中輸入主要分為6部分：neutronProperties、topoSetDict、controlDict、decomposeParDict、blockMeshDict及phi.group，用于初始條件及控制條件等參數(shù)的輸入。網(wǎng)格生成用于將整個(gè)計(jì)算空間離散為網(wǎng)格，每個(gè)網(wǎng)格內(nèi)選擇1個(gè)點(diǎn)作為此網(wǎng)格的代表。由于蒙特卡羅均勻化理論已將復(fù)雜的幾何模型了進(jìn)行簡(jiǎn)化，這里使用正方形網(wǎng)格劃分即可滿足要求。離散后的擴(kuò)散方程求解使用OpenFOAM提供的函數(shù)庫實(shí)現(xiàn)。

圖1 基于OpenFOAM求解時(shí)空動(dòng)力學(xué)方程的基本流程Fig.1 Basic process of solving spatiotemporal dynamics equation based on OpenFOAM

三維擴(kuò)散時(shí)空動(dòng)力學(xué)程序的輸出主要分為3部分，分別為keff、通量分布以及功率分布。瞬態(tài)過程的計(jì)算首先需要求解穩(wěn)態(tài)擴(kuò)散方程，將穩(wěn)態(tài)計(jì)算通量分布用于瞬態(tài)的計(jì)算條件。因此，內(nèi)核方程的求解分為兩種情形：不含時(shí)間項(xiàng)的穩(wěn)態(tài)中子擴(kuò)散方程的求解與含時(shí)間項(xiàng)的中子擴(kuò)散時(shí)空動(dòng)力學(xué)方程(式(4)、(5))的求解。

2.3 熱工反饋與耦合的實(shí)現(xiàn)

式(6)、(7)分別為熱工反饋采用非穩(wěn)態(tài)傳熱模型和輻照裂解氣泡模型。

(6)

式中：ρ為介質(zhì)密度，kg/m3；c為比熱容，J·kg-1·K-1；T為溫度，K；λ為介質(zhì)導(dǎo)熱系數(shù)，W·m-1·K-1；S為內(nèi)熱源，對(duì)核燃料溶液系統(tǒng)即為功率密度，W/m3；h為軸向上的坐標(biāo)；r為徑向上的坐標(biāo)。

(7)

式中：V(t，z)為t時(shí)刻、軸向位置z處的空泡體積份額；ζ為氫氣產(chǎn)量與其他氣體產(chǎn)量之比；G(H2)為單位能量的氫氣產(chǎn)額，mol/J；P(t，z)為t時(shí)刻、z位置處的功率密度，W/m3；R為摩爾氣體常數(shù)；Tg為氣泡內(nèi)部溫度，K；p0為環(huán)境大氣壓強(qiáng)，Pa；g為重力加速度，m2/s；H(z)為z位置上方溶液的高度，m；σ為溶液的表面張力系數(shù)，N/m；r為氣泡半徑，m；v(z)為z位置處的氣泡向上運(yùn)動(dòng)速度，m/s，參考CRITEX的兩速度模型，設(shè)定功率上升時(shí)，氣泡的上升速度為4.0×10-2m/s，功率下降時(shí)，氣泡的上升速度為1.5×10-2m/s。

基于OpenFOAM平臺(tái)，對(duì)上述兩項(xiàng)熱工反饋模型進(jìn)行離散求解[9]，并實(shí)現(xiàn)物理計(jì)算與熱工反饋計(jì)算的網(wǎng)格匹配。由于熱工反饋計(jì)算得到的參數(shù)為溫度與空泡份額，GETAC-S通過設(shè)置溫度與空泡反饋系數(shù)將其傳遞到反應(yīng)性的變化中。在升級(jí)的程序中，由于使用三維擴(kuò)散時(shí)空動(dòng)力學(xué)模型代替點(diǎn)堆模型，而多群截面、動(dòng)力學(xué)參數(shù)等作為擴(kuò)散軟件的輸入代替作為點(diǎn)堆模型的反應(yīng)性輸入。因此，使用原來的反饋系數(shù)傳遞反應(yīng)性的方法不能充分體現(xiàn)出三維擴(kuò)散時(shí)空動(dòng)力學(xué)的優(yōu)勢(shì)。

采用以溫度與空泡份額為自變量，將截面預(yù)先制成相關(guān)函數(shù)，形成函數(shù)庫(或擬合成曲線)，在每一時(shí)間步長(zhǎng)內(nèi)熱工反饋及計(jì)算結(jié)束后，得到相應(yīng)的溫度與空泡份額，然后調(diào)用相關(guān)函數(shù)庫(或擬合曲線)得到對(duì)應(yīng)溫度與空泡份額下的少群截面與動(dòng)力學(xué)參數(shù)，將其代入到時(shí)空動(dòng)力學(xué)方程中，實(shí)現(xiàn)每一時(shí)間步長(zhǎng)內(nèi)的迭代求解。

3 TRACY實(shí)驗(yàn)驗(yàn)證

使用日本TRACY瞬態(tài)實(shí)驗(yàn)裝置的模擬臨界事故實(shí)驗(yàn)數(shù)據(jù)對(duì)升級(jí)的GETAC-S程序進(jìn)行可靠性驗(yàn)證[10]。TRACY裝置的主體是一個(gè)帶有中間孔道的環(huán)形不銹鋼容器，如圖2所示。實(shí)驗(yàn)時(shí)可通過拔出中間孔道內(nèi)的碳化硼控制棒從而模擬不同反應(yīng)性引入所引起的瞬時(shí)功率上升事故。TRACY裝置可記錄瞬態(tài)事故發(fā)展過程中功率、溫度等參數(shù)的變化，本文對(duì)其中的RUN72實(shí)驗(yàn)[11-12]進(jìn)行了驗(yàn)證，表1列出RUN72實(shí)驗(yàn)?zāi)M參數(shù)，表2列出RUN72實(shí)驗(yàn)計(jì)算結(jié)果及其相對(duì)偏差。由表2可看到，程序計(jì)算的功率峰值、總能量和最終溫度與實(shí)驗(yàn)結(jié)果的相對(duì)偏差均在10%以內(nèi)，結(jié)果均符合良好，驗(yàn)證了升級(jí)后的GETAC-S在計(jì)算結(jié)果上有較大改進(jìn)(改進(jìn)前GETAC-S計(jì)算結(jié)果與實(shí)驗(yàn)的相對(duì)偏差均在10%以上)。圖3示出RUN72實(shí)驗(yàn)功率與溫度模擬結(jié)果，二者的程序計(jì)算值與RUN72實(shí)驗(yàn)值均符合較好。

圖2 TRACY裝置結(jié)構(gòu)[12]Fig.2 TRACY device structure[12]

表1 RUN72實(shí)驗(yàn)?zāi)M參數(shù)Table 1 RUN72 experiment simulation parameter

表2 RUN72實(shí)驗(yàn)計(jì)算結(jié)果及相對(duì)偏差Table 2 RUN72 experiment calculation result and relative deviation

圖3 RUN72實(shí)驗(yàn)功率與溫度模擬結(jié)果Fig.3 Power and temperature simulation results of RUN72 experiment

4 瞬態(tài)分析程序在事故反演上的應(yīng)用

核燃料瞬態(tài)分析程序的一個(gè)重要作用在于對(duì)已發(fā)生的事故進(jìn)行反演，了解事故發(fā)生的過程，從而對(duì)事故的預(yù)防、屏蔽設(shè)計(jì)等提供指導(dǎo)。本文利用改進(jìn)后的GETAC-S對(duì)日本JCO臨界事故進(jìn)行了反演[13]。JCO事故模擬參數(shù)列于表3[14-15]。

表3 JCO事故模擬參數(shù)Table 3 JCO accident simulation parameter

JCO事故的裂變功率震蕩持續(xù)時(shí)間很長(zhǎng)，由于事故發(fā)生后工廠采取了多種干預(yù)措施，對(duì)事故全過程的模擬非常困難。臨界事故發(fā)生后人員受到的輻照程度主要由臨界事故的第一裂變功率峰產(chǎn)額決定，因此本文主要對(duì)第一裂變功率峰裂變次數(shù)進(jìn)行了模擬驗(yàn)證，其計(jì)算值和估計(jì)值分別為4.54×1016和5×1016，二者符合較好。

圖4示出程序計(jì)算得到的功率和溫度。由圖4a可見，程序模擬結(jié)果能較好地體現(xiàn)臨界事故的發(fā)展過程，即功率在短時(shí)間內(nèi)上升到一較高峰值后迅速下降，最終趨向于一穩(wěn)定的較低值。在這一過程中，由于輻照裂解氣體逸出水面等因素的影響，在第一個(gè)功率峰值之后還可能會(huì)出現(xiàn)幾個(gè)較低的峰值。

圖4 JCO事故功率與溫度的變化Fig.4 Power and temperature variations of JCO accident

5 結(jié)論

本文針對(duì)傳統(tǒng)的基于點(diǎn)堆模型的溶液瞬態(tài)分析程序難以適應(yīng)復(fù)雜幾何及復(fù)雜物料情況的事實(shí)，基于蒙特卡羅均勻化理論與有限體積方法建立了適用于瞬態(tài)分析問題的三維擴(kuò)散時(shí)空動(dòng)力學(xué)模型，然后將該模型與非穩(wěn)態(tài)傳熱模型、輻照裂解氣泡模型相耦合，對(duì)GETAC-S程序進(jìn)行升級(jí)。改進(jìn)后的GETAC-S從截面生成端到時(shí)空動(dòng)力學(xué)的求解，均實(shí)現(xiàn)了對(duì)復(fù)雜幾何、材料及能譜條件下模型的適配，從而使程序可針對(duì)任意條件下的溶液系統(tǒng)進(jìn)行瞬態(tài)特性分析。采用日本TRACY瞬態(tài)實(shí)驗(yàn)裝置的臨界事故模擬實(shí)驗(yàn)數(shù)據(jù)對(duì)改進(jìn)后的GETAC-S的可靠性進(jìn)行了驗(yàn)證，結(jié)果在10%內(nèi)符合。針對(duì)日本JCO事故進(jìn)行了事故進(jìn)程反演，取得了較好結(jié)果。驗(yàn)證結(jié)果表明改進(jìn)后的GETAC-S具備對(duì)復(fù)雜條件下溶液系統(tǒng)的臨界事故后果進(jìn)行評(píng)價(jià)與反演的能力，可為今后溶液核系統(tǒng)的設(shè)計(jì)起到一定的指導(dǎo)作用。