基于蒙特卡羅均勻化理論與有限體積方法的溶液系統臨界事故分析方法

孫 旭，周 琦,*，余慧鶯，朱慶福，夏兆東，寧 通，馬驍笛

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

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

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

1 基礎理論

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

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

1.1 蒙特卡羅均勻化理論

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

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

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

(1)

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

內積函數的統一性，如空間均勻化與能量積分通量可表示為：

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

(2)

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

(3)

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

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

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

g=1,…,G

(4)

i=1,…,ND

(5)

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

2 軟件實現

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

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

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

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

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

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

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

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

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

2.3 熱工反饋與耦合的實現

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

(6)

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

(7)

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

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

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

3 TRACY實驗驗證

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

圖2 TRACY裝置結構[12]Fig.2 TRACY device structure[12]

表1 RUN72實驗模擬參數Table 1 RUN72 experiment simulation parameter

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

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

4 瞬態分析程序在事故反演上的應用

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

表3 JCO事故模擬參數Table 3 JCO accident simulation parameter

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

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

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

5 結論

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