基于高階CRAM的燃耗程序開發(fā)與驗(yàn)證

張彬航，楊森權(quán)，陳云龍，洪 鋒,*，袁顯寶，張永紅，唐海波，馮虹瑛

(1.三峽大學(xué) 機(jī)械與動(dòng)力學(xué)院，湖北 宜昌 443002;2.中核武漢核電運(yùn)行技術(shù)股份有限公司，湖北 武漢 430074)

燃耗計(jì)算是反應(yīng)堆設(shè)計(jì)及燃料管理的重要環(huán)節(jié)，也是反應(yīng)堆燃耗信任制技術(shù)應(yīng)用及乏燃料后處理的基礎(chǔ)[1]。燃耗計(jì)算的核心是求解燃耗方程，由于堆芯燃料在中子輻照過程中產(chǎn)生大量的裂變產(chǎn)物，且不同核素間的半衰期差別極大，導(dǎo)致燃耗方程組具有大型、稀疏、剛性強(qiáng)的特點(diǎn)。目前求解的方法主要包括線性鏈解析方法、常微分方程組差分方法和矩陣指數(shù)法。線性鏈解析方法將堆內(nèi)復(fù)雜的核素轉(zhuǎn)換關(guān)系拆分成一系列線性反應(yīng)鏈后進(jìn)行求解，具有較高的計(jì)算精度，但為了平衡計(jì)算精度和效率，往往需要引入截?cái)嚅撝瞪釛壷匾暂^低的反應(yīng)鏈。同時(shí)對(duì)于深燃耗問題，已有研究表明線性鏈解析方法存在較大的計(jì)算誤差[2]。常微分方程組差分方法則是通過采取特殊的數(shù)值計(jì)算技巧處理燃耗方程組的剛性問題，計(jì)算速度較快，但計(jì)算精度對(duì)時(shí)間步長的選取較為敏感，代表程序?yàn)镕ISPACT[3]。矩陣指數(shù)法主要包括泰勒展開與截?cái)喾椒ā⑶斜妊┓蛴欣斫品椒?CRAM)等。ORIGEN程序采用泰勒展開與截?cái)喾椒ㄟM(jìn)行燃耗方程求解，具有速度快、步長包容性好的特點(diǎn)，但需對(duì)短壽命核素單獨(dú)處理后近似求解[4]。CRAM具有速度快、精度高且不需單獨(dú)對(duì)短壽命核素進(jìn)行處理的特點(diǎn)，在很多燃耗程序中得到了廣泛的發(fā)展與應(yīng)用，包括Serpent[5]、RMC[6]和JMCT[7]等。但目前已有研究表明，CRAM對(duì)于在長時(shí)間步長下核素系統(tǒng)的衰變計(jì)算精度較差，難以滿足高保真燃耗計(jì)算的發(fā)展需求[8-9]。

本文基于高階CRAM研制點(diǎn)燃耗程序ICRAM，同時(shí)內(nèi)耦合于蒙特卡羅輸運(yùn)程序OpenMC，形成燃耗計(jì)算分析程序OPICE，通過一系列基準(zhǔn)題的計(jì)算，驗(yàn)證OPICE的計(jì)算精度和準(zhǔn)確性。

1 高階切比雪夫有理近似方法

在中子輻照下，堆內(nèi)燃料中核素成分隨時(shí)間的變化可用燃耗方程組描述，其矩陣形式如下：

(1)

式中：N(t)為t時(shí)刻核素的核子密度向量；A為系數(shù)矩陣，由系統(tǒng)中所有核素的反應(yīng)率和衰變率組成。因此，可得矩陣方程(式(1))的解：

N(t)=eAtN(0)

(2)

(3)

(4)

N(t)≈α0N(0)-

(5)

式中，I為單位矩陣。

目前基于CRAM發(fā)展的大部分燃耗程序皆使用了PFD形式的有理近似展開，當(dāng)階數(shù)k取到14階或16階時(shí)就能獲得較為理想的計(jì)算結(jié)果。但如果核素核子密度在某一燃耗步長內(nèi)急劇下降，那么PFD形式的CRAM計(jì)算精度就會(huì)下降或?qū)е掠?jì)算錯(cuò)誤[12]。因此，本文基于CRAM發(fā)展了不完全局部分解(IPF)形式的高階有理近似展開方法來滿足高保真燃耗計(jì)算的精度需求。

ml=1,2,…，μ

(6)

式中：ml為低階有理函數(shù)的階數(shù)；μ為整數(shù)。

(7)

(8)

N(t)=eAtN(0)≈

(9)

基于上述高階CRAM研制了點(diǎn)燃耗程序ICRAM，其中階數(shù)k的取值為16、32和48。由式(5)和(9)可看出，在相同階數(shù)下，PFD和IPF形式的CRAM計(jì)算量相當(dāng)，且在求解過程中涉及多次系數(shù)矩陣的求逆，對(duì)發(fā)展大規(guī)模高保真燃耗計(jì)算的效率會(huì)存在一定影響，因此對(duì)ICRAM的效率優(yōu)化將是未來研究的重點(diǎn)。

2 蒙特卡羅輸運(yùn)-燃耗耦合計(jì)算方法

蒙特卡羅輸運(yùn)-燃耗耦合計(jì)算實(shí)質(zhì)上是蒙特卡羅中子輸運(yùn)計(jì)算與點(diǎn)燃耗計(jì)算過程的相互耦合。按照耦合方式的不同，分為外耦合和內(nèi)耦合兩種。外耦合方法通過第三方接口程序和外部存儲(chǔ)數(shù)據(jù)來實(shí)現(xiàn)蒙特卡羅輸運(yùn)程序和點(diǎn)燃耗程序之間的數(shù)據(jù)傳遞，導(dǎo)致計(jì)算效率低，計(jì)算規(guī)模受限。新近開發(fā)的蒙特卡羅燃耗程序多采用內(nèi)耦合的方式，即在蒙特卡羅程序中嵌入點(diǎn)燃耗計(jì)算模塊，能減少用戶工作量，提高計(jì)算效率。本文在完成點(diǎn)燃耗程序ICRAM研制的基礎(chǔ)上，進(jìn)一步將OpenMC(V0.10.0)與ICRAM進(jìn)行內(nèi)耦合開發(fā)，形成了燃耗計(jì)算程序OPICE。

OPICE的簡(jiǎn)要計(jì)算流程如圖1所示。首先通過中子輸運(yùn)計(jì)算得到各燃耗區(qū)中子通量密度、單群反應(yīng)率等中子學(xué)參數(shù)。然后經(jīng)過加工處理后傳遞給燃耗計(jì)算模塊并對(duì)各燃耗區(qū)執(zhí)行點(diǎn)燃耗計(jì)算，得到步長末各燃耗區(qū)的核素種類及核子密度。最后更新相應(yīng)燃耗區(qū)核素的核子密度并傳遞給輸運(yùn)計(jì)算進(jìn)行下一時(shí)間步長的耦合計(jì)算，依次迭代直至燃耗壽期末。OPICE是基于OpenMC進(jìn)行輸運(yùn)計(jì)算來獲得不同燃耗區(qū)的中子通量密度和單群反應(yīng)截面。而在蒙特卡羅輸運(yùn)計(jì)算中，中子通量的計(jì)數(shù)結(jié)果對(duì)源粒子進(jìn)行了歸一化處理，并不能真實(shí)反映燃耗區(qū)的中子通量密度。因此需根據(jù)當(dāng)前燃耗步長內(nèi)的熱功率進(jìn)行換算得到每個(gè)燃耗區(qū)的真實(shí)中子通量密度。假設(shè)堆芯共有N個(gè)燃耗區(qū)，每個(gè)燃耗區(qū)內(nèi)有M個(gè)核素，熱功率為P，計(jì)算公式如下：

(10)

式中：φreal為真實(shí)中子通量密度；φstatic為歸一化中子通量密度；Vj為燃耗區(qū)j的體積；Ri和Qi分別為核素i的裂變反應(yīng)率及可利用能。為進(jìn)一步提高耦合計(jì)算精度，OPICE采用了預(yù)估-校正和子步法兩種耦合策略。

圖1 OPICE燃耗計(jì)算流程圖Fig.1 Flow chart of burnup calculation of OPICE

同時(shí)OPICE還支持點(diǎn)燃耗計(jì)算功能，此時(shí)OPICE采用精細(xì)燃耗數(shù)據(jù)庫進(jìn)行點(diǎn)燃耗計(jì)算。該數(shù)據(jù)庫整合了ORIGEN-S和ORIGEN-2的最新燃耗數(shù)據(jù)庫，包含了1 487個(gè)核素、23種中子截面反應(yīng)、11種衰變反應(yīng)及30種重核的裂變產(chǎn)額。對(duì)于特定的24種錒系核素和32種裂變產(chǎn)物核素給出了精確的裂變能和中子俘獲能，其余核素則分別認(rèn)為裂變能近似為200 MeV，中子俘獲能近似為5 MeV，以滿足中子通量密度與功率的精確轉(zhuǎn)換[15]。

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

3.1 高階CRAM的數(shù)值精度

為驗(yàn)證分析高階CRAM的數(shù)值精度，測(cè)試用例選取了富集度為3.4%的UO2燃料棒進(jìn)行測(cè)試，計(jì)算該燃料棒輻照至燃耗深度為1 MW·d/kg(U)后，再冷卻100 d的核素核子密度的變化情況。參考解選用了基于MATLAB符號(hào)工具箱和高精度計(jì)算工具箱研制的TTA算法的計(jì)算結(jié)果。在OPICE中，分別使用了16、32、48階PFD和IPF兩種算法進(jìn)行計(jì)算。為保證計(jì)算效率，核素核子密度的截?cái)嘀翟O(shè)為1.0×10-30cm-3，兩種算法相對(duì)于TTA算法計(jì)算結(jié)果的相對(duì)偏差如圖2所示。

如圖2a所示，當(dāng)階數(shù)k同取16時(shí)，PFD形式的計(jì)算精度比IPF形式的計(jì)算精度稍差。隨著階數(shù)k提高至32，由圖2b可見IPF形式的計(jì)算精度相比于PFD形式有了顯著提高，大部分核素核子密度的相對(duì)偏差在10-5以內(nèi)，其中核素核子密度大于10-6cm-3的所有核素的相對(duì)偏差均在10-13以內(nèi)。當(dāng)階數(shù)k取至48時(shí)，IPF形式計(jì)算得到的核素核子密度大于10-20cm-3的所有核素的相對(duì)偏差在10-10以內(nèi)，其中大部分核素核子密度的相對(duì)偏差在10-14以內(nèi)，其精度優(yōu)勢(shì)進(jìn)一步擴(kuò)大，表明高階CRAM的計(jì)算結(jié)果是足夠精確的。值得注意的是，對(duì)于傳統(tǒng)PFD形式，隨階數(shù)的提高，計(jì)算精度并不會(huì)顯著提高。其中一個(gè)重要原因可能是在PFD形式中，相應(yīng)的留數(shù)會(huì)隨階數(shù)的提高而增大，導(dǎo)致計(jì)算結(jié)果對(duì)舍入誤差更敏感，計(jì)算精度下降[11]。

為評(píng)估高階CRAM的步長包容性，將上述測(cè)試?yán)}中UO2燃料棒的冷卻時(shí)間步長分別取至10、103和106a。OPICE分別采用了16、32、48階的PFD和IPF形式的CRAM進(jìn)行計(jì)算，核素核子密度截?cái)嘀等?.0×10-20cm-3，相對(duì)偏差及單步計(jì)算時(shí)間列于表1。由表1可見，在相同時(shí)間步長下，IPF形式的計(jì)算精度優(yōu)于PFD形式。同時(shí)隨時(shí)間步長的增加，高階CRAM的計(jì)算結(jié)果充分體現(xiàn)了其高精度性和優(yōu)異的步長包容性，單步時(shí)間步長可達(dá)到106a。注意到表1中對(duì)于48階PFD形式的大步長計(jì)算，核素相對(duì)偏差小于10-10的核素百分比顯著低于16階和32階的計(jì)算結(jié)果。其可能原因如前文所述，隨階數(shù)的提高，高階PFD形式中留數(shù)的增大會(huì)引起其計(jì)算結(jié)果對(duì)舍入誤差的敏感性增強(qiáng)，導(dǎo)致數(shù)值精度下降。同時(shí)隨時(shí)間步長的增長，更多核素的核子密度低于截?cái)嘀?，?duì)表1的統(tǒng)計(jì)結(jié)果沒有貢獻(xiàn)。后續(xù)工作還需進(jìn)一步采用更多例題進(jìn)行測(cè)試分析。在計(jì)算效率方面，單步計(jì)算時(shí)間隨階數(shù)的提高近似成倍增加。計(jì)算效率主要由兩方面因素影響：一是隨階數(shù)的提高，計(jì)算過程中系數(shù)矩陣的多次求逆導(dǎo)致了耗時(shí)略長；二是為保證計(jì)算精度，ICRAM采用了多精度浮點(diǎn)數(shù)據(jù)庫MPFR進(jìn)行浮點(diǎn)數(shù)運(yùn)算，使得計(jì)算時(shí)間增加[16]。因此計(jì)算效率優(yōu)化將是下一步研究工作重點(diǎn)。

a——16階；b——32階；c——48階圖2 IPF和PFD形式計(jì)算結(jié)果的相對(duì)偏差 Fig.2 Relative deviation of calculation result of IPF and PFD forms

表1 高階CRAM的步長包容性計(jì)算結(jié)果Table 1 Calculation result of step size tolerance of high-order CRAM

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

OECD/NEA壓水堆柵元燃耗基準(zhǔn)題來自于歐洲經(jīng)濟(jì)合作與發(fā)展組織核能署(OECD/NEA)發(fā)布的燃耗信任制系列基準(zhǔn)題之一，其計(jì)算目標(biāo)是壓水堆組件中的單個(gè)柵元，并比較不同程序系統(tǒng)對(duì)燃料中同位素成分進(jìn)行燃耗計(jì)算和分析的能力。OECD/NEA發(fā)布了16個(gè)機(jī)構(gòu)共計(jì)21套程序的計(jì)算結(jié)果，包括對(duì)乏燃料有效增殖因數(shù)有著重要影響的核素質(zhì)量等結(jié)果。該基準(zhǔn)題包含4個(gè)完整的運(yùn)行循環(huán)，各循環(huán)的燃耗時(shí)間、冷卻時(shí)間及硼濃度列于表2。根據(jù)功率和最終燃耗深度的不同，分為工況A(27.35 GW·d/tHM)、工況B(37.12 GW·d/tHM)和工況C(44.34 GW·d/tHM)。具體的幾何與材料參數(shù)詳見文獻(xiàn)[17]。

表2 運(yùn)行循環(huán)參數(shù)Table 2 Parameter of burnup operation

OPICE采用ENDF/B-Ⅶ數(shù)據(jù)庫進(jìn)行計(jì)算，表3列出工況B下核素質(zhì)量的計(jì)算結(jié)果。由表3可見，OPICE對(duì)大部分核素質(zhì)量計(jì)算結(jié)果的相對(duì)偏差在5%以內(nèi)，吻合良好。另外部分核素如237Np、149Sm與實(shí)驗(yàn)值的相對(duì)偏差較大，但通過對(duì)比可發(fā)現(xiàn)，這2個(gè)核素的計(jì)算結(jié)果落在其他程序的結(jié)果范圍內(nèi)。由表3還可看出，21套程序的計(jì)算結(jié)果存在一定范圍，主要原因是不同程序進(jìn)行輸運(yùn)計(jì)算時(shí)所使用的方法不同，對(duì)中子通量密度和反應(yīng)率的計(jì)算會(huì)造成一定的誤差。同時(shí)參考程序所使用的數(shù)據(jù)庫主要來自ENDF/B-Ⅴ、ENDF/B-Ⅵ等，在進(jìn)行燃耗計(jì)算時(shí)，大部分程序直接使用少群截面庫或ORIGEN-2的數(shù)據(jù)庫，對(duì)最終的計(jì)算結(jié)果也會(huì)造成一定的誤差?？傮w來看，OPICE具有良好的燃耗計(jì)算精度。

表3 工況B下OPICE核素質(zhì)量的計(jì)算結(jié)果Table 3 Calculation result of isotopic mass of OPICE under condition B

3.3 快堆燃耗基準(zhǔn)題

該基準(zhǔn)題是OECD/NEA發(fā)布的一個(gè)二氧化钚燃料的快堆基準(zhǔn)題，其計(jì)算目標(biāo)是1個(gè)快中子增殖堆的等效模型。該快堆R-z幾何結(jié)構(gòu)如圖3所示，包括堆芯內(nèi)部區(qū)域、堆芯外部區(qū)域、控制棒孔道區(qū)域和反射層區(qū)域，其中不同富集度的MOX燃料組件被布置在堆芯內(nèi)、外兩個(gè)區(qū)域。該堆的熱功率為1 500 MW，負(fù)載因子為0.8，具體的幾何、材料和功率參數(shù)詳見文獻(xiàn)[18]。

圖3 快堆基準(zhǔn)題R-z幾何結(jié)構(gòu)圖Fig.3 R-z geometry sketch for fast reactor benchmark

不同研究機(jī)構(gòu)發(fā)布的計(jì)算結(jié)果包括了有效增殖因數(shù)、核素核子密度等物理量隨燃耗的變化。此外從文獻(xiàn)[19]中找到了OMCB的計(jì)算結(jié)果。OPICE采用ENDF/B-Ⅶ數(shù)據(jù)庫進(jìn)行計(jì)算，在燃耗過程中設(shè)置了3個(gè)有效滿功率天(EFPD)作為重要時(shí)間點(diǎn)，分別為0、375和625 EFPD。表4列出不同機(jī)構(gòu)反應(yīng)率虧損率的計(jì)算結(jié)果。由表4可見，OPICE的計(jì)算結(jié)果落在了其他機(jī)構(gòu)發(fā)布的計(jì)算結(jié)果之間。若以發(fā)布結(jié)果的均值作為參考值，則在375 EFPD和625 EFPD時(shí)，OPICE對(duì)應(yīng)的差異為-2.28%和-1.48%，最大差異來自PNC(J3.2)的-2.81%和-3.10%。表5進(jìn)一步列出堆芯燃耗區(qū)中重要核素在燃耗周期內(nèi)核素質(zhì)量的變化情況。以各機(jī)構(gòu)的平均值作為參考值，可看出OPICE的計(jì)算結(jié)果同樣在各機(jī)構(gòu)的范圍內(nèi)，且具有較小的相對(duì)偏差，證明了OPICE對(duì)快譜燃耗系統(tǒng)具有較好的計(jì)算精度。

表4 快堆基準(zhǔn)題中初始keff和反應(yīng)率虧損率的對(duì)比Table 4 Comparison of initial keff and reactivity loss ratio for fast reactor benchmark

表5 快堆基準(zhǔn)題中核素質(zhì)量的變化對(duì)比 Table 5 Comparison of isotopic mass variation for fast reactor benchmark

4 結(jié)論

本文基于高階CRAM研制了點(diǎn)燃耗程序ICRAM，并與蒙特卡羅輸運(yùn)程序OpenMC內(nèi)耦合形成了一套燃耗計(jì)算分析程序OPICE。在點(diǎn)燃耗算法的研究中，高階IPF形式的CRAM相比于傳統(tǒng)PFD形式的CRAM，具有數(shù)值穩(wěn)定性好、計(jì)算精度高、步長包容性更好等特點(diǎn)，更符合高保真燃耗計(jì)算的需求。通過數(shù)值驗(yàn)證與分析可發(fā)現(xiàn)，高階的CRAM單步時(shí)長可達(dá)到106a，并仍可保持很高的計(jì)算精度。在耦合方式選擇中，OPICE采用內(nèi)耦合的方法避免了傳統(tǒng)外耦合導(dǎo)致計(jì)算效率低、計(jì)算規(guī)模受限等問題，并采用了預(yù)估-校正和子步法兩種耦合策略保證計(jì)算精度。選用了OECD/NEA壓水堆柵元燃耗基準(zhǔn)題和快堆燃耗基準(zhǔn)題對(duì)OPICE進(jìn)行驗(yàn)證，程序計(jì)算結(jié)果與實(shí)驗(yàn)值及各程序的參考值吻合良好，初步證明了程序的正確性和有效性。后續(xù)工作中將重點(diǎn)對(duì)高階CRAM中的系數(shù)矩陣求逆效率、數(shù)據(jù)存儲(chǔ)結(jié)構(gòu)等方面進(jìn)行優(yōu)化，同時(shí)開展大規(guī)模并行燃耗計(jì)算方法研究，以保證OPICE在大規(guī)模高保真燃耗計(jì)算中的求解效率，使其具有更加普適的工程應(yīng)用價(jià)值。