基于切比雪夫有理逼近和矩陣自適應(yīng)降階的活化計(jì)算方法

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

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

基于切比雪夫有理逼近和矩陣自適應(yīng)降階的活化計(jì)算方法

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

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

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

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

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

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

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

1 活化方程與CRAM方法

(1)

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

若定義矩陣元：

(2)

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

(3)

n(t)=eAtn(0)

(4)

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

(5)

(6)

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

2 活化計(jì)算中的關(guān)鍵問題研究

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

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

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

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

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

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

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

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

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

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

3 測試與驗(yàn)證

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

3.1 IAEA-ACB國際基準(zhǔn)例題

IAEA-ACB國際基準(zhǔn)例題旨在驗(yàn)證活化計(jì)算程序的四個方面：(1)正確讀取活化數(shù)據(jù)庫的能力；(2)在多步長的計(jì)算中，活化程序的計(jì)算結(jié)果與基準(zhǔn)值偏差不超過5%；(3)正確處理輕質(zhì)量核素(H、He等)；(4)正確處理數(shù)據(jù)庫中的同質(zhì)異能核素。該國際基準(zhǔn)例題包括Fe和Cr兩個子例題，覆蓋了常見的中子活化反應(yīng)和衰變反應(yīng)類型。根據(jù)初始條件，F(xiàn)e的質(zhì)量為1kg，Cr的核子密度為1.0E+20cm-3，輻照時間均為一年，中子通量譜參見文獻(xiàn)[22]。本文的活化計(jì)算結(jié)果與基準(zhǔn)值的相對誤差參見圖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)構(gòu)建方法得到的核素種類與基準(zhǔn)值列出的核素種類一致。對于Fe和Cr兩個子例題，經(jīng)大規(guī)模矩陣自適應(yīng)降階方法處理后的矩陣求解效率相對于處理前的矩陣求解效率分別提高了1.82倍和2.11倍。同時本文在同質(zhì)異能核素(52mMn、58mMn、44mSc、45mSc)、輕質(zhì)量核素(1H、2H、3H等)、中等質(zhì)量核素(Fe、Cr、Mn、Ni等)等活化計(jì)算中常見核素的計(jì)算結(jié)果與基準(zhǔn)值的相對誤差均低于0.5%，驗(yàn)證了本文發(fā)展的大規(guī)模矩陣自適應(yīng)降階方法的精確性。59Co和59Ni的計(jì)算結(jié)果與基準(zhǔn)值相對偏差較大，分別為0.75%和0.71%，但中子活化過程中這兩種核素的密度非常小，與初始核素密度相差1015～1018，對于實(shí)際工程中的影響非常小?？傮w來看，本文的計(jì)算結(jié)果與基準(zhǔn)值的相對偏差遠(yuǎn)低于基準(zhǔn)例題中要求的5%，滿足基準(zhǔn)例題中對活化計(jì)算程序的精度要求。

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

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

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

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

圖5 新型鋯合金N18活度計(jì)算結(jié)果Fig.5 The activity results of new zirconium alloy N18

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

4 結(jié)語

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

致謝

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

[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

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

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

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

TL329.2

A

0258-0918(2017)01-0065-08