張彬航，郝麗娟，葛 鵬，宋 婧，何 鵬

(1.中國科學院核能安全技術研究所，中國科學院中子輸運理論與輻射安全重點實驗室，安徽 合肥 230031；2.中國科學技術大學，安徽 合肥 230027)

基于切比雪夫有理逼近和矩陣自適應降階的活化計算方法

張彬航1，2，郝麗娟1，葛 鵬1，宋 婧1，何 鵬1

(1.中國科學院核能安全技術研究所，中國科學院中子輸運理論與輻射安全重點實驗室，安徽 合肥 230031；2.中國科學技術大學，安徽 合肥 230027)

在核反應堆運行過程中，產(chǎn)生的大量中子對結構材料、回路中的腐蝕產(chǎn)物有很強的活化作用，從而對工作人員在運行、檢修、退役等多個環(huán)節(jié)造成輻射危害。因此高精度、高效率的活化計算在反應堆設計和安全分析研究中有著重要的作用。切比雪夫有理逼近方法(Chebyshev rational approximation method，CRAM) 相比于傳統(tǒng)的活化計算方法，不需要單獨處理短壽命核素，具有計算速度快、精度高、步長包容性好的優(yōu)點。本文基于超級蒙卡核模擬軟件系統(tǒng)SuperMC，采用切比雪夫有理逼近方法，并發(fā)展了基于深度優(yōu)先搜索的活化鏈動態(tài)構建方法和大規(guī)模矩陣自適應降階方法，進行了活化計算的初步研究。通過選用IAEA-ACB(International Atomic Energy Agency-Activation Calculation Benchmark)國際基準例題及壓水堆燃料包殼材料例題，初步驗證了該活化計算方法的正確性，且測試結果表明，本文發(fā)展的大規(guī)模矩陣自適應降階方法能夠有效的提高活化系數(shù)矩陣的求解效率。

活化計算；切比雪夫有理逼近方法；SuperMC；活化鏈

反應堆在運行過程中，堆芯內部產(chǎn)生的大量中子會對結構部件、回路中的腐蝕產(chǎn)物有很強的活化作用，所形成的放射性活化產(chǎn)物是反應堆運行及檢修人員職業(yè)照射的主要來源。因此，精確計算活化問題對于反應堆的安全分析、輻射防護和退役方案的制訂具有重要意義。

活化計算主要關心的是材料內的輕核素和中等質量核素受中子輻照后引起的放射性情況。在實際活化計算中，通常涉及的核素種類及反應類型繁多，核素間的轉換關系復雜，同時會產(chǎn)生大量的短壽命核素，給活化計算帶來困難。目前活化計算方法大致分為解析方法和數(shù)值計算方法兩類[1,2]。解析方法是將復雜的核素系統(tǒng)轉換為一系列的線性子鏈進行求解，計算速度慢，同時需要引入截斷誤差。數(shù)值計算方法主要包含常微分方程差分解法和矩陣指數(shù)法。常微分方程組差分解法對時間步長包容性差，矩陣指數(shù)法需要對短壽命核素進行特殊處理，影響計算精度。近年來發(fā)展的切比雪夫有理逼近方法(Chebyshev rational approxima-tion method，CRAM)[3]，作為一種新的矩陣指數(shù)法，具有計算速度快，精度高、步長包容性好等特點，且不需要單獨處理短壽命核素。

本文基于超級蒙卡核模擬軟件系統(tǒng)SuperMC(Super Monte Carlo Simulation Program for Nuclear and Radiation Process)[4-6]，采用切比雪夫有理逼近方法，進行了活化計算的研究。SuperMC是FDS團隊自主研發(fā)的通用、智能、精準的核系統(tǒng)設計與安全評價軟件，可用于反應堆的設計與分析[7-11]，輻照屏蔽分析[12-16]等領域。SuperMC目前已通過2000多個國際基準模型與實驗的校驗[17-20]。

1 活化方程與CRAM方法

(1)

式中，φ(t)表示t時刻的中子通量密度，λi表示核素i的衰變常數(shù)，σi,k表示核素i生成核素k單群截面，bj,i表示核素j生成核素i的分支比，yi,j,k表示核素j發(fā)生k反應生成核素i的產(chǎn)額。從而問題簡化為關于單一變量t的一階線性微分方程，描述了核反應系統(tǒng)中核素隨時間的變化規(guī)律。

若定義矩陣元：

(2)

則可將一階線性微分方程組轉換為矩陣形式并得到其指數(shù)解的形式如下：

(3)

n(t)=eAtn(0)

(4)

式中，n(t)表示t時刻所有核素密度組成的向量，n(0)表示t=0時刻的初始值，A為系數(shù)矩陣。

(5)

(6)

式中，I為單位矩陣，k為展開階數(shù)，通過求解線性方程組(6)得到t時刻核核子密度。

2 活化計算中的關鍵問題研究

在活化計算中，首先建立活化計算的物理模型，包括活化計算涉及的核素、活化反應鏈及反應率等；其次，采用切比雪夫有理逼近方法對物理模型進行數(shù)值求解；最后，根據(jù)數(shù)值求解得到的核子密度進行活化物理量的計算，包括：活度、余熱、接觸劑量、衰變光子譜等。本文針對活化計算中活化鏈構建方式、求解效率等關鍵問題開展了研究。

2.1 基于深度優(yōu)先搜索的活化鏈動態(tài)構建方法

在建立活化計算的物理模型時，由于初始材料核素信息通常不包括活化產(chǎn)物的核素信息，因此需要充分考慮材料中核素在中子輻照過程中產(chǎn)生的新核素。FISPACT[21]、ORIGEN等計算程序中采取預先定義的中子活化鏈進行活化產(chǎn)物添加，但往往覆蓋的核素不夠全面，擴展性較差。

由于活化涉及的核素種類繁多且變化復雜，需要尋求一種既能覆蓋所有核素反應路徑又能根據(jù)實際問題自動搜索活化產(chǎn)物的算法。本文基于深度優(yōu)先搜索算法，發(fā)展了活化鏈動態(tài)構建方法。從核素的中子反應信息和衰變信息出發(fā)，基于深度優(yōu)先搜索算法跟蹤初始核素及其子核的所有反應路徑，獲得當前問題的解空間，同時在搜索過程中引入截斷條件(反應截面、重復核素等)來避免無效搜索，提高搜索效率。當約束條件判斷當前反應鏈終止時，則搜索返回上級子核素開始下一條反應鏈的搜索，因此在搜索過程就形成了樹狀結構，如圖1所示，從而完成系統(tǒng)內所有核素的添加和活化反應鏈的動態(tài)構建。該方法能夠根據(jù)初始計算條件自動搜索活化產(chǎn)物并構建活化鏈，具有擴展性強、精度高的特點。

圖1 活化鏈動態(tài)構建示意圖Fig.1 Schematic diagram of dynamic construction of activation chain

2.2 大規(guī)模矩陣自適應降階方法

在利用切比雪夫有理逼近方法對活化方程進行數(shù)值求解時，大量的短壽命核素會導致系數(shù)矩陣規(guī)模大，剛性強，直接影響數(shù)值求解的效率和精度?？紤]到實際的輻照中子能譜，活化鏈中的部分活化反應可能不會發(fā)生，從而對應的活化產(chǎn)物產(chǎn)額為零。這些產(chǎn)額為零的核素會增加系數(shù)矩陣的規(guī)模，使得求解效率和精度降低。本文發(fā)展了大規(guī)模矩陣自適應降階方法，有效降低了系數(shù)矩陣的規(guī)模和剛性，提高活化計算的求解效率。

在活化計算中，系數(shù)矩陣A如圖2所示，對角線矩陣元代表核素的總反應消失率，同時對于A中的任意一行，除去對角線矩陣元以外，其余矩陣元分別代表其他行核素生成該核素的反應率。在確定系統(tǒng)內核素與各矩陣元對應關系的條件下，本方法依次對A中的每一個核素的總反應產(chǎn)生率進行判斷，若該核素的總產(chǎn)生率為零，則認為該核素對于系統(tǒng)中其他核素的生成沒有貢獻。同時對該核素進行標記，并將該核素對其他核素的產(chǎn)生反應率置為零。若標記核素的總個數(shù)為n，則由該方法處理后得到的系數(shù)矩陣規(guī)模由原先的m階降低為(m-n)階，從而降低了矩陣規(guī)模，提高了活化計算的求解效率。

圖2 大規(guī)模矩陣自適應降階方法示意圖

Fig.2 Schematic diagram of the mass matrix adaptive order reduction method

3 測試與驗證

為了驗證所發(fā)展的活化計算方法的正確性，本文選用IAEA-ACB國際基準例題[22]以及壓水堆包殼材料例題[23]進行了測試，選取FISPACT-II作為參考程序。

3.1 IAEA-ACB國際基準例題

IAEA-ACB國際基準例題旨在驗證活化計算程序的四個方面：(1)正確讀取活化數(shù)據(jù)庫的能力；(2)在多步長的計算中，活化程序的計算結果與基準值偏差不超過5%；(3)正確處理輕質量核素(H、He等)；(4)正確處理數(shù)據(jù)庫中的同質異能核素。該國際基準例題包括Fe和Cr兩個子例題，覆蓋了常見的中子活化反應和衰變反應類型。根據(jù)初始條件，F(xiàn)e的質量為1kg，Cr的核子密度為1.0E+20cm-3，輻照時間均為一年，中子通量譜參見文獻[22]。本文的活化計算結果與基準值的相對誤差參見圖3和圖4。

圖3 Fe的活化產(chǎn)物核子密度相對誤差Fig.3 The relative difference of activation of Fe calculation results

圖4 Cr的活化產(chǎn)物核子密度相對誤差Fig.4 The relative difference of activation of Cr calculation results

從圖3和圖4可知，由活化鏈動態(tài)構建方法得到的核素種類與基準值列出的核素種類一致。對于Fe和Cr兩個子例題，經(jīng)大規(guī)模矩陣自適應降階方法處理后的矩陣求解效率相對于處理前的矩陣求解效率分別提高了1.82倍和2.11倍。同時本文在同質異能核素(52mMn、58mMn、44mSc、45mSc)、輕質量核素(1H、2H、3H等)、中等質量核素(Fe、Cr、Mn、Ni等)等活化計算中常見核素的計算結果與基準值的相對誤差均低于0.5%，驗證了本文發(fā)展的大規(guī)模矩陣自適應降階方法的精確性。59Co和59Ni的計算結果與基準值相對偏差較大，分別為0.75%和0.71%，但中子活化過程中這兩種核素的密度非常小，與初始核素密度相差1015～1018，對于實際工程中的影響非常小?？傮w來看，本文的計算結果與基準值的相對偏差遠低于基準例題中要求的5%，滿足基準例題中對活化計算程序的精度要求。

3.2 壓水堆燃料包殼材料例題

壓水堆燃料包殼材料是反應堆的重要材料之一，其材料的活化水平直接影響到燃料包殼的壽命與堆芯安全設計。本文以壓水堆燃料包殼材料鋯合金N18為例，首先通過蒙卡輸運計算得到燃料包殼的中子通量分布，然后開展N18材料的活化計算分析與驗證。N18材料初始元素組成如表1所示，輻照方案參見文獻[23]。

表1 新型鋯合金N18的元素組成Table 1 The composition of new zirconium alloy N18

圖5給出了本文和FISPACT-II分別計算包殼材料在停堆后12個冷卻時段的活度分布，其計算結果隨時間的衰減趨勢完全一致。同時經(jīng)大規(guī)模矩陣自適應降階方法處理后的矩陣求解效率相對于處理前的矩陣求解效率提高了1.56倍。在停堆初期，包殼材料的活度約為1015Bq。在停堆后的一個月內，包殼材料的活度衰減緩慢，未出現(xiàn)量級變化，在此期間，包殼材料的活度主要受到核素89mY(T1/2=1.56e×10-1s)和89Zr(T1/2=2.82×105s)的影響。在停堆一個月至104年內，包殼材料的活度依次受到核素117mSn(T1/2=1.17×106s)、119mSn(T1/2=2.53×107s)和117mSn(T1/2=1.17×109s)的影響而出現(xiàn)量級衰減。停堆104年后，包殼材料的活度衰減程度趨于平緩。由圖6可知，本文的計算結果在12個冷卻時段與FISPACT-II的相對偏差均低于0.04%，驗證了本文發(fā)展的活化計算方法的正確性。

圖5 新型鋯合金N18活度計算結果Fig.5 The activity results of new zirconium alloy N18

圖6 新型鋯合金N18活度計算結果相對誤差Fig.6 The relative difference of new zirconium alloy N18 calculation results

4 結語

本文基于超級蒙卡核模擬軟件系統(tǒng)SuperMC，采用切比雪夫有理逼近方法，進行了活化計算的研究，并對活化鏈構建等關鍵問題的處理方法進行了優(yōu)化，發(fā)展了基于深度優(yōu)先搜索的活化鏈動態(tài)構建方法和大規(guī)模矩陣自適應降階方法。通過對IAEA-ACB國際基準例題及壓水堆燃料包殼材料例題驗證，初步證明了本文所實現(xiàn)活化計算方法的有效性與正確性，并且測試結果表明，本文發(fā)展的大規(guī)模矩陣自適應降階方法能夠有效提高矩陣的求解效率。本文活化計算方法的更多驗證與應用研究將在后續(xù)工作中開展。

致謝

本文開展研究工作中，得到了FDS團隊其他成員的大力幫助和支持，在此深表感謝！

[1] Cetnar J.General solution of Baterman equations for nuclear transmutations[J].Annals of nuclear Energy，2006，33：640-645.DOI：10.1016/j.anucene.2006.02.004.

[2] Croff A G.A user’s manual for ORIGEN2 computer code[R].Oak Ridge National Laboratory，1980.DOI：10.2172/5285077.

[3] Pusa M，Leppanen J.Computing the matrix exponential in burnup calculations[J].Nuclear science and engineering，2010，164(2)：140-150.

[4] Wu Y，Song J，Zheng H，et al.CAD-based Monte Carlo program for integrated simulation of nuclear system SuperMC[J].Annals of Nuclear Energy，2015，82(1)：161-168.DOI：10.1016/j.anucene.2014.08.058.

[5] Wu Y，F(xiàn)DS Team.CAD-Based Interface Programs for Fusion Neutron Transport Simulation.Fusion Engineering and Design，2009，84(7-11)：1987-1992.DOI：10.1016/j.fusengdes.2008.12.041.

[6] Y.Li，L.Lu，A.Ding，H.Hu，Q.Zeng，S.Zheng，Wu Y.Benchmarking of MCAM4.0 with the ITER 3D Model.Fusion Engineering and Design，2007，82：2861-2866.DOI：10.1016/ j.fusengdes.2007.02.022.

[7] Wu Y，F(xiàn)DS Team.Conceptual Design Activities of FDS Series Fusion Power Plants in China.Fusion Engineering and Design，2006，81(23-24)：2713-2718.DOI：10.1016/j.fusengdes.2006.07.068.

[8] L.Qiu，Wu Y，B.Xiao，Q.Xu，Q.Huang，B.Wu，Y.Chen，W.Xu，Y.Chen，X.Liu.A Low Aspect Ratio Tokamak Transmutation System.Nuclear Fusion，2000，40：629-633.DOI：10.1088/0029-5515/40/3Y/325.

[9] Wu Y，J.Jiang，M.Wang，M.Jin，F(xiàn)DS Team.A Fusion-Driven Subcritical System Concept Based on Viable Technologies.Nuclear Fusion，2011，51(10)：103036.1-7.DOI：10.1088/0029-5515/ 51/10/103036.

[10] Wu Y，F(xiàn)DS Team.Conceptual Design and Testing Strategy of a Dual Functional Lithium-Lead Test Blanket Module in ITER and EAST.Nuclear Fusion，2007，47(11)：1533-1539.DOI：10.7538/yzk.2015.49.S0.0007.

[11] Wu Y，F(xiàn)DS Team.Fusion-Based Hydrogen Production Reactor and Its Material Selection.Journal of Nuclear Materials，2009，386-388：122-126.DOI：10.1016/j.jnucmat.2008.12.075.

[12] Wu Y，F(xiàn)DS Team.Design Status and Development Strategy of China Liquid Lithium-Lead Blankets and Related Material Technology，2007，367-370：1410-1415.DOI：10.1016/j.jnucmat.2007.04-031.

[13] Wu Y，J.Qian，J.Yu.The Fusion-Driven Hybrid System and Its Material Selection.Journal of Nuclear Materials，2002，307-311：1629-1636.DOI：10.1016/S0022-3115(02)01272-2.

[14] Q.Huang，Wu Y，J.Li，et al.Status and Strategy of Fusion Materials Development in China.Journal of Nuclear Materials，2009，386-388：400-404.DOI：10.1016/j.jnucmat.2008.12.158.

[15] Q.Huang，J.Li，Y.Chen.Study of Irradiation Effects in China Low Activation Martensitic Steel CLAM.Journal of Nuclear Materials，2004，329：268-272.DOI：10.1016/j.jnucmat.2001.04.056.

[16] Y.Li，Q.Huang，Wu Y，T.Nagasaka，T. Muroga. Mechanical Properties and Microstructures of China Low Activation Martensitic Steel Compared with JLF-1.Journal of Nuclear Materials，2007，367-370：117-121.DOI：10.1016/j.jnucmat.2007.03.012.

[17] Wu Y，F(xiàn)DS Team.Design Analysis of the China Dual-Functional Lithium Lead (DFLL) Test Blanket Module in ITER.Fusion Engineering andDesign，2007，82：1893-1903.DOI：10.1016/j.fusengdes.2007.08.012.

[18] Wu Y，Xie Z，F(xiàn)ischer U.A discrete ordinates nodal method for one-dimensional neutron transport calculation in curvilinear geometries[J].Nuclear Science and Engineering，1999，133(3)：350-357.DOI：10.1016/j.anucene.2013.10.018.

[19] Wu Y，Zheng S，Zhu X，Wang W，et al.Conceptual Design of the Fusion-Driven Subcritical System FDS-I[J].Fusion Engineering and Design，2006，81 (Part B)：1305-1311.DOI：10.1016/j.fusengdes.2005.10.015.

[20] Wu Y，F(xiàn)DS Team.Conceptual design of the china fusion power plant FDS-II[J].Fusion Engineering and Design，2008，83(1)：1683-1689.DOI：10.1016/j.fusengdes.2008.06.048.

[21] Jean-Christphe C.Sublet，et al.The FISPACT-II User Manual[R].Culham Centre Fusion Energy，UK，2014.

[22] E.T.Cheng，et al.Report on the second international activation calculation benchmark comparison study[R].Austria，1994.

[23] Han Wenjing.Numerical Study and Preliminary Application of Activation Method Based on CRAM [D].North China Electric Power University，2015.

Researchand Verification of Activation Calculations Basedon Chebyshev Rational Approximation Method

ZHANG Bin-hang1，2，HAO Li-juan1，GE Peng1，SONG Jing1，PENG He1

(1.Institute of Nuclear Energy Safety Technology，Key Laboratory of Neutronics and Radiation Safety，
Chinese Academy of Sciences，Hefei of Anhui Prov.230031，China，2.University of Science & Technology of China，Hefei of Anhui Prov.230027，China)

A large number of neutrons are produced in nuclear reactor operation.These neutrons cause strong activation on structural materials and corrosion products，which is main radiation source on persons working on operation，maintenance and decommission.It is essential to develop activation calculation with high precision and high efficiency for the reactor design and safety analysis.Compared with the traditional activation calculation algorithm，Chebyshev rational approximation method has the advantages of computing efficiency，high precision，longtime step calculation and no need to deal with short-lived nuclides separately.In this paper，based on SuperMC，Cheyshev rational approximation method was adopted to study activation calculation，and dynamic construction method of activation chain and adaptive order reduction on massive matrix method were developed.Benchmark cases of IAEA-Activation Calculation Benchmark and PWR-Cladding Materials Case were used to validate the method.The calculation results have good agreements with reference values，which demonstrated the correctness of the activation calculation method of this work.In addition，the reduction on massive matrix method can effectively improve the efficiency of the matrix calculation.

Activation; Chebyshev rational approximation method; SuperMC; Activation chain;

2016-12-06

國家自然科學基金(11305203)、國家磁約束核聚變能發(fā)展研究專項(2014GB1120001)、中國科學院國防科技創(chuàng)新基金項目(CXJJ-16Q 231、CXJJ-15M037)，產(chǎn)業(yè)化基金

張彬航(1989—)，男，湖北荊州人，博士研究生，主要從事活化計算方法及程序研發(fā)工作

何 鵬：peng.he@fds.org.cn

TL329.2

A

0258-0918(2017)01-0065-08