反應(yīng)堆燃耗方程數(shù)值解法研究摘""要：在反應(yīng)堆燃耗計算中，每種參與燃耗反應(yīng)的核素都具有與時間相關(guān)的燃耗方程，從而可得到整體的燃耗矩陣方程。由于燃耗矩陣方程通常為大型稀疏剛性的病態(tài)指數(shù)矩陣，直接求解難度較高，因此一般采用數(shù)值處理方法化指數(shù)矩陣計算為矩陣乘法計算等，在保持計算精度的基礎(chǔ)上提高計算效率。主要探討了泰勒級數(shù)展開、切比雪夫有理近似、Krylov子空間法、龍格庫塔方法以及改進歐拉方法的計算能力?；贛ATLAB程序，開展了5種方法對典型燃耗鏈的計算驗證對比，結(jié)果表明CRAM在計算精度、計算效率和計算穩(wěn)定性上均有較為良好的表現(xiàn)，適用于大型精細燃耗計算。關(guān)鍵詞：燃耗計算"""數(shù)值求解"""切比雪夫有理近似"""Krylov子空間法中圖分類號：TL327Study"of"Numerical"Methods"of"Reactor"Burnup"EquationHU"Yuying1""LIAO"Hongkuan2，3*""YU"Yingrui3""ZHOU"Bingyan1""LI"Tianya2，31.National"Key"Laboratory"of"Reactor"Technology，"Nuclear"Power"Institute"of"China;"2.Chengdu"Hezong"Nuclear"Power"Engineering"Co.，"Ltd.;"3.Nuclear"Power"Institute"of"China，"Chengdu，"Sichuan"Province，"610041"ChinaAbstract："In"the"burnup"calculation，"each"nuclide"in"the"burnup"reactor"has"a"burnup"equation"with"the"time，"so"that"the"whole"burnup"equation"could"be"gained."Due"to"the"difficulty"in"directly"solving"the"exponential"matrix，nbsp;numerical"methods"which"translate"the"exponential"matrix"calculationnbsp;into"the"matrix"multiplication"are"generally"used"to"improve"computational"efficiency."In"this"study，"Taylor"series"expansion，"Chebyshev"Rational"Approximation"Method（CRAM），"Krylov"Subspace"Method，"Runge-Kutta"Method"and"Improved"Euler"Method"are"discussed."Through"the"five"methods"with"MATLAB"code"and"verification"of"some"typical"burnup"chains，"it"is"shown"that"CRAM"performs"the"best"in"accuracy，"efficiency"and"stability，"which"is"suitable"in"the"large"fine"depletion"calculation.Key"Words："Burnup"calculation;"Numerical"solution;"Chebyshev"Rational"Approximation"Method;"Krylov"sub-space"method燃耗計算是反應(yīng)堆物理計算的重要組成部分，它直接關(guān)系到整個反應(yīng)堆計算的計算精度。因此為保證反應(yīng)堆在整個計算過程中各核素構(gòu)成的準(zhǔn)確性，燃耗方程求解方法的研究十分重要。對反應(yīng)堆燃耗系統(tǒng)來說，每一個核素都可寫為與時間（t）相關(guān)的燃耗方程，其中燃耗方程求解的重點在于其中指數(shù)矩陣的求解。由于燃耗矩陣A通常為大型稀疏剛性的病態(tài)矩陣，因此通常采用數(shù)值近似對指數(shù)矩陣進行降階處理或化指數(shù)矩陣計算為矩陣乘法計算等，最終得到高精度的計算結(jié)果，同時減少計算時間，提高計算效率。本文通過分析各類數(shù)值求解方法，針對燃耗矩陣的自身特點，確定適用于燃耗矩陣求解的幾類主要數(shù)值分析方法，采用基準(zhǔn)例題等對各種方法的計算精度、效率和穩(wěn)定性進行分析。通過燃料鏈解析解驗證方法的適用性。1"數(shù)值計算方法指數(shù)矩陣數(shù)值求解方法很多，主要分為級數(shù)法、常微分方程求解法、多項式方法以及矩陣分解法等，下面介紹幾種常見的數(shù)值求解方法。泰勒級數(shù)展開[1]是處理指數(shù)函數(shù)的一種常見辦法，根據(jù)展開公式可得到：。假設(shè)，在實際計算過程中，可以根據(jù)，得到N值，達到計算截斷的目的，就是的近似解，即在時，核子密度：該方法[2-5]是在復(fù)平面上逼近指數(shù)的一種數(shù)值近似方法，也屬于一種級數(shù)方法。通過該方法可以將寫為：為切比雪夫有理近似（Chebyshev"Rational"Approximation"Method，CRAM）的階數(shù)，該方法的收斂精度約為，計算時間隨值線性增加。其中的和為給定參數(shù)，根據(jù)階數(shù)不同對應(yīng)不同的值。該方法[1，6]主要是對指數(shù)矩陣進行近似多項式表達，滿足：。為m-1階多項式，從Krylov子空間選?。骸Ｊ紫龋?biāo)準(zhǔn)正交基，得到。令，則；從而可以得到最終公式。由于，因此。子空間方法的主要功效在于降低矩陣階數(shù)，通過近似后的新指數(shù)矩陣依舊可以通過泰勒展開求解，此時由于矩陣階數(shù)大幅度減小，泰勒級數(shù)展開的計算時間也會大為縮減。龍格庫塔方法[7]也是一種常微分方法，四階該方法可寫為以下形式，該方法具有四階精度。該方法屬于常微分方法，是一種較為簡單的方法，從程序設(shè)計的角度上來說也是最快的一種方法。但該方法沒有利用到特定問題的線性、常系數(shù)特性。結(jié)合燃耗求解的問題，求解方法可寫為：該方法[1]是MATLAB自帶的指數(shù)矩陣求解函數(shù)使用的。尺度平方方法主要是為了規(guī)避在泰勒級數(shù)展開或帕德近似中時間隨或矩陣特征值增大而增加的問題。該方法主要利用，尋找一個，使其滿足。接著再使用泰勒級數(shù)方法或帕德近似方法求解，最后再對結(jié)果進行次方計算即可。一般來說，先尋找一個整數(shù)，使其能夠滿足，那么。在時，核子密度通過尺度平方方法和帕德近似的結(jié)合使用，可以得到：除以上介紹的方法以外，目前比較常用的方法還有拉蓋爾多項式逼近、求積組有理近似、逆拉普拉斯變換和拉格朗日插值法等。前兩種方法與切比雪夫有理近似相類似；采用拉普拉斯變換方法需要先將整體燃耗反應(yīng)分解為各條單鏈再進行計算；而拉格朗日插值等多項式方法對求解矩陣的要求較高，針對實際燃耗矩陣并不適用。2"數(shù)值計算方法分析與評價2.1"概述根據(jù)上述各數(shù)值計算方法和燃耗矩陣求解的特點，本文選擇了泰勒展開、切比雪夫有理近似、Krylov子空間方法、四階龍格庫塔和改進歐拉法5種數(shù)值計算方法進行分析研究，從數(shù)值計算方法的計算精度、計算效率和計算穩(wěn)定度角度進行分析評價。主要內(nèi)容如下。（1）采用MATLAB程序編碼，實現(xiàn)上述五種數(shù)值求解方法，CRAM求解階數(shù)為14階。MATLAB中EXPM（）函數(shù)進行指數(shù)矩陣的求解時采用尺度平方加帕德近似的數(shù)值方法；如果求解的矩陣含有全套的特征值向量和特征值，則采用特征值法求解。（2）采用隨機生成的計算矩陣，分析評價各數(shù)值分析方法對于不同階數(shù)下隨機矩陣的適用性。（3）根據(jù)燃耗矩陣病態(tài)問題比較嚴(yán)重的特點，采用隨機生成的大型病態(tài)矩陣，分析評價各數(shù)值計算方法的適用性。2.2"隨機矩陣計算分析計算所需的矩陣A和初始值由MATLAB程序隨機生成，矩陣階數(shù)從10階不斷增加到2"200階，計算結(jié)果如下。圖1中橫坐標(biāo)表示矩陣階數(shù)，縱坐標(biāo)為各數(shù)值方法計算結(jié)果與MATLAB自帶指數(shù)函數(shù)求解結(jié)果最大偏差的對數(shù)值，其中實線表示各種數(shù)值方法的最大相對偏差，虛線表示平均相對偏差。由圖中可以看出，泰勒展開法的精度最高，效果最好，2"000階矩陣的計算精度約在10-16～10-15量級；Krylov子空間方法除了開頭有較大偏差外其余值基本和泰勒展開的結(jié)果相近。這主要是由于泰勒展開的階數(shù)很高，因此精度高。同樣的子空間方法在完成矩陣降階后新生成的指數(shù)矩陣求解也是通過泰勒展開完成的，只是在矩陣降階過程中由于矩陣近似產(chǎn)生了一些誤差，因此比泰勒展開的計算結(jié)果差一些，2"000階矩陣的計算精度約在10-15～10-14量級。由低階矩陣的結(jié)果看，四階龍格庫塔方法的計算精度較高，但隨著矩陣階數(shù)增加，計算精度逐漸下降，最終精度約在10-7～10-6量級。由于該方法也是由泰勒展開推導(dǎo)而來，但相比泰勒法展開階數(shù)較少，因此精度在矩陣階數(shù)較高時大幅下降。同樣問題出現(xiàn)在改進歐拉法中，作為一個二階精度的常微分方程解法，其計算精度更不如四階龍格庫塔法，2"000階矩陣的計算精度約在0.1%量級。CRAM的最大特點在于計算精度較高，總體在10-12～10-11量級，并且在矩陣階數(shù)變化過程中始終維持這一精度，隨階數(shù)增加的精度下降量很小。從CRAM的平均偏差和最大偏差變化趨勢可以看出，兩者之間的差距很小，這也進一步說明了該方法的計算穩(wěn)定性。對燃耗方程來說，很適合求解精細燃耗鏈下的燃耗指數(shù)矩陣。根據(jù)MATLAB程序驗證，泰勒展開方法在低階矩陣處理上計算效率和其他方法相差不大，但處理高階矩陣時速度遠慢于其他5種方法，千階以上計算時間會達到數(shù)十秒，2"000階隨機矩陣的指數(shù)計算求解時間達到了140"s以上。這主要是由于泰勒展開的計算過程中會進行大量的矩陣相乘計算，隨矩陣階數(shù)增加，其中元素進行的相乘計算也會極大增加。圖2中為除泰勒展開法以外的其他數(shù)值方法的計算時間結(jié)果，如圖所示，矩陣階數(shù)較小時，計算速度都很快，后期CRAM耗時較多，但2"000階矩陣單次計算結(jié)果也為2.491"s；改進歐拉法由于計算步驟少計算效率很高，Krylov子空間法優(yōu)勢在于精簡了矩陣階數(shù)，除去降階所耗時間外，新的低階矩陣泰勒展開計算效率也很高。結(jié)合圖和計算結(jié)果可以看出，在矩陣階數(shù)在千階級以下時，計算效率都還不錯，階數(shù)繼續(xù)增長所需時間則大幅度增加。2.3"隨機病態(tài)矩陣計算分析在實際燃耗矩陣中，病態(tài)問題比較嚴(yán)重，與MATLAB生成的普通隨機矩陣有一定的差別。而Hilbert矩陣是一種典型的病態(tài)矩陣，矩陣中數(shù)值的微小變動即會帶來結(jié)果的巨大變化，與燃耗矩陣有一定的相似性，因此，利用MATLAB生成10階到2"100階的隨機Hilbert矩陣進行數(shù)值驗證，結(jié)果如圖3所示。從圖中可以看出，引入病態(tài)矩陣計算后，其他方法的計算精度偏差與隨機矩陣的情況基本類似，但Krylov子空間方法的計算精度非常差，2"000階矩陣的計算誤差達到了60%以上。這主要是因為Krylov方法進行數(shù)值計算的第一步是求出原矩陣的標(biāo)準(zhǔn)正交基進行矩陣降階，而病態(tài)矩陣在求解正交基過程中會引起極大的誤差，因此最終結(jié)果的偏差也非常大。根據(jù)圖4，從計算效率上看，各數(shù)值方法計算時間隨矩陣階數(shù)變化的趨勢同隨機矩陣的情況基本相同，依舊是泰勒展開方法耗時最多，千階矩陣計算時間達到了200"s以上。其次就是CRAM，在剩下的幾種方法中計算效率偏低一些，但2"000階矩陣也僅需2.586"s。Krylov子空間法和改進歐拉法耗時最少，這主要是由實際計算的矩陣階數(shù)和計算步驟決定的。總體來說，通過病態(tài)矩陣的數(shù)值驗證，大多數(shù)數(shù)值求解方法的計算精度、計算效率及計算穩(wěn)定性與隨機矩陣驗證相比都沒有太大變化。只有Krylov子空間方法在求解標(biāo)準(zhǔn)正交基降低矩陣階數(shù)的過程中由于病態(tài)矩陣的引入導(dǎo)致了極大的誤差。因此，如果燃耗矩陣的病態(tài)問題嚴(yán)重時，Krylov子空間法的適用性會極大的降低。2.4"結(jié)果分析通過MATLAB程序針對3種類型的隨機矩陣采用6種不同數(shù)值計算方法進行驗證，計算精度的比較以MATLAB程序自帶的指數(shù)矩陣函數(shù)EXPM（）為基準(zhǔn)進行比較，最終可以得到以下結(jié)論。（1）從計算精度和計算穩(wěn)定性上看，CRAM和泰勒展開方法的效果最好。CRAM的計算精度穩(wěn)定在10-12左右，隨矩陣階數(shù)增加計算精度基本保持不變；泰勒展開方法的計算精度在10-16左右，隨矩陣階數(shù)增加計算精度略微下降。（2）從計算效率看，400階及以下的矩陣在進行計算時，各種數(shù)值方法所耗費的時間都很短；矩陣階數(shù)增加后，泰勒展開方法的計算耗時大量增加，遠大于其他幾種方法。其余5種方法中，CRAM耗時最多，改進歐拉法耗時最少，但2"000階矩陣計算中CRAM的耗時也僅為3"s左右。（3）結(jié)合各方面因素可以看出，CRAM最適合于各類型的指數(shù)矩陣計算；而Krylov子空間法在非病態(tài)矩陣的計算中效果比CRAM更好，因此在非病態(tài)問題中該方法也很適用。在燃耗計算中，可以進一步就這兩種方法進行實際燃耗鏈的驗證。3"燃耗鏈驗證3.1"概述如前文所述，針對5種數(shù)值求解方法的隨機矩陣驗證分析，可以看出CRAM和Krylov子空間方法從計算精度和效率上都表現(xiàn)良好，因此，本文采用CRAM和Krylov子空間方法編寫燃耗計算程序，對部分典型燃耗[8]計算進行進一步驗證，與解析解[5]結(jié)果相對比。3.2"衰變單鏈計算分析本文選取211Pb衰變單鏈分析算法的適用性，最終生成的207Pb為穩(wěn)定核素，其燃耗鏈如下：式（6）中，211Pb、211Bi和207Tl的初始濃度分別為1×1010、1×104和1×101，由于該燃耗鏈中核素的半衰期都較短，因此驗證時間步長為1×104"s。根據(jù)表1中結(jié)果可以看出，兩種方法的計算精度都較好，精細燃耗步下Krylov子空間法的計算精度更高一些。由表中可以發(fā)現(xiàn)，Krylov子空間法受時間步長的影響非常大，100步和1"000步計算的計算精度相差較大。"CRAM的結(jié)果是一步計算得到的，而一步計算時Krylov子空間方法同樣無法計算。相同燃耗時間下，燃耗步增加則子空間法的計算精度提高，并且在粗燃耗步下由于燃耗矩陣的很大，子空間法將無法計算。因此，雖然計算結(jié)果表明子空間法計算精度高，但其計算效率和計算穩(wěn)定性很差。與此同時，CRAM受燃耗步劃分的影響極小，計算穩(wěn)定性很好，計算時間很少。3.3"復(fù)雜錒系鏈計算分析針對復(fù)雜燃耗鏈的計算分析，本文選取錒系鏈的計算分析評價方法的適用性。改燃耗鏈從238U的衰變和輸運開始，含有復(fù)雜的衰變和中子反應(yīng)，最大特點在于燃耗鏈中含有從244Cm到240Pu的衰變閉環(huán)，如圖5所示。燃料鏈從238U開始，初始濃度為1×1012，燃耗時間為8.64×104"s，計算結(jié)果如表2所示，Krylov子空間法為劃