張彬航，張 聰，畢彥釗，張永紅，袁顯寶，唐海波，馮虹瑛

(三峽大學(xué) 機械與動力學(xué)院，湖北 宜昌 443002)

反應(yīng)堆燃耗計算是反應(yīng)堆設(shè)計與分析的重要環(huán)節(jié)[1]。在實際計算中，其核心是求解燃耗方程，并通常假設(shè)反應(yīng)率在單個時間步長內(nèi)恒定，因此燃耗方程可被視為一階線性方程組。其矩陣方程和指數(shù)解的形式如下：

(1)

(2)

式中：N(t)∈Rn為t時刻的核素密度向量；A∈Rn×n為燃耗系數(shù)矩陣，由系統(tǒng)中所有核素的反應(yīng)率和衰變率組成。

由于堆內(nèi)核素繁多，且不同核素的衰變常量數(shù)量級跨度極大，使得燃耗方程組具有大型、稀疏、剛性的特性，常規(guī)解法無法精確計算此問題。為精確計算矩陣指數(shù)eAt，近年來諸多基于近似理論和矩陣?yán)碚摰挠嬎惴椒ǖ玫搅藦V泛的發(fā)展與應(yīng)用，如切比雪夫有理近似方法(CRAM)和圍道積分有理近似方法(QRAM)[2-3]。CRAM基于復(fù)平面將指數(shù)函數(shù)的最優(yōu)有理近似式擴展到矩陣指數(shù)上以逼近eAt：

(3)

QRAM則是基于柯西積分公式和圍道求積組將式(2)轉(zhuǎn)換為復(fù)平面上的積分表達式：

(4)

式中：C為包圍矩陣At所有本征值的閉合圍道；zk為求積組的插值點；ck為相應(yīng)插值點的權(quán)重；N為展開階數(shù)。點燃耗程序MODEC采用32階QRAM和優(yōu)化后的拋物線圍道求積組用于燃耗計算中，取得了較好的計算精度和效率[5]。

由式(3)、(4)可發(fā)現(xiàn)，為求解t時刻的核素密度，CRAM和QRAM需分別進行k/2和N次大規(guī)模稀疏矩陣的求逆計算，且所有計算皆為復(fù)數(shù)運算，計算效率有限。

本文研究一種基于最佳一致逼近多項式(MMPA)的燃耗計算方法，只需進行1次矩陣求逆計算即可求解燃耗方程，具有較高的計算效率。并研制點燃耗程序AMAC，與蒙特卡羅輸運程序OpenMC進行耦合，通過計算衰變和固定輻照例題、OECD/NEA壓水堆柵元燃耗基準(zhǔn)題和沸水堆組件燃耗基準(zhǔn)題初步驗證該方法的正確性和有效性。

1 最佳一致逼近多項式方法

基于近似理論和矩陣?yán)碚揫6]，矩陣指數(shù)eAt可由某一確定矩陣函數(shù)f(At)進行逼近。若設(shè)矩陣At的特征值為{λ1,λ2，…，λn}∈C,則矩陣At可通過非奇異矩陣P轉(zhuǎn)換為Jordan標(biāo)準(zhǔn)型：

(5)

式中：J(λi)為對應(yīng)矩陣特征值λi的Jordan塊；矩陣P和P-1對應(yīng)于每個Jordan塊J(λi)，可轉(zhuǎn)換為：

(6)

對于矩陣At的特征值{λ1，λ2，…，λn}，若可為每個特征值λi定義函數(shù)f(λi)及其到ki-1階的導(dǎo)數(shù)fki-1(λi)，則矩陣函數(shù)f(At)的一般表達式[7]為：

(7)

式中：ki為最大Jordan塊J(λi)的大?。籲為矩陣At的特征值個數(shù)；Gi為(At-λiI)j在其廣義本征空間上的投射矩陣。由式(7)可得被逼近矩陣指數(shù)eAt的一般表達式[7]為：

(8)

由式(7)、(8)可知，對于矩陣At的每個特征值λi，若函數(shù)f(λi)及其導(dǎo)數(shù)f(j)(λi)在逼近指數(shù)eλi時有足夠小的誤差，則矩陣指數(shù)eAt可由矩陣函數(shù)f(At)進行最佳近似。根據(jù)文獻[8]可知，燃耗矩陣At的本征值{λ1，λ2，…，λn}位于復(fù)平面上的負實軸附近。因此，在矩陣函數(shù)f(At)對矩陣指數(shù)eAt進行最佳逼近時，應(yīng)滿足如下兩個條件：1) |f(At)-eAt|在負實軸附近足夠?。?) |f(j)(λi)-eλi|/j!在負實軸附近足夠小。為求出滿足上述條件的函數(shù)f(At)，本文基于MMPA對矩陣指數(shù)eAt進行了最佳近似，定義多項式函數(shù)f(t)為：

(9)

式中：a0、ai∈R為多項式系數(shù)；n為展開階數(shù)。

(10)

則通過變量代換后可得到指數(shù)函數(shù)ex的近似表達式：

(11)

式中，x∈(-∞,0]。根據(jù)文獻[9]可知，當(dāng)常數(shù)c取24.1并進行32階展開時，最大近似誤差的最小值為2.2×10-14，能滿足燃耗計算的高精度要求。基于Remez算法[10]可得32階MMPA方法的實數(shù)系數(shù)ai。由于系數(shù)ai皆為實數(shù)，因此在進行燃耗計算時不需進行復(fù)矩陣運算，這是MMPA方法相對于CRAM、QRAM等有理近似方法的優(yōu)點之一。由式(2)和(11)可得t時刻的核素密度N(t)為：

(At-cI)-1}iN(0)

(12)

式(12)中，矩陣(At+cI)(At-cI)-1的特征值λ′i為：

(13)

由于矩陣At的特征值{λ1,λ2，…，λn}∈(-∞,0]，因此特征值λ′i∈[-1,1)，則在進行多項式逼近過程中{(At+cI)(At-cI)-1}i的特征值始終在區(qū)間[-1，1]上，具有良好的數(shù)值穩(wěn)定性。如式(12)所示，為了求得t時刻的核素密度，只需獲得矩陣(At+cI)(At-cI)-1后再進行31次迭代相乘即可。相對于目前廣泛應(yīng)用于燃耗計算的16階CRAM和32階QRAM分別所需的8次和32次復(fù)矩陣求逆相比，MMPA方法只需進行1次矩陣的求逆運算，在計算效率與方法實現(xiàn)上具有一定的優(yōu)勢。因此，本文進一步基于32階MMPA方法開發(fā)了點燃耗程序AMAC，并耦合蒙特卡羅輸運程序OpenMC以驗證該方法的正確性和有效性。

2 程序開發(fā)

AMAC采用了精細燃耗數(shù)據(jù)庫進行燃耗計算。該燃耗數(shù)據(jù)庫整合了ORIGEN-2與ORIGEN-S的最新燃耗數(shù)據(jù)庫，共包含1 487種不重復(fù)核素、11種衰變反應(yīng)、23種中子截面反應(yīng)和30種重核裂變產(chǎn)額數(shù)據(jù)[11-12]。同時AMAC也能直接讀取ORIGEN-2或ORIGEN-S的燃耗數(shù)據(jù)庫進行燃耗計算。在基于MMPA方法求解燃耗方程的過程中，由式(12)可知只需對矩陣(At-cI)進行1次求逆運算，考慮到該矩陣具有大型、稀疏、剛性的特性，直接采用高斯消元法會引入較大的舍入誤差。為保證求解的精度和效率，AMAC首先對該矩陣進行符號LU分解，獲得相應(yīng)的填入矩陣，然后對填入矩陣進行高斯消元，從而完成矩陣(At-cI)的求逆[13]。另外考慮到矩陣的稀疏性，為節(jié)省內(nèi)存及計算時間，AMAC采用三元組方法對稀疏矩陣進行壓縮儲存。

AMAC的簡要程序架構(gòu)如圖1所示。程序主要分為前處理、求解器及后處理3部分。前處理部分主要讀取燃耗計算中所需的數(shù)據(jù)，包括衰變常量、微觀截面、裂變產(chǎn)額、通量、時間步長等信息，并將這些數(shù)據(jù)處理成燃耗矩陣的形式并存儲，由求解器直接調(diào)用計算。求解器部分主要是根據(jù)不同的計算方法進行燃耗方程的求解。除32階MMPA方法外，AMAC還實現(xiàn)了16階CRAM和32階QRAM的燃耗計算方法。后處理部分將計算結(jié)果進行格式化輸出，輸出數(shù)據(jù)包括核素密度、放射性活度、衰變熱等。在完成點燃耗程序AMAC研制的基礎(chǔ)上，進一步基于Python完成了與OpenMC的耦合以便計算真實的燃耗問題，并引入預(yù)估-校正和子步法兩種耦合策略以保證燃耗計算的精度與效率。

圖1 AMAC程序框架簡圖Fig.1 Brief structure of AMAC code

3 數(shù)值驗證與分析

對于MMPA方法的正確性與有效性的驗證，主要基于AMAC來完成，主要分為兩個方面：一是在相同數(shù)據(jù)庫下，計算了237Np的衰變和UO2燃料的固定輻照問題，并與ORIGEN-2和RMC程序進行對比，以驗證MMPA方法在求解燃耗方程時的精度與效率；二是通過耦合程序OpenMC-AMAC完成OECD/NEA壓水堆柵元燃耗基準(zhǔn)題和沸水堆組件燃耗基準(zhǔn)題的計算與驗證分析，以評估MMPA方法計算實際燃耗問題的能力。

3.1 237Np的衰變計算

衰變例題計算了初始含量為1 mol的237Np衰變106a后各核素成分的變化。參考程序選取了ORIGEN-2和RMC，皆基于ORIGEN-2的燃耗庫完成計算，計算結(jié)果和相對偏差列于表1。以RMC的線性子鏈解析方法的計算結(jié)果作為基準(zhǔn)解，分別列出了ORIGEN-2和AMAC相對于基準(zhǔn)解的相對偏差，可看到MMPA方法的計算結(jié)果沿燃耗鏈長度方向的相對偏差沒有增大趨勢，且與基準(zhǔn)解吻合良好，最大相對偏差僅為0.03%。對于ORIGEN-2的計算結(jié)果，由于時間步長為1.0×106a，根據(jù)ORIGEN-2的處理機制，當(dāng)核素半衰期小于燃耗時間步長的1/10時，就判定該核素為短壽命核素，需從燃耗矩陣中移除后進行平衡性假設(shè)處理，而燃耗矩陣中僅保留長壽命核素，以解決燃耗矩陣的剛性問題。在該例題中，沿燃耗鏈上的核素233Pa、225Ra、225Ac、221Fr、217At、213Bi、213Po、209TI、209Pb被判定為短壽命核素且需進行平衡性假設(shè)處理，導(dǎo)致這些核素“過度”衰變，使得這些核素的計算結(jié)果小于解析方法的計算結(jié)果，并最終使得燃耗鏈末端核素209Bi的計算結(jié)果大于解析解。MMPA方法計算時直接對燃耗矩陣At進行處理，不需單獨處理短壽命核素，因此燃耗鏈長度的增長不會導(dǎo)致計算精度的明顯下降。同時基于Remez算法求解多項式系數(shù)時嚴(yán)格遵守了對矩陣指數(shù)eAt進行最佳近似的兩個約束條件，能保證求解精度。因此，MMPA方法對于長燃耗鏈的計算具有良好的精度和穩(wěn)定性。在計算效率方面，ORIGEN-2耗時0.002 s，而MMPA方法的耗時僅為0.001 s。

3.2 UO2燃料的固定輻照計算

固定輻照例題計算了初始含量為1 mol、富集度為5%的UO2燃料，在3.0×1014cm-2·s-1的中子注量率下輻照1 a再冷卻3 a后的各核素成分變化情況。輻照時間和冷卻時間各劃分10個子步長，參考程序選取RMC，使用其點燃耗計算功能完成了該例題的計算[14]。表2列出重要錒系核素和部分代表性核素的計算結(jié)果。與RMC所采用的16階CRAM的計算結(jié)果相比，AMAC所采用MMPA方法的計算結(jié)果與其最大相對偏差不超過0.02%，具有較高的數(shù)值精度和穩(wěn)定性。在計算效率方面，AMAC還使用了16階CRAM與32階QRAM對該例題進行了測試，在保留4位有效數(shù)字的情況下，這兩種計算方法的結(jié)果與32階MMPA方法完全一致。各方法的計算時間列于表3，可看到MMPA方法相對于其他方法在計算效率上有一定的優(yōu)勢。

表1 237Np衰變鏈的計算結(jié)果Table 1 Calculation result of 237Np decay chain

表2 5% UO2燃料的計算結(jié)果Table 2 Calculation result of 5% UO2 fuel

表3 不同計算方法的計算時間Table 3 Computing time of different methods

3.3 OECD/NEA壓水堆柵元燃耗基準(zhǔn)題

該基準(zhǔn)題的計算目標(biāo)是通過計算壓水堆組件中的單個柵元以比較不同程序?qū)θ剂现型凰爻煞值挠嬎闩c分析能力。根據(jù)功率和最終燃耗深度的不同，該基準(zhǔn)題包括工況A(27.35 GW·d/tHM)、工況B(37.12 GW·d/tHM)和工況C(44.34 GW·d/tHM)，詳細幾何與材料參數(shù)參考文獻[15]。

耦合程序OpenMC-AMAC采用MMPA方法和ORIGEN-2燃耗數(shù)據(jù)庫完成了工況C的計算。同時基于數(shù)據(jù)庫ENDF/B-Ⅶ，Serpent程序采用16階CRAM也完成了計算[16]，計算結(jié)果列于表4。由表4可見，OpenMC-AMAC計算得到的大部分核素質(zhì)量結(jié)果與實驗值的相對偏差在5%以內(nèi)，表明MMPA方法具有良好的數(shù)值精度和穩(wěn)定性。其中對于裂變產(chǎn)額的計算，除核素149Sm外，其余核素的計算結(jié)果均在合理范圍之內(nèi)。149Sm的相對偏差為-26.493 6%，可能是149Sm的產(chǎn)生量太少或消失量太多所導(dǎo)致。149Sm從其先驅(qū)核149Nd經(jīng)過兩次β-衰變而來，而其消失途徑主要為(n，γ)反應(yīng)，所以O(shè)RIGEN-2的衰變數(shù)據(jù)和裂變產(chǎn)額會影響149Sm的產(chǎn)生量，OpenMC給出的(n，γ)截面會影響149Sm的消失量。因此，OpenMC計算的截面數(shù)據(jù)或ORIGEN-2的衰變數(shù)據(jù)的準(zhǔn)確度會導(dǎo)致149Sm產(chǎn)生較大偏差。進一步發(fā)現(xiàn)Serpent的計算結(jié)果與AMAC吻合較好，但同樣與測量值相差很大，因此該基準(zhǔn)題的149Sm測量值可能存在問題[17]。對于錒系核素的計算，238Pu的相對偏差相對于其他一些錒系核素較大，為-4.828 9%。在燃耗計算中，238Pu、237Np這些核素主要由其他核素衰變產(chǎn)生，衰變庫的數(shù)據(jù)準(zhǔn)確度直接決定了這些核素的計算精度。同時OpenMC的統(tǒng)計誤差也會影響核素的最終計算精度。為精確比較不同燃耗計算方法的耗時情況，統(tǒng)計了OpenMC-AMAC與Serpnet僅在進行燃耗計算時的耗時情況，其中AMAC耗時26.45 s，Serpent耗時26.41 s，兩者效率相當(dāng)。

3.4 OECD/NEA沸水堆組件燃耗基準(zhǔn)題

該基準(zhǔn)題由OECD/NEA發(fā)布，旨在比較不同燃耗程序在計算含可燃毒物組件時的能力。該沸水堆組件采用了8×8的排列方式進行燃料棒布置，其中包含8根含釓燃料棒以吸收壽期初的剩余反應(yīng)性?？偟娜己纳疃葹?0 MW·d/kgHM，功率密度為25.6 W/gHM，詳細的幾何與材料參數(shù)參考文獻[18]。

表5列出核素密度的16套參考程序的計算結(jié)果、平均參考值和OpenMC-AMAC的計算結(jié)果。由表5可看出，本文計算結(jié)果均在參考值范圍內(nèi)，且大部分核素的計算結(jié)果與平均參考值的相對偏差在5%以內(nèi)。對于相對偏差較大的核素，如243Am及Gd同位素，由于缺乏實驗值和各程序所采用的數(shù)據(jù)庫、輸運方法及燃耗計算方法的不同，難以評價其計算結(jié)果的精確性。尤其對于Gd同位素，由于其中子吸收截面很大，對于燃耗步長、可燃毒物棒的徑向劃分等因素十分敏感。為進一步驗證MMPA方法的正確性和穩(wěn)定性，在燃耗步長相同的情況下，對8根含釓燃料棒進行了徑向分區(qū)和計算，表6列出徑向均分為3區(qū)、5區(qū)、7區(qū)、10區(qū)4種情況下Gd同位素的相對偏差。由表6可見，隨徑向分區(qū)數(shù)目的增加，Gd同位素的相對偏差顯著下降。由于Gd同位素的吸收截面很大，這一特性將引起燃耗過程中中子注量率分布在短時間內(nèi)的顯著變化，因此進行合理的徑向分區(qū)能有效提高釓棒的計算精度，得到合理的燃耗計算結(jié)果。

表4 工況C下OpenMC-AMAC核素質(zhì)量的計算結(jié)果Table 4 Calculation result of OpenMC-AMAC under condition C

表5 OpenMC-AMAC核素密度的計算結(jié)果Table 5 Calculation result of nuclear density of OpenMC-AMAC

續(xù)表5

表6 不同徑向分區(qū)下OpenMC-AMAC計算核素Gd的相對偏差Table 6 Relative deviation of Gd calculated by OpenMC-AMAC in different radial zones

表7列出組件的無限增殖因數(shù)k∞的計算結(jié)果，包括參考值、平均參考值、Serpent和OpenMC-AMAC隨燃耗深度變化的計算結(jié)果。由表7可見，MMPA方法得到的k∞在參考值的區(qū)間內(nèi)，且隨燃耗深度變化的趨勢相同，并與平均值的相對偏差較小。在燃耗計算效率方面，OpenMC-AMAC的燃耗計算耗時為1 664 s，Serpent耗時為1 710 s?？傮w來看，MMPA方法在進行燃耗計算時具有良好的計算精度與效率。

3.5 計算效率分析

通過上述基準(zhǔn)題的驗證與分析，可看到MMPA方法在燃耗計算中具有良好的計算精度和數(shù)值穩(wěn)定性。在計算效率方面，對于32階MMPA方法，由式(12)可知，為了求得t時刻的核素密度，只需獲得矩陣(At+cI)(At-cI)-1后再進行31次矩陣的迭代相乘和相加。對于16階CRAM方法，由式(3)可知，需對復(fù)矩陣進行8次求逆和7次相加即可求得t時刻的核素密度。理論上大型稀疏矩陣的1次求逆計算耗時高于矩陣的1次相乘或相加計算，因此MMPA方法的耗時應(yīng)比CRAM方法少，這一點通過上述例題的耗時分析得到了驗證，只是時間上的縮短效果有限。

表7 OpenMC-AMAC的k∞計算結(jié)果Table 7 k∞ calculation result of OpenMC-AMAC

為了分析實際計算效率與理論分析的差異，本文進一步統(tǒng)計了AMAC中涉及矩陣求逆、乘法、加法等計算的耗時情況，發(fā)現(xiàn)32階MMPA方法中計算矩陣(At-cI)-1的耗時與16階CRAM方法中進行矩陣αi(At+θiI)-1的計算耗時相當(dāng)，但在矩陣{(At+cI)(At-cI)-1}i的迭代相乘計算中，耗時較為嚴(yán)重，約6次矩陣(At+cI)(At-cI)-1的自迭代相乘與矩陣αi(At+θiI)-1的計算耗時相當(dāng)，從而導(dǎo)致32階MMPA方法在實際計算效率上沒有顯著優(yōu)于16階CRAM方法。本文嘗試使用C++矩陣處理庫Eigen替換AMAC中自編寫的相關(guān)矩陣計算函數(shù)后，重新計算了OECD/NEA沸水堆組件燃耗基準(zhǔn)題，發(fā)現(xiàn)32階MMPA方法相比于16階CRAM方法能加速約30%，加速效果較為顯著，但應(yīng)用于實際燃耗計算中還需大量的測試和程序優(yōu)化，后續(xù)工作中對于AMAC中矩陣的計算效率需重點研究與優(yōu)化。

4 結(jié)論

本文研究了一種基于最佳一致逼近多項式(MMPA)的燃耗計算方法。相比于CRAM和QRAM在求解中的多次復(fù)矩陣求逆，該方法僅需1次矩陣求逆即可求解燃耗方程，且所有計算皆為實數(shù)運算，數(shù)值穩(wěn)定性好且求解效率高。同時基于MMPA方法開發(fā)了點燃耗程序AMAC，并與OpenMC進行耦合。為驗證MMPA方法的正確性和有效性，基于AMAC對衰變問題和固定輻照問題進行了點燃耗測試，并與ORIGEN-2和RMC的計算結(jié)果進行了比較，結(jié)果吻合良好且計算效率較高?；贠penMC-AMAC對OECD/NEA壓水堆柵元燃耗基準(zhǔn)題和沸水堆組件燃耗基準(zhǔn)題分別進行了計算分析，計算結(jié)果與實驗值及各參考值吻合良好，且與Serpent所采用的16階CRAM的燃耗計算效率和精度相當(dāng)。上述例題及基準(zhǔn)例題的計算驗證了MMPA方法在理論和數(shù)值上的正確性和有效性。后續(xù)工作將重點優(yōu)化程序中的矩陣計算效率，進一步提高MMPA方法的計算效率。