楊萬奎，劉耀光，馬紀敏，楊 鑫，，王冠博，佘 頂

（1.中國工程物理研究院 核物理與化學研究所，四川 綿陽 621900；2.清華大學 工程物理系，北京 100084）

反應堆燃耗計算關系到反應堆核燃料管理、乏燃料后處理、燃耗信任制等，是反應堆領域的一重要分支。實現(xiàn)該目的的基本方法為中子輸運計算和核素燃耗計算耦合。中子輸運計算部分，基于確定論方法的一、二維幾何處理能力的輸運計算方式已不能滿足反應堆燃耗計算的需求，隨著計算機硬件的發(fā)展，基于蒙特卡羅方法的三維輸運計算程序MCNP 得到了廣泛應用［1］；核素的燃耗計算部分，基于矩陣指數(shù)法的單群點燃耗計算程序ORIGEN2占統(tǒng)治地位［2］。

因此，國內(nèi)外均開展了基于MCNP 和ORIGEN 耦合系統(tǒng)的開發(fā)。一種為外耦合方式，將兩個程序的計算結果相互提取作為自身的輸入文件，如美國愛達荷國家實驗室的MOCUP［3］、阿 拉 莫 斯 國 家 實 驗 室 的Monteburns［4］、清 華 大 學 的MCBurn［5］等；另 一 種 直接在接口程序中加入核素轉換、衰變計算功能，取代ORIGEN2程序計算，系統(tǒng)僅由MCNP和接口程序組成，減少了讀取文件的時間，直接通過變量傳遞，如美國Sandia國家實驗室開發(fā)的BURNCAL程 序［6］、清 華 大 學 的RMC［7］。本文采用第2種耦合方法，直接在程序中加入核素轉換和衰變的計算功能，采用MPI并行方式（國內(nèi)做過相關研究對比［8］），研制并行化蒙特卡羅燃耗程序MCBMPI，并以壓水堆柵元燃耗基準題對程序進行驗證。

1 程序原理

燃耗方程的求解采用3種方法結合，即截斷泰勒展開的矩陣指數(shù)法、TTA 線性子鏈解析法和高斯-賽德爾迭代法。核素的分類參考ORIGEN2，分為錒系核素及其子體、裂變產(chǎn)物和活化產(chǎn)物3類，且初始數(shù)據(jù)庫（衰變數(shù)據(jù)庫、截面文件庫）均采用該點燃耗程序配套庫，后續(xù)燃耗步實時計算生成點截面庫。

截斷泰勒展開的矩陣指數(shù)法的適用對象為長壽命核素。指數(shù)函數(shù)的級數(shù)表示矩陣指數(shù)，并將級數(shù)展開至足夠多項以達到指定的精度要求。其中，最關鍵的級數(shù)項的計算可用遞歸算法方便實現(xiàn)。TTA 線性子鏈解析法的適用對象為無長壽命先驅核的短壽命核素。此類核素導致分離本征值，不適于矩陣指數(shù)法求解，且在1個時間步長內(nèi)，可認為此類核素已達平衡，可用解析解直接表示。高斯-賽德爾迭代法的適用對象為具有長壽命先驅核的短壽命核素，通過迭代法直接求解方程組。

該燃耗計算方法的另一特點是非空系數(shù)矩陣的存儲技術。該方程組的系數(shù)矩陣維數(shù)可達1 700 多，且為大型稀疏矩陣。理論上可通過存儲所有的方程組系數(shù)，求解方程組得到所有核素的濃度。但隨核素的增多，系數(shù)矩陣的存儲量將急劇增大（若以雙精度浮點數(shù)8字節(jié)存儲1個系數(shù)，則1 700×1 700的系數(shù)矩陣存儲容量為22.049 MB），對系數(shù)索引極不利。因此采用系數(shù)索引方法，只存儲非零元素。從數(shù)據(jù)庫文件中讀取各核素的半衰期、衰變分支比、截面、裂變產(chǎn)額等計算相關數(shù)據(jù)后，確定各核轉換類型及其產(chǎn)物，用1個數(shù)組存放轉換率及產(chǎn)物標識；通過查詢該數(shù)組以獲取每個子代核素Xi的父代核素Xj標識，將Xi的所有父代核素Xj的轉換率存于一維浮點數(shù)組中；將父代核素標識存放于一維整型數(shù)組中；用計數(shù)器確定索引位置，實現(xiàn)非零系數(shù)矩陣的存儲。

MCBMPI程序采用模塊化思想，其中蒙特卡羅輸運采用MCNP 程序的并行化版本MCNP5MPI，自編并行子程序MPIORI2實現(xiàn)核素截面替換、核素轉變和衰變計算、核素信息的提取等功能，并行策略為對多燃耗區(qū)采用區(qū)域分解的MPI消息傳遞并行。程序整體框架如圖1所示。

1.1 區(qū)域分解策略

區(qū)域分解通過自動讀取輸運計算中建模中的燃耗區(qū)幾何柵元，并依次編號實現(xiàn)，之后對編號進行計算核心分配。一種方式為采用常用等間距分配（圖2a），按計算核心1～N，接收燃耗區(qū)域1～TotalBurnupZone的計算任務。這種方式的區(qū)域分解方法實現(xiàn)起來很直接。但每個計算核心得到的結果在匯總時，材料序號的排列不是順序排列的。雖可再排序，但增加了額外的計算任務。因此，需考慮一種可直接得到材料號順序排列的區(qū)域分解方法：即讓每個計算核心所計算的燃耗區(qū)域序號是連續(xù)的（圖2b）。

圖1 MCBMPI程序流程圖Fig.1 Calculation flow chart of MCBMPI

圖2 區(qū)域分解方式Fig.2 Type of zone decomposition

按此區(qū)域分解方法，程序實現(xiàn)過程中為包含所有的燃耗區(qū)，需用到向上取整函數(shù)ceiling（）。至此，完成了燃耗區(qū)的區(qū)域分解，將燃耗計算任務分配給各計算核心，通過獲取相關截面信息后可并行調(diào)用燃耗方程求解器實現(xiàn)燃耗的并行計算。

1.2 燃耗并行存儲策略

通常的燃耗計算并行存儲方式為對各燃耗區(qū)分別建立對應的文件夾來存儲，但在燃耗區(qū)較多時，文件夾的數(shù)量過多，管理起來較凌亂。出于以上考慮，在保存完整信息的前提下，根據(jù)MPI相關的文件I／O 并行處理方法，采取一新的并行存儲策略，即根據(jù)計算核心的數(shù)量來建立相應的文件夾數(shù)量存儲單個計算核心計算生成的結果。

燃耗并行計算中所需存儲的結果主要包括燃耗計算的輸出文件及輸出文件中提取出的核素信息。這些核素信息可分為3類：下一次輸運計算所需的核素成分；下一次核素轉變和衰變計算的繼承核素；下一次燃耗計算所需的替換截面核素。

通過采取上述的燃耗計算并行存儲策略可對所有必需的結果實現(xiàn)并行化處理，在全部燃耗區(qū)都完成核素的轉變與衰變計算后，可提供輸運計算所需的全部參數(shù)。

1.3 單群截面的生成

通過MCNP5MPI輸運計算可得到當前燃耗深度下系統(tǒng)的有效增殖因數(shù)、各燃耗區(qū)的中子通量密度以及可通過中子通量密度響應得到的量，如反應率、功率等［9］。另一方面，大部分錒系核素和裂變產(chǎn)物的各反應截面與具體裝置對象的能譜相關，在進行1個新的燃耗計算前需對重要核素的截面進行修正或直接替換。一種有效的截面替換方法是通過輸運計算實時地產(chǎn)生。

所需替換的截面類型根據(jù)不同的核素類型有所不同，但均只需考慮以下反應截面［2］：錒系核素的（n，γ）、（n，2n）、（n，3n）、（n，f）；裂變產(chǎn)物的（n，γ）、（n，2n）、（n，α）、（n，p）。各反應類型在MCNP5輸運截面庫中的對應關系列于表1。

表1 ORIGEN2反應類型與MCNP5截面庫之間的對應關系Table 1 Relationship between ORIGEN2reaction type and MCNP5cross section library

本燃耗計算子程序中采用的核素截面為單群微觀截面，相對于多群微觀截面，更有利于MCNP5輸運計算實時生成相應的替換截面，并在輸運計算時的誤差估計方面也放寬了要求，即為達到規(guī)定的誤差要求所需的迭代次數(shù)減少了。單群微觀截面的生成需保證反應率守恒，通常采用能譜權重得到，如式（1）所示。

輸運計算中，只能分別統(tǒng)計得到式（1）中分子和分母。分母為全部能量區(qū)間的總通量，MCNP5中只需應用體通量計數(shù)即可得到。分子為微觀反應率，MCNP5中可通過Fq乘子卡結合體通量計數(shù)得到。最終通過對MCNP5結果文件中通量和反應率的提取，得到替換核素的單群微觀截面。

2 程序驗證

2.1 OECD／NEA 壓水堆柵元基準題

采用OECD／NEA 發(fā)布的壓水堆柵元計算基準題［10］，對該程序進行驗證。該基準題的目的在于對不同研究機構提供的針對該簡化柵元模型的核素成分的計算結果進行比較?；鶞暑}結果發(fā)布時，包含了世界范圍內(nèi)16家研究機構的21套計算結果，主要是對系統(tǒng)反應性影響較顯著的核素在壽期末的濃度進行對比，以驗證燃耗計算程序及數(shù)據(jù)庫。對比的核素包括12種 錒 系 核 素（234U、235U、236U、238U、237Np、238Pu、239Pu、240Pu、241Pu、242Pu、241Am、243Am）和16 種 裂 變 產(chǎn) 物（95Mo、99Tc、101Ru、103Rh、109Ag、133Cs、135Cs、143Nd、145Nd、147Sm、149Sm、150Sm、151Sm、152Sm、153Eu、155Gd）。

該燃料樣品的化學分析結果來自于CE對14×14 的組件設計過程。為滿足該基準題的目的，使用實際棒尺寸，但為使鈾水比與實際的二維排列相當，對燃料芯體的尺寸進行了修改。燃料樣品經(jīng)歷4 個完整的燃耗循環(huán)，各反應堆循環(huán)的輻照時間和衰變時間列于表2。該基準題包含3 種工況分別對應于參考燃料棒3個軸向位置處的燃料樣品，各樣品的總燃耗深度不同。

表2 柵元計算基準題的運行歷史參數(shù)Table 2 Operation history data for benchmark problem pin-cell calculation

2.2 驗證結果及分析

通過對3個工況的燃料樣品應用MCBMPI程序進行計算，得到壽期末3個工況下的樣品核素濃度，為評價MCBMPI程序的計算結果，分別給出了MCBMPI計算值與測量值和計算平均值兩者的相對誤差，將3個工況的相對誤差集中示于圖3。

根據(jù)各核素濃度的量級范圍與核素濃度的測量精度，基準題中所采用的相對誤差限為10%，即相對誤差在±10%以內(nèi)即認為理論與實驗符合得很好。由圖3a可知：239Pu、240Pu、150Sm、152Sm、153Eu 5 種 核 素 濃 度 的MCBMPI計算值與測量值的相對誤差在10%附近；只有149Sm 的相對誤差較大；其余核素的相對誤差均在誤差限以內(nèi)，符合較好。

圖3b 中，237Np、239Pu、240Pu、150Sm、152Sm 5種核素濃度的計算值與平均值的相對誤差在10%附近；243Am、149Sm、155Gd 3種 核 素 濃 度 的相對誤差較大，但其相對誤差與基準題中給出的計算平均值的標準偏差相當，且所有計算值均在其他計算值的波動范圍內(nèi)；其余核素符合得很好，尤其是U 元素的各同位素和95Mo、99Tc、101Ru、103Rh、109Ag 等 核 素，相 對 誤 差 在3%以內(nèi)。

圖3 濃度計算值與測量值（a）和計算平均值（b）的相對誤差Fig.3 Relative error between calculation and experiment（a）and averaged calculation（b）concentrations

另外，濃度的誤差較大的核素大部分微量，如149Sm 和155Gd，其濃度在給出的28種核素中最低，此處的濃度是1g UO2中目標核素的質(zhì)量，3 個工況中這兩種核素的濃度均在μg·g-1量級，243Am 的則為10-4μg·g-1量級，很小的偏差就能引起較大的相對誤差；另一方面，截面及裂變產(chǎn)額的因素也不可忽視，這也是基準 題 中 給 出 的243Am、149Sm、155Gd 3 種 核 素 的計算平均值存在較大標準偏差的原因［10］。導致截面差異的原因除基準庫不同外，燃耗步長是重要因素。因為這3 種核素的熱中子吸收截面較大，對中子通量和能譜的影響較大，反過來又影響自身核素的燃耗計算，故燃耗步長對其影響較大。程序對該基準題兼顧計算時間與計算精度，設置的燃耗步數(shù)為4，計算結果也在基準題給出的計算值范圍內(nèi)。另外，程序中計算參數(shù)的選取亦可導致計算結果的偏差，需根據(jù)試算結果進行優(yōu)化，如核素選取時的截斷份額（核素繼承截斷份額設置為-25 次方，核素吸收截斷份額設置為-6次方）、核素轉換計算時內(nèi)步長的確定（內(nèi)步長設置為10d）。

3 加速比與并行效率

以工況A 為例說明采用并行計算在計算時間方面的優(yōu)勢。計算硬件條件為Intel?CoreTM2 Quad CPU Q8300 ＠ 2.50 GHz，3.25G 內(nèi)存，計算耗時及加速情況列于表3。其中，加速比為多核計算耗時與單核計算耗時之比，并行效率為加速比與計算核數(shù)之比。

表3 MCBMPI的并行耗時及計算加速Table 3 Collapsed time and parallel speedup of program MCBMPI

由于該程序系統(tǒng)是由兩個部分組成的，整體加速效果是兩部分的綜合作用，因此需分別討論其加速情況。對MCNP5MPI而言，可從單次MCNP5MPI的計算耗時得到該程序對于該基準題的并行加速情況，其結果列于表4。

表4 MCNP5MPI的并行耗時及計算加速Table 4 Collapsed time and parallel speedup of program MCNP5MPI

從表3、4的加速比及并行效率可知，該并行燃耗計算程序系統(tǒng)的并行效率低于MCNP5MPI的并行效率。這是由于程序系統(tǒng)融入了輸運計算和燃耗計算兩部分，MCNP5MPI只包含輸運計算，而燃耗計算采用區(qū)域分解并行。燃耗區(qū)小于計算核數(shù)時，并行效果不能完全體現(xiàn)出來，但若應用到多燃耗區(qū)的燃耗計算則能完全體現(xiàn)并行加速，這將在后續(xù)的實際應用中得以體現(xiàn)。

4 小結

本文開發(fā)的并行蒙特卡羅燃耗計算程序MCBMPI與壓水堆柵元燃耗基準題的對比結果表明：各核素的計算值是可信的，且并行計算結果與串行結果完全一致，該程序加速了蒙特卡羅燃耗計算，可用于多燃耗區(qū)的燃耗計算。尤其是U 元素的各同位素的計算值與實驗值與基準計算值均符合很好，表明該程序在燃料管理與換料分析方面都具有應用前景。

［1］ 鄧力，李剛.粒子輸運蒙特卡羅模擬現(xiàn)狀概述［J］.計算物理，2010，27（6）：791-798.DENG Li，LI Gang.A summarization on Monte Carlo simulation［J］.Chinese Journal of Computational Physics，2010，27（6）：791-798（in Chinese）.

［2］ CROFF A G.A user’s manual for ORIGEN2 computer code，ORNL／TM-7175［R］.US：Oak Ridge National Laboratory，1980.

［3］ MOORE R L，SCHNITZLER B G，WEMPLE C A，et al.MOCUP：MCNP-ORIGEN2coupled utility program，INEL-95／0523US［R］.US：Idaho National Engineering Laboratory，1995.

［4］ TRELLUE H R.Development of Monteburns：A code that links MCNP and ORIGEN2in an automated fashion for burnup calculations［R］.US：Los Alamos National Laboratory，1999.

［5］ 余綱林，王侃，王煜宏.MCBurn：MCNP和ORIGEN耦合程序系統(tǒng)［J］.原子能科學技術，2003，37（3）：250-254.YU Ganglin，WANG Kan，WANG Yuhong.MCBurn：A coupling package of program MCNP and ORIGEN［J］.Atomic Energy Science and Technology，2003，37（3）：250-254（in Chinese）.

［6］ PARMA E J.BURNCAL：A nuclear reactor burnup code using MCNP tallies，SAND2002-3868［R］. US：Sandia National Laboratory，2002.

［7］ 佘頂，王侃，余綱林.堆用蒙卡程序燃耗計算功能開發(fā)［J］.核動力工程，2012，33（3）：1-5.SHE Ding，WANG Kan，YU Ganglin.Development of burnup calculation function in reactor Monte Carlo code RMC［J］.Nuclear Power Engineering，2012，33（3）：1-5（in Chinese）.

［8］ 鄧力，張文勇，徐涵，等.蒙特卡羅程序MCNP-Ⅱ與MCNP-5并行效率比較［J］.計算機工程與科學，2009，31（A1）：185-187.DENG Li，ZHANG Wenyong，XU Han，et al.The comparison of parallel efficiency between Monte Carlo code MCNP-Ⅱand MCNP-5［J］.Computer Engineering ＆ Science，2009，31（A1）：185-187（in Chinese）.

［9］ X-5 Monte Carlo Team. MCNP：A general Monte Carlo N-particle transport code，Version 5，LA-UR-03-1987［R］.US：Los Alamos National Laboratory，2003.

［10］de HART M D，BRADY M C，PARKS C V.OECD／NEA burnup credit calculational criticality benchmark phase Ⅰ-B results，ORNL-6901［R］.US：Oak Ridge National Laboratory，1996.