呂煒楓 陳明亮 黃倩倩 宋 譽 洪 浩

（中廣核工程有限公司核電安全監(jiān)控技術(shù)與裝備國家重點實驗室 深圳 518172）

在壓水堆核電站功率運行期間，來自堆芯的中子將會對堆外設(shè)備和墻體材料進行活化，從而產(chǎn)生活化產(chǎn)物源項，進而對停堆期間和退役期間核電站內(nèi)工作人員產(chǎn)生輻射影響，而且隨著核電站運行年限的增長，中子活化產(chǎn)生的長壽命核素（如60Co）的活度將逐漸增加，因此有必要建立堆外設(shè)備和墻體材料內(nèi)活化產(chǎn)物源項的精細(xì)評估手段，從而為后端的屏蔽設(shè)計和廢物處理提供輸入。

美國橡樹嶺實驗室基于矩陣指數(shù)法開發(fā)了求解燃耗方程的燃料組件中裂變產(chǎn)物核素和活化產(chǎn)物核素源項的計算程序 ORIGEN-S［1］。ORIGEN-S 程序基于點燃耗進行計算，其已廣泛應(yīng)用于國內(nèi)外核電站設(shè)計。

歐洲原子能共同體和英國原子能管理局共同開發(fā)了歐洲壓水堆活化產(chǎn)物源項計算程序FISPACT-2007［2］。FISPACT-2007基于歐拉指數(shù)法開發(fā)，可計算高能中子裝置的材料活化以及壓水堆和聚變堆中材料的活化情況。其計算配套的核數(shù)據(jù)庫為EAF-2007。

華北電力大學(xué)韓文靜等［3-5］分別基于歐拉指數(shù)法和切比雪夫有理近似法（Chebyshev Rational Approximation Method，CRAM）對燃耗方程進行求解，采用EAF-2007數(shù)據(jù)庫開發(fā)了活化產(chǎn)物源項計算程序EuACT和ABURN。其中，EuACT程序為基于歐拉指數(shù)法開發(fā)，ABURN程序為基于CRAM方法開發(fā)。其基于ABURN程序計算了聚變堆第一壁活化例題和UO2燃料燃耗例題，并將ABURN程序的計算結(jié)果與歐洲活化產(chǎn)物源項計算程序FISPACT-2007進行對比。結(jié)果表明：ABURN程序可達到FISPACT-2007程序同等精度，并且由于采用了CRAM方法，程序在燃耗步設(shè)置方面具有高度的靈活性，初步驗證了ABURN程序的可用性與準(zhǔn)確性。

中國核動力研究設(shè)計院胡建軍等［6］采用二維離散縱標(biāo)輸運程序DORT及ENDF/B-VI庫計算壓水堆壓力容器內(nèi)中子注量率分布，用自行研制的活化產(chǎn)物源項計算程序計算冷卻劑16N和17N源項，并驗證堆芯徑向與軸向功率分布對計算結(jié)果的影響。其自行研制的活化源項計算程序名稱為SCACT。

中國科學(xué)院核能安全技術(shù)研究所張彬航等［7］基于超級蒙特卡羅核模擬軟件系統(tǒng)SuperMC，采用CRAM方法發(fā)展了基于深度優(yōu)先搜索的活化鏈動態(tài)構(gòu)建方法和大規(guī)模矩陣自適應(yīng)降階方法，進行了活化產(chǎn)物源項計算的初步研究。通過選用IAEA-ACB（InternationalAtomic Energy Agency-Activation Calculation Benchmark）國際基準(zhǔn)例題及壓水堆燃料包殼材料例題，初步驗證了該活化產(chǎn)物源項計算方法的正確性，且測試結(jié)果表明發(fā)展的大規(guī)模矩陣自適應(yīng)降階方法能夠有效地提高活化系數(shù)矩陣的求解效率。

綜上所述，國內(nèi)外各研究機構(gòu)針對壓水堆核電站活化產(chǎn)物源項計算程序進行了大量的研究工作，但因程序與工程設(shè)計脫節(jié)以及軟件認(rèn)證的問題，除了國外開發(fā)的相關(guān)程序，尚未有自主研發(fā)的活化產(chǎn)物源項計算程序在工程設(shè)計中大規(guī)模應(yīng)用。在華龍一號“走出去”的大背景下，有必要研發(fā)可應(yīng)用于工程實際的具有自主知識產(chǎn)權(quán)的壓水堆活化產(chǎn)物源項計算程序。

1 核電站內(nèi)中子活化產(chǎn)物源項計算模型

1.1 計算方程組的建立

核電站內(nèi)受中子輻照的非裂變材料內(nèi)放射性核素產(chǎn)生途徑有以下兩種：先驅(qū)核素衰變產(chǎn)生和先驅(qū)核素中子俘獲產(chǎn)生。核電站內(nèi)受中子輻照的非裂變材料內(nèi)放射性核素消失途徑有以下兩種：放射性衰變消失以及與中子反應(yīng)消失。可分別建立放射性核素i的產(chǎn)生率和消失率計算公式如下：

式中：Si(t)為t時刻放射性核素i的產(chǎn)生率，s-1；Pi(t)為t時刻放射性核素i的消失率，s-1；Φk為k群中子注量率，m-2·s-1；Ni(t)為t時刻中子輻照材料內(nèi)放射性核素i的積存量；V為中子輻照材料體積，m3；nj為先驅(qū)核素 j的核子密度，m-3；σɑ，i，k為 k群中子對先驅(qū)核素i的微觀吸收截面，m2；σx，j，k為k群中子對先驅(qū)核素j的微觀x反應(yīng)截面，x反應(yīng)指（n，r）、（n，2n）、（n，p）、（n，ɑ）、（n，t）等中子反應(yīng)，m2；λi為 i核素的衰變常數(shù)，s-1；Ki，j為衰變分支比，無量綱；Γx，i，j為先驅(qū)核素 j通過與中子的x反應(yīng)產(chǎn)生核素i的產(chǎn)額，無量綱。

根據(jù)以上放射性核素產(chǎn)生率和消失率的計算公式，可建立放射性核素i積存量的微分方程如下：

求解以上微分方程可得中子輻照材料內(nèi)放射性核素i的積存量。如需計算n種核素，則以上方程可寫成矩陣的形式：

其中：

在以上方程中，非放射性核素也在n種核素之內(nèi)，由此先驅(qū)核素中子俘獲產(chǎn)生項可包含放射性核素和非放射性核素中子俘獲活化的貢獻。

對于式（3）中方程，其解為：

式中：X(t)是核素積存量矩陣；A是包含衰變和中子反應(yīng)的轉(zhuǎn)換矩陣（活化矩陣）。X(0)是初始核素積存量矩陣。

1.2 矩陣指數(shù)法求解

對式（8）中指數(shù)函數(shù)部分進行Taylor級數(shù)展開，有：

式中：I表示單位矩陣。矩陣指數(shù)法在理論上可以精確求解燃耗問題，但是在具體的例子中，求矩陣冪的精確值時會產(chǎn)生各種困難。在這種情形下，就需要借助Bateman方程或者使用Gauss-Seidel迭代法來解決。

燃耗矩陣A的規(guī)模通常在2000×2000量級，但是在實際的活化源項計算中，往往僅關(guān)注于幾個主要的核素，極少出現(xiàn)一次性計算2000×2000的情況，考慮到核電廠實際情況，這個規(guī)模將大大縮小。

1.3 CRAM方法求解

CRAM方法的核心是在[0，∞)上對e-x進行有理逼近。定義πk，l表示所有如下形式實有理函數(shù)rk，l(x)的集合：

式中：pj是不高于j階的多項式。根據(jù)逼近理論可知存在唯一的 r^k，l∈ πk，l滿足：

對于給定的k和l，可通過Remez算法或Carathéodory-Fejér方法得到逼近表達式。由式（10）可知，有理函數(shù)r^k，k的極點{θ1，θ2，...，θk}是分立極點，因此可以通過部分分式展開得到：

式中：α0是函數(shù)在無窮遠的極限；αi表示極點θi處的留數(shù)。式（12）也可以通過留數(shù)定理得到。另外，注意到αi和θi共軛出現(xiàn)，因此，式（12）可以簡化成：

對于CRAM方法計算矩陣指數(shù)，主要的難點在于對給定的k值求出相應(yīng)的有理函數(shù)的系數(shù)。原則上來說，p^k和q^k的系數(shù)都可以通過Remez型方法得到，但是計算過程需要相當(dāng)有技巧的算法和非常高精度的計算過程。文獻調(diào)研顯示14階系數(shù)的計算結(jié)果已經(jīng)具有較高的精度［5，8］。基于此，CPAP程序采用16階系數(shù)以進一步提高精度，采用的CRAM系數(shù)見表1。

2 活化產(chǎn)物源項計算程序CPAP的開發(fā)

基于以上中子活化產(chǎn)物源項計算模型，采用面向?qū)ο蟮木幊趟季S以及當(dāng)前主流的圖形程序界面，開發(fā)了具有良好人機界面的各堆型壓水堆核電站通用的中子活化源項計算程序CPAP。其中，程序圖形化界面和計算模塊均采用C++語言開發(fā)。CPAP程序的計算流程圖見圖1和圖2。

CPAP程序有如下特點：

1）各堆型壓水堆核電站在計算中子活化產(chǎn)物源項時，一般采用矩陣指數(shù)法或者CRAM方法，由此CPAP程序中提供了矩陣指數(shù)法和CRAM方法兩種選項；

2）CPAP程序中采用的CRAM方法基于16階系數(shù)，計算精度進一步提高；

3）CPAP程序采用的中子反應(yīng)多群截面數(shù)據(jù)庫基于ENDF/B-VII.0中的中子反應(yīng)子庫（NSUB=10）及EAF-2010中子活化反應(yīng)數(shù)據(jù)庫加工得出，涉及860種核素；

4）CPAP程序的截面數(shù)據(jù)庫采用47群結(jié)構(gòu)，能群劃分參照BUGLE-96多群數(shù)據(jù)庫的能群邊界；加工多群截面數(shù)據(jù)庫時，采用傳統(tǒng)核截面加工處理軟件NJOY中的相關(guān)運行模塊，最終制成多群數(shù)據(jù)文件時采用壓水堆堆芯外的中子能譜作為權(quán)重譜；截面數(shù)據(jù)庫涉及的反應(yīng)類型共11種，包括總反應(yīng)、彈性散射反應(yīng)、非彈性散射反應(yīng)、輻射俘獲反應(yīng)、（n，2n）、（n，3n）、（n，α）、（n，p）、（n，t）、（n，2p）和裂變反應(yīng)；

表116 階CRAM系數(shù)Table 1 16th order coefficient of CRAM method

圖1 CPAP程序計算流程圖1Fig.1 Calculation flowchart 1 of CPAP code

5）具有良好的人機界面，根據(jù)用戶需求預(yù)先嵌入大量常見材料數(shù)據(jù)庫。

CPAP程序界面見圖3。

3 活化產(chǎn)物源項計算程序CPAP的測試驗證

對于活化產(chǎn)物源項計算程序CPAP的測試包括算法的測試以及數(shù)據(jù)庫的測試兩部分。以下分別說明。

3.1 算法的測試

選取了部分國際上通用的基準(zhǔn)題以及自主建立的典型例題對CPAP程序的兩種算法進行測試。算法測試時CPAP程序所使用的核素豐度、反應(yīng)類型、截面數(shù)據(jù)和衰變常數(shù)均取自EAF-2007數(shù)據(jù)庫，與FISPACT-2007程序一致。本文以兩個典型的活化算例說明程序的算法測試結(jié)果。

3.1.1 AISI304活化計算測試

在本測試算例中，AISI304的初始質(zhì)量為1 kg，其密度為8.04 g·cm-3。元素質(zhì)量成分見表2。中子通量密度為 3.641×1015cm-2·s-1，中子輻照時間為1 a。中子能譜為100群的GAM-II格式能譜［9］。

AISI304活化例題的計算結(jié)果以及其與FISPACT-2007程序計算結(jié)果之間的對比如表3所示。由于計算所得核素較多，僅選取典型核素的結(jié)果展示。

圖2 CPAP程序計算流程圖2Fig.2 Calculation flowchart 2 of CPAP code

圖3 CPAP程序界面展示Fig.3 Interface display of CPAP code

表2 AISI304的元素成分Table 2 Elemental content ofAISI304

表3中，Nf表示使用FISPACT-2007程序計算出的核素的原子數(shù)量；Nc表示使用CPAP程序中的CRAM求解器計算出的核素的原子數(shù)量；Nm表示使用CPAP程序中的矩陣指數(shù)求解器計算出的核素的 原 子 數(shù) 量 ；delta_c=（Nc-Nf）/Nf，delta_m=（Nm-Nf）/Nf。

3.1.2 混凝土活化計算測試

在本測試算例中，混凝土的初始質(zhì)量為1 kg，其密度為2.35 g·cm-3。元素質(zhì)量成分見表4。中子通量密度為3.641×1015cm-2·s-1，中子輻照時間為1 a。中子能譜為100群的GAM-II格式能譜［9］。

混凝土活化例題的計算結(jié)果以及其與FISPACT-2007程序計算結(jié)果之間的對比如表5所示。由于計算所得核素較多，僅選取典型核素的結(jié)果進行展示。

3.1.3 測試結(jié)論

經(jīng)測試，對于測試算例，CPAP程序的兩種算法在活化計算時大部分核素的計算結(jié)果與FISPACT-2007程序的計算結(jié)果之間的誤差范圍均在1%以下。但在H和He的計算時CPAP程序的計算結(jié)果與FISPACT-2007程序有較大誤差，這主要是因為FISPACT-2007程序在計算H和He時實際上是合并計算了（n，p）和（n，α）反應(yīng)產(chǎn)生的質(zhì)子和α粒子，而CPAP程序不把（n，p）和（n，α）反應(yīng)產(chǎn)生的質(zhì)子和α粒子計算為H和He元素，所以FISPACT-2007程序計算結(jié)果高于CPAP程序的計算結(jié)果。

表3 算法測試時AISI304活化例題的計算結(jié)果Table 3 Calculation result ofAISI304 activation during algorithm test

表4 混凝土的元素成分Table 4 Elemental content of concrete

3.2 截面庫和衰變數(shù)據(jù)庫的測試

3.2.1 測試方法

CPAP程序的截面庫和衰變數(shù)據(jù)庫的測試采用與ORIGEN-S程序?qū)Ρ确治龅姆绞?。具體方法為首先使用CPAP程序的截面庫和衰變數(shù)據(jù)庫替換ORIGEN-S程序自帶的截面庫和衰變數(shù)據(jù)庫，然后基于替換截面庫和衰變數(shù)據(jù)庫后的ORIGEN-S程序構(gòu)建典型算例，并與原版ORIGEN-S程序的計算結(jié)果進行對比。

由于CPAP程序的截面數(shù)據(jù)庫為47群的多群截面庫，而ORIGEN-S程序最終讀取和使用的中子反應(yīng)截面為單群截面，由此需要對CPAP程序的截面數(shù)據(jù)庫進行并群處理。并群所基于的中子能譜為使用蒙特卡羅輸運計算程序?qū)Φ湫蛪核训亩研具M行輸運計算后按照BUGLE-96庫中能群邊界統(tǒng)計得出的分為47群的中子能譜。

表5 算法測試時混凝土活化例題的計算結(jié)果Table 5 Calculation result of concrete activation during algorithm test

3.2.2 測試算例

對于CPAP程序的截面庫和衰變數(shù)據(jù)庫的測試算例見表6。其中混凝土和AISI304的元素成分分別見表2和表4。

選取各算例中典型核素的計算結(jié)果與ORIGEN-S程序原版計算結(jié)果進行對比，對比結(jié)果分別見表7和表8。由對比結(jié)果可知，部分核素存在小于10%的偏差，其可能原因為截面數(shù)據(jù)庫替換時并群處理考慮的中子譜的差異以及截面加工所基于的評價庫版本的差異。

表6 CPAP數(shù)據(jù)庫的測試算例Table 6 List of test cases for database of CPAPcode

表7和表8中，Ns表示使用ORIGEN-S程序原數(shù)據(jù)庫計算出的核素的原子密度，cm-3；Nc表示使用CPAP程序中的數(shù)據(jù)庫替ORIGEN-S程序數(shù)據(jù)庫后計算出的核素的原子數(shù)量，cm-3。delta_c=（Nc-Ns）/Ns。

表7 數(shù)據(jù)庫測試時AISI304活化例題的計算結(jié)果Table 7 Calculation result ofAISI304 activation during database test

3.2.3 測試結(jié)論

由表7和表8中測試結(jié)果可知，對于所測試的算例，使用CPAP程序的截面庫和衰變數(shù)據(jù)庫所得絕大部分核素的計算結(jié)果與ORIGEN-S程序自帶數(shù)據(jù)庫所得計算結(jié)果之間的偏差在1%以內(nèi)，少量核素存在不超過10%的偏差?？紤]到測試時截面數(shù)據(jù)庫替換過程中并群處理考慮的中子譜的差異以及截面加工所基于的評價庫版本的差異，以上偏差可認(rèn)為合理。

4 結(jié)語

本文基于堆外設(shè)備和墻體材料中活化產(chǎn)物的產(chǎn)生和遷移的機理建立了堆外設(shè)備和墻體材料中活化產(chǎn)物源項的計算模型，并編制了計算程序CPAP1.0版。基于典型材料的活化例題，采用國外同類型軟件對比計算的方式對CPAP程序的算法和數(shù)據(jù)庫分別進行了測試，測試結(jié)果顯示：CPAP程序計算精度可滿足工程設(shè)計的需要。CPAP程序具有人機界面友好以及求解器可選的優(yōu)點，可廣泛應(yīng)用于壓水堆核電站的設(shè)計、運行和退役階段，解決堆外設(shè)備和墻體材料中活化源項計算的相關(guān)問題。

表8 數(shù)據(jù)庫測試時混凝土活化例題的計算結(jié)果Table 8 Calculation result of concrete activation during database test