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

(1.中國原子能科學研究院 反應堆工程技術研究所，北京 102413；2.國家電投集團科學技術研究院有限公司，北京 102209)

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

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

1 基礎理論

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

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

1.1 蒙特卡羅均勻化理論

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

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

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

(1)

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

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

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

(2)

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

(3)

1.2 多群時空中子動力學擴散方法

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

多群時空中子動力學擴散方程如式(4)所示，本文使用有限體積方法對其進行離散求解[6]。

g=1,…,G

(4)

i=1,…,ND

(5)

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

2 軟件實現(xiàn)

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

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

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

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

2.2 基于OpenFOAM的三維擴散時空動力學程序

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

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

圖1 基于OpenFOAM求解時空動力學方程的基本流程Fig.1 Basic process of solving spatiotemporal dynamics equation based on OpenFOAM

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

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

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

(6)

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

(7)

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

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

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

3 TRACY實驗驗證

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

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

表1 RUN72實驗模擬參數(shù)Table 1 RUN72 experiment simulation parameter

表2 RUN72實驗計算結(jié)果及相對偏差Table 2 RUN72 experiment calculation result and relative deviation

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

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

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

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

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

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

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

5 結(jié)論

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