譜虧損修正法在點(diǎn)燃耗方程中的應(yīng)用

摘要：點(diǎn)燃耗方程是反應(yīng)堆物理最重要的微分方程之一，用于分析反應(yīng)堆內(nèi)核素隨時間的變化關(guān)系。在精細(xì)燃耗計(jì)算中，需要考慮1 000多種核素，且由于燃耗方程的剛性很強(qiáng)，使其數(shù)值求解難度大。譜虧損校正（Spectral Deferred Correction，SDC）法是一種新的用于復(fù)雜核素系統(tǒng)點(diǎn)燃耗計(jì)算的方法。該方法基于Gauss-Radau II-A 求積點(diǎn)，適用于剛性問題，且具備良好的計(jì)算效率和精度。最后，采用了3個算例檢驗(yàn)SDC的計(jì)算精度和效率，并和國際知名算法進(jìn)行對比。

關(guān)鍵詞：點(diǎn)燃耗 譜虧損修正法 微分方程 剛性

Application of Spectral Deferred Correction Method on the Point Burnup Equation

CAI Yun*" LI Qing" LI Tianya YANG Lingfang

Key Laboratory of Nuclear Reactor System Design Technology， Nuclear Power Institute of China， Chengdu， Sichuan Province， 610213 China

Abstract： The point burnup equation is one of the most important differential equations in reactor physics， which is used to analyze the relationship between reactor nuclides and time. In the fine burnup calculation， more than 1000 nuclides need to be considered， and the stiffness of the burnup equations makes it difficult to be solved. In this study， a new Spectral Deferred Correction （SDC）" method for point burnup equation in complex nuclide systems is proposed. This method is based on Gauss-Radau II-A to find the quadrature points， which is suitable for stiff problems and has good computational efficiency and accuracy. Finally， three examples are used to test the calculation accuracy and efficiency of SDC， compared with internationally renowned algorithms.

Key Words： Point burnup; Spectral deferred correction; Differential equation; Stiffness

點(diǎn)燃耗方程描述核反應(yīng)系統(tǒng)中核素隨時間的變化規(guī)律。 （1）為第i個核素的有效衰變常數(shù)，為第i個核素反應(yīng)生成第j個核素的有效份額。和分別由下式計(jì)算得到： （2）為第i個核素的衰變常數(shù)，為中子通量，為第i個核素反應(yīng)生成第j個核素的反應(yīng)截面。在實(shí)際燃耗計(jì)算過程中，通常假設(shè)每個燃耗步長內(nèi)的截面、通量等參數(shù)為常量，從而可將點(diǎn)燃耗方程當(dāng)作一階線性常微分方程處理。 （3），則式（3）的解可以寫成含矩陣指數(shù)的形式： （4）復(fù)雜核素系統(tǒng)覆蓋的同位素多達(dá)千余種，其中涉及大量短半衰期核素。例如：ORIGEN-S燃耗數(shù)據(jù)庫包含約1 500種核素，其中，某些核素（如212Po）的半衰期在1E-7 s量級。因此，復(fù)雜核素系統(tǒng)的點(diǎn)燃耗方程組不僅數(shù)量龐大，而且剛性很強(qiáng)，給點(diǎn)燃耗計(jì)算方法帶來困難。

求解點(diǎn)燃耗方程的方法主要包括3類：常微分方程組差分方法、線性子鏈分析（Transmutation Trajectory Analysis， TTA）方法和矩陣指數(shù)方法。

由于點(diǎn)燃耗方程剛性很強(qiáng)，因此，常微分方程組差分方法需要采用特殊迭代技巧才能獲得較好的收斂精度與計(jì)算效率。例如：MC21燃耗模塊采用變系數(shù)向后差分方法；FISPACT點(diǎn)燃耗程序采用指數(shù)歐拉差分方法，并通過平衡假設(shè)處理短半衰期核素。

在1910年，Bateman提出TTA方法并將它用于衰變鏈的計(jì)算。CINDER點(diǎn)燃耗程序和MCB等蒙卡程序的燃耗模塊將TTA方法作為核心算法。BALA K A等人[1]基于TTA方法開發(fā)了CNUCTRAN點(diǎn)燃耗程序，其具有較高的計(jì)算精度。

ORIGEN程序采用的泰勒展開法屬于多項(xiàng)式逼近方法，但它需要預(yù)先將短半衰期核素從燃耗系數(shù)矩陣中移出以降低矩陣病態(tài)性。ZHAO Z Y等人[2]基于矩陣指數(shù)方法開發(fā)了適用于加速器驅(qū)動的次臨界裝置，其中，矩陣指數(shù)采用Chebyshev有理近似方法（Chebyshev Rational Approximation Method，CRAM）求解。ZHAO Z L[3]分析了組件程序Dragon5對VVER堆型燃耗的適用性，其中，燃耗計(jì)算是關(guān)鍵因素。

譜虧損校正法（Spectral Deferred Correction， SDC）[4]是一種L穩(wěn)定的微分方程計(jì)算方法，具有高精度、高穩(wěn)定性的特點(diǎn)。基于隱式Euler的SDC方法有很好的穩(wěn)定性和收斂速率[5]，目前，該方法已經(jīng)被廣泛應(yīng)用到各個領(lǐng)域，如點(diǎn)堆動力學(xué)方程[6]、聚變堆[7]等。

1" 點(diǎn)燃耗方程

1.1" 點(diǎn)燃耗方程的一般形式

對于燃耗鏈中的每一個核素，都可以寫出其與時間相關(guān)的點(diǎn)燃耗方程：

式（1）中：t為時刻；ni為第i個核素的核素密度；為第i個核素的有效衰變常數(shù)；為第i個核素反應(yīng)生成第j個核素的有效份額。和分別由下式計(jì)算得到：

式（2）中：為第i個核素的衰變常數(shù)；為中子通量；為第i個核素反應(yīng)生成第j個核素的反應(yīng)截面；為第i個核素衰變生成第j個核素的份額。在實(shí)際燃耗計(jì)算過程中，通常假設(shè)每個燃耗步長內(nèi)的截面、通量等參數(shù)為常量，則可將點(diǎn)燃耗方程當(dāng)作一階線性常微分方程處理。

將點(diǎn)燃耗方程組寫成矩陣形式，得到以下更加簡潔的表達(dá)式：

式（3）中：矩陣A為N階點(diǎn)燃耗方程系數(shù)矩陣；N為核素系統(tǒng)中的核素?cái)?shù)量。在實(shí)際燃耗計(jì)算過程中，系數(shù)矩陣A在單個燃耗步長內(nèi)可視為常數(shù)矩陣。若給定初值邊界條件，則式（3）的解可以寫成含矩陣指數(shù)的形式：

為引述方便，將矩陣A稱為燃耗矩陣。

1.2" 矩陣指數(shù)法

矩陣指數(shù)方法將對點(diǎn)燃耗方程的求解轉(zhuǎn)化為對式（4）中矩陣指數(shù)eAt的計(jì)算。由于燃耗矩陣A具有大型、稀疏、剛性強(qiáng)的特點(diǎn)，數(shù)學(xué)上的大多數(shù)計(jì)算矩陣指數(shù)的方法不適用于燃耗問題。現(xiàn)有計(jì)算燃耗矩陣指數(shù)的方法主要基于數(shù)值逼近理論，其中包括有理逼近方法和多項(xiàng)式逼近方法。

式（4）給出了矩陣指數(shù)的點(diǎn)燃耗時間離散格式。在一個時間步內(nèi)[tn，tn+1]，點(diǎn)燃耗方程的解可以寫為

式（5）中：h為時間步長。矩陣指數(shù)方法的核心在于如何計(jì)算指數(shù)矩陣。

通常，矩陣指數(shù)函數(shù)eA采用有理函數(shù)或Taylor展開。矩陣指數(shù)函數(shù)的一般有理函數(shù)形式如式（6）所示，又稱為Padé近似。展開系數(shù)ak和bl，可以通過Taylor使與ex的展開系數(shù)盡可能的保持一致。因此，當(dāng)時間步足夠小時，Padé近似收斂到式（5）。常見的低階Padé近似有式（7）形式。

但由于A的剛性非常強(qiáng)，Padé近似并不能獲得較好的效果。但A的特征值通常為負(fù)數(shù)，因此，可以采用CRAM方法。CRAM方法可通過如下最佳逼近問題獲得：尋找有理近似，使它在負(fù)實(shí)軸上和指數(shù)函數(shù)之間的最大絕對差值最小。其一般形式如式（8）所示，隨著展開項(xiàng)M增加，收斂非常快，其收斂速率為。式（8）的展開系數(shù)由Carathéodory-Fejér算法給出。

式（8）中，ck和zk分別為展開系數(shù)和極點(diǎn)。

另一種獲得指數(shù)矩陣函數(shù)的有理近似方法是通過復(fù)平面的圍道積分獲得，見式（9）、式（10），這類方法是求積組有理近似方法（Quadrature Rational Approximation Method，QRAM）。

QRAM方法的收斂速度為，比CRAM收斂慢，但展開階不低于24項(xiàng)時，也可以獲得非常高的精度。

矩陣指數(shù)函數(shù)求解方法的缺點(diǎn)在于式（5）的時間離散格式是二階精度的，只能通過減少時間步來獲得高精度。SDC方法的時間離散格式是高階的，下面將針對該方法詳細(xì)展開。

2 譜虧損修正理論模型

2.1" 積分形式的燃耗方程

SDC的基本思想是用低階方法，如Euler方法，通過修正得到高階解。該方法將微分方程變成等價(jià)的Picard積分方程，并用基于Euler方法實(shí)施虧損修正過程。非剛性問題用顯式Euler預(yù)估和修正，剛性問題則用隱式Euler預(yù)估和修正。數(shù)值實(shí)驗(yàn)表明SDC方法可以獲得任意高的精度。

式（3）的等價(jià)Picard積分方程為

式（11）中，時刻a的解記作。假設(shè)有了一個近似解，那么原方程的殘差函數(shù)為

式（12）中，k表示第k次迭代的解。

記真實(shí)解與近似解之間的誤差為

那么，誤差函數(shù)滿足如下方程

式（12）～式（14）構(gòu)成了求解燃耗方程的關(guān)鍵。利用高精度的求積公式可以求解式（12），式（14）則采用低階方法求解。通過交替求解式（12）和式（14），解的精度可以得到提高。

2.2" 基于隱式Euler方法的SDC

在求解過程中，首先，將整個求解時間段分為小的時間步，記某一個時間步的區(qū)間為[a，b]。其次，將[a，b] 細(xì)分為s段子時間步，記作 。這里的子時間節(jié)點(diǎn)采用Legendre-Gauss-Radau II-A節(jié)點(diǎn)。在時間步[a，b]上，采用Lagrange插值多項(xiàng)式近似函數(shù)，那么，式（12）可以進(jìn)一步簡化為

式（15）中，L是Lagrange插值基函數(shù)。

誤差函數(shù)和初始的近似解在時間節(jié)點(diǎn)上的值采用隱式Euler方法求解。那么，式（14）可以寫為

最后解更新為

一個時間步內(nèi)的計(jì)算流程如下：

（1）采用隱式Euler方法，估計(jì)每個子步上的核素濃度

（2）Do k=1，…，K

①根據(jù)式（15）計(jì)算得到殘差函數(shù)

②根據(jù)式（16）計(jì)算誤差函數(shù)

③根據(jù)式（17）更新解

EndDo

該方法不僅可以用來求解線性問題，也可以用來求解非線性問題，其差別在于計(jì)算誤差函數(shù)和初始近似解的求解變成一個非線性問題，而不是一個線性問題。

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

考慮到燃耗系數(shù)矩陣的稀疏性，本文采用稀疏矩陣存儲方式存儲矩陣A。

3.1" 237Np純衰變計(jì)算

使用SDC和ORIGEN-2程序計(jì)算237Np衰變106年后所有主要產(chǎn)物的核密度。為方便比較，兩個程序均使用ORIGEN-2的衰變數(shù)據(jù)庫。

表1給出了SDC和CRAM對237Np衰變問題的計(jì)算結(jié)果對比，可以看到，SDC與CRAM核素密度計(jì)算值的相對誤差基本小于1×10-6%，即使對于半衰期為3.72×10-6秒的213Po，計(jì)算結(jié)果也符合很好。實(shí)際上，對于該問題，共花費(fèi)10步且采用定步長。可以看到，SDC計(jì)算效率和計(jì)算精度都比較好。

3.2" UO2點(diǎn)燃耗計(jì)算

使用SDC和CRAM分別計(jì)算以下含中子通量的燃耗問題：UO2以3.0×1014 cm-2s-1的通量經(jīng)10個燃耗步焚燒1年，再經(jīng)10個燃耗步衰變3年。UO2的初始成分分別為：235U， 0.03；238U， 0.97；16O， 2，單位為mol。

表2給出了部分重要錒系核素和裂變產(chǎn)物核素密度的計(jì)算結(jié)果。SDC給出的結(jié)果與CRAM使用矩陣指數(shù)方法給出的計(jì)算結(jié)果符合很好，并且與CRAM給出的計(jì)算結(jié)果符合很好：基本上，核素的相對誤差低于1×10-5 %。

由于該問題是先焚燒再衰變，因此，在建模時，將該問題分成兩個連續(xù)的問題：焚燒和衰變。每個過程都是采用10個定步長求解，可以看到，計(jì)算結(jié)果仍然符合很好。

3.3" 通量線性變化

構(gòu)造通量變化的算例，這樣，式（3）中A矩陣是隨時間線性變化的。計(jì)算條件和第3.2節(jié)一樣，除了通量從3.0×1014 n/（cm2.s）下降到3.0×1013 n/（cm2.s），且不考慮純衰變過程。用精細(xì)的CRAM算法得到的結(jié)果（等分為40步）作為參考結(jié)果，相應(yīng)結(jié)果見表3。將CRAM分為M步，可以看到步長減小1倍（N增加一倍），CRAM相對誤差就減小大約1倍。這表明CRAM具有2階精度。SDC選用定步長且一共用10個步長，可以看到，定步長的SDC用比較少的步數(shù)目就可以實(shí)現(xiàn)高精度。

3.4 參數(shù)s和K的影響

這里，將通過一個數(shù)值實(shí)驗(yàn)證明SDC的收斂速度。選取K=s，這是由于SDC的收斂速度為min（K，s）。這里選擇第3.3節(jié)的算例，保持10個步長不變。在其他計(jì)算條件不變的情況下，圖1給出了SDC相對于參考解的誤差百分比隨參數(shù)s的變化。從圖1可以看到，基本上每當(dāng)s增加1時，計(jì)算精度就提高一個量級。這表明了SDC的收斂速度是相當(dāng)快的。

4 結(jié)論

本文研究了用于復(fù)雜核素系統(tǒng)點(diǎn)燃耗計(jì)算的SDC。數(shù)值驗(yàn)證結(jié)果表明，SDC計(jì)算精度高、計(jì)算效率好，能夠準(zhǔn)確高效地計(jì)算復(fù)雜核素系統(tǒng)中的各類點(diǎn)燃耗問題。

參考文獻(xiàn)

BALA K A ， OMAR M R， SOO J Y H，et al. CNUCTRAN： A program for computing final nuclide concentrations using a direct simulation approach[J].Computer Physics Communications， 2024， 302（1）： 109258.

ZHAO Z Y， YANG Y W， GAO Q Y.Development and validation of burn-up calculation code IMPC-Burnup2.0 for accelerator-driven sub-critical system[J].Computer Physics Communications， 2021， 261（3）： 107343.

ZHAO Z L， GUO Y P， ZOU Y S，et al.Validation and application of the Dragon5 lattice code for neutronics and burnup analysis of VVER-1000 pin cell and assembly model[J].Nuclear Engineering and Design， 2023，407（12）：112279.

YAO L， XIA Y H， XU Y.L-stable spectral deferred correction methods and applications to phase field models[J]. Applied Numerical Mathematics， 2024， 197（6）：288-306.

CAUSLEY M F， SEAL D C.On the convergence of spectral deferred correction methods[J]. Communications in Applied Mathematics and Computational Science， 2019， 14（1）： 33-64.

SHEN Q C， KOCHUNAS B.High-order accurate solutions of the point kinetics equationswith the spectral deferred correction method [J]. Nuclear Science and Engineering， 2023， 197（7）， 1364-1385.

TRETIAK K， BUCHANAN J， AKERS R， et al. Performance of the BGSDC integrator for computing fast ion trajectories in nuclear fusion reactors [J]. Computer Physics Communications， 2021， 264（4）， 107876.