基于最佳一致逼近多項式的燃耗計算方法研究

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

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

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

(1)

(2)

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

由于堆內核素繁多，且不同核素的衰變常量數量級跨度極大，使得燃耗方程組具有大型、稀疏、剛性的特性，常規解法無法精確計算此問題。為精確計算矩陣指數eAt，近年來諸多基于近似理論和矩陣理論的計算方法得到了廣泛的發展與應用，如切比雪夫有理近似方法(CRAM)和圍道積分有理近似方法(QRAM)[2-3]。CRAM基于復平面將指數函數的最優有理近似式擴展到矩陣指數上以逼近eAt：

(3)

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

(4)

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

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

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

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

基于近似理論和矩陣理論[6]，矩陣指數eAt可由某一確定矩陣函數f(At)進行逼近。若設矩陣At的特征值為{λ1,λ2，…，λn}∈C,則矩陣At可通過非奇異矩陣P轉換為Jordan標準型：

(5)

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

(6)

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

(7)

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

(8)

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

(9)

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

(10)

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

(11)

式中，x∈(-∞,0]。根據文獻[9]可知，當常數c取24.1并進行32階展開時，最大近似誤差的最小值為2.2×10-14，能滿足燃耗計算的高精度要求。基于Remez算法[10]可得32階MMPA方法的實數系數ai。由于系數ai皆為實數，因此在進行燃耗計算時不需進行復矩陣運算，這是MMPA方法相對于CRAM、QRAM等有理近似方法的優點之一。由式(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的特征值始終在區間[-1，1]上，具有良好的數值穩定性。如式(12)所示，為了求得t時刻的核素密度，只需獲得矩陣(At+cI)(At-cI)-1后再進行31次迭代相乘即可。相對于目前廣泛應用于燃耗計算的16階CRAM和32階QRAM分別所需的8次和32次復矩陣求逆相比，MMPA方法只需進行1次矩陣的求逆運算，在計算效率與方法實現上具有一定的優勢。因此，本文進一步基于32階MMPA方法開發了點燃耗程序AMAC，并耦合蒙特卡羅輸運程序OpenMC以驗證該方法的正確性和有效性。

2 程序開發

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

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

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

3 數值驗證與分析

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

3.1 237Np的衰變計算

衰變例題計算了初始含量為1 mol的237Np衰變106a后各核素成分的變化。參考程序選取了ORIGEN-2和RMC，皆基于ORIGEN-2的燃耗庫完成計算，計算結果和相對偏差列于表1。以RMC的線性子鏈解析方法的計算結果作為基準解，分別列出了ORIGEN-2和AMAC相對于基準解的相對偏差，可看到MMPA方法的計算結果沿燃耗鏈長度方向的相對偏差沒有增大趨勢，且與基準解吻合良好，最大相對偏差僅為0.03%。對于ORIGEN-2的計算結果，由于時間步長為1.0×106a，根據ORIGEN-2的處理機制，當核素半衰期小于燃耗時間步長的1/10時，就判定該核素為短壽命核素，需從燃耗矩陣中移除后進行平衡性假設處理，而燃耗矩陣中僅保留長壽命核素，以解決燃耗矩陣的剛性問題。在該例題中，沿燃耗鏈上的核素233Pa、225Ra、225Ac、221Fr、217At、213Bi、213Po、209TI、209Pb被判定為短壽命核素且需進行平衡性假設處理，導致這些核素“過度”衰變，使得這些核素的計算結果小于解析方法的計算結果，并最終使得燃耗鏈末端核素209Bi的計算結果大于解析解。MMPA方法計算時直接對燃耗矩陣At進行處理，不需單獨處理短壽命核素，因此燃耗鏈長度的增長不會導致計算精度的明顯下降。同時基于Remez算法求解多項式系數時嚴格遵守了對矩陣指數eAt進行最佳近似的兩個約束條件，能保證求解精度。因此，MMPA方法對于長燃耗鏈的計算具有良好的精度和穩定性。在計算效率方面，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列出重要錒系核素和部分代表性核素的計算結果。與RMC所采用的16階CRAM的計算結果相比，AMAC所采用MMPA方法的計算結果與其最大相對偏差不超過0.02%，具有較高的數值精度和穩定性。在計算效率方面，AMAC還使用了16階CRAM與32階QRAM對該例題進行了測試，在保留4位有效數字的情況下，這兩種計算方法的結果與32階MMPA方法完全一致。各方法的計算時間列于表3，可看到MMPA方法相對于其他方法在計算效率上有一定的優勢。

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

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

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

3.3 OECD/NEA壓水堆柵元燃耗基準題

該基準題的計算目標是通過計算壓水堆組件中的單個柵元以比較不同程序對燃料中同位素成分的計算與分析能力。根據功率和最終燃耗深度的不同，該基準題包括工況A(27.35 GW·d/tHM)、工況B(37.12 GW·d/tHM)和工況C(44.34 GW·d/tHM)，詳細幾何與材料參數參考文獻[15]。

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

3.4 OECD/NEA沸水堆組件燃耗基準題

該基準題由OECD/NEA發布，旨在比較不同燃耗程序在計算含可燃毒物組件時的能力。該沸水堆組件采用了8×8的排列方式進行燃料棒布置，其中包含8根含釓燃料棒以吸收壽期初的剩余反應性。總的燃耗深度為40 MW·d/kgHM，功率密度為25.6 W/gHM，詳細的幾何與材料參數參考文獻[18]。

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

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

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

續表5

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

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

3.5 計算效率分析

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

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

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

4 結論

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