朱金輝，陶應(yīng)龍，卓 俊，謝紅剛，左應(yīng)紅，商 鵬，韋 源，牛勝利

（西北核技術(shù)研究所，西安710024；強(qiáng)脈沖輻射環(huán)境模擬與效應(yīng)國(guó)家重點(diǎn)實(shí)驗(yàn)室，西安710024）

經(jīng)過(guò)半個(gè)多世紀(jì)的發(fā)展，特別是近年來(lái)隨著高性能計(jì)算技術(shù)的飛速提升，蒙特卡羅（MC）方法已成為解決粒子輸運(yùn)問(wèn)題的首選工具［1］。目前，比較成熟的粒子輸運(yùn)大型MC計(jì)算程序主要有MCNP、EGS、PENELOPE 和 GEANT4等［2－4］。這些程序功能強(qiáng)大、通用性好，被廣泛應(yīng)用于核技術(shù)相關(guān)領(lǐng)域的模擬計(jì)算中，為核工程技術(shù)及輻射環(huán)境的研究與發(fā)展發(fā)揮了重要的支撐作用［5－6］。但是針對(duì)一些特定模擬計(jì)算問(wèn)題，如粒子次級(jí)反應(yīng)過(guò)程中的特定參數(shù)獲取、物理過(guò)程之間的耦合計(jì)算等，這些程序皆無(wú)法滿足需求，且由于這些程序結(jié)構(gòu)復(fù)雜和龐大，在其基礎(chǔ)上進(jìn)行二次開(kāi)發(fā)實(shí)現(xiàn)特定功能比較困難。

為此，本文采用模塊化的程序設(shè)計(jì)方法，在前期開(kāi)發(fā)的中子－光子輸運(yùn)程序TOPAN［7］的基礎(chǔ)上，開(kāi)發(fā)了能夠模擬光子、質(zhì)子、電子和中子與物質(zhì)相互作用的主要過(guò)程，粒子在幾何中輸運(yùn)及輸運(yùn)結(jié)果統(tǒng)計(jì)記錄的模塊化程序，形成了模塊化程序包PHEN（the transportation program of photon，hadron，electron，and neutron）。該程序包中的各個(gè)程序模塊能夠較方便地被調(diào)用，可根據(jù)用戶需要開(kāi)發(fā)特定功能程序，也能夠較方便地實(shí)現(xiàn)與其他物理過(guò)程的耦合計(jì)算。

1 PHEN程序包基本功能介紹

PHEN程序包采用FORTRAN語(yǔ)言編制，基于MC方法模擬各個(gè)要素［8］，能夠模擬光子、質(zhì)子、電子和中子與物質(zhì)相互作用及其輸運(yùn)問(wèn)題，包含粒子輸運(yùn)所需的各基本功能模塊。與常見(jiàn)的MCNP及GEANT4程序相比，PHEN程序包具有3個(gè)特點(diǎn)：

1）清晰的樹(shù)結(jié)構(gòu)幾何建模功能

在幾何建模方面，PHEN程序包采用樹(shù)結(jié)構(gòu)描述幾何模型［4］。樹(shù)結(jié)構(gòu)幾何模型的基本思想是利用子體和母體的關(guān)系建立具有代（母子體）關(guān)系：幾何結(jié)構(gòu)的樹(shù)狀關(guān)系圖，在幾何相關(guān)計(jì)算中只需處理本代和具有直接子母體關(guān)系的部分，從而減小復(fù)雜幾何結(jié)構(gòu)的幾何算法計(jì)算量，理論上可從O（N）量級(jí)降低到O（logkN）量級(jí)，其中，N為幾何體總數(shù)，k為樹(shù)結(jié)構(gòu)總代數(shù)。

樹(shù)結(jié)構(gòu)幾何的基本思想如圖1所示。根據(jù)圖1左側(cè)的結(jié)構(gòu)可以建立圖1右側(cè)的樹(shù)狀結(jié)構(gòu)，其中，⑨是③和⑧的母體，⑧是⑦和④的母體……，當(dāng)粒子位置在體⑧中時(shí)，求粒子前向與結(jié)構(gòu)界面交點(diǎn)的程序，只需處理⑧的子體（⑦和④）及其與母體（⑨）的交界面，不需要遍歷模型內(nèi)所有的結(jié)構(gòu)交界面，理論上用于計(jì)算粒子輸運(yùn)前方界面交點(diǎn)的時(shí)間將減少約2／3。

圖1 樹(shù)結(jié)構(gòu)幾何示例Fig.1Sketch of tree structure

PHEN程序包的粒子輸運(yùn)幾何建模模塊能夠建立各種標(biāo)準(zhǔn)體模型（包含各種二次曲面圍成的體）及其各種組合模型。采用PHEN程序的樹(shù)結(jié)構(gòu)幾何建模方法，建模過(guò)程中從最小的子體開(kāi)始，逐級(jí)建模，幾何結(jié)構(gòu)清晰，可以較容易地建立較為復(fù)雜的幾何模型。圖2為利用PHEN幾何模塊建立的計(jì)算模型示意圖，圖中不同顏色代表不同的材料。

圖2 利用PHEN程序幾何模塊建立的計(jì)算模型Fig.2Acalculation model established by PHEN geographic module

PHEN程序包的幾何模塊，還可以用于射線追蹤模擬計(jì)算，比如可見(jiàn)光的輸運(yùn)模擬等［9］。

2）帶電粒子輸運(yùn)的功能拓展

在帶電粒子輸運(yùn)的MC模擬中，低能帶電粒子的輸運(yùn)一直是模擬的難點(diǎn)之一。在PHEN程序包的研發(fā)過(guò)程中，研究了低能電子直接MC模擬方法［10］，并將電子輸運(yùn)的直接模擬方法和壓縮歷史方法相結(jié)合，建立了適合于能量在50eV～1GeV范圍的電子輸運(yùn)MC方法，即混合模擬方法。把低能電子與物質(zhì)的相互作用分為能量損失較少的“軟碰撞”和能量損失較大的“硬碰撞”。其中，“軟碰撞”發(fā)生次數(shù)特別多，采用壓縮歷史方法模擬；“硬碰撞”發(fā)生次數(shù)較少，采用直接模擬方法模擬。

在低能質(zhì)子與物質(zhì)的相互作用模擬中，非彈性散射物理過(guò)程的模擬采用預(yù)復(fù)合模型［5］，能夠模擬能量為10keV的質(zhì)子在常用材料中的輸運(yùn)過(guò)程。另外，針對(duì)中子與含氫材料產(chǎn)生反沖質(zhì)子的物理現(xiàn)象，在中子彈性散射模塊中添加了產(chǎn)生次級(jí)質(zhì)子的模型，能夠模擬由彈性碰撞產(chǎn)生的反沖質(zhì)子，結(jié)合質(zhì)子輸運(yùn)模塊，能夠?qū)崿F(xiàn)中子輸運(yùn)－反沖質(zhì)子－質(zhì)子輸運(yùn)過(guò)程的一體化模擬。

3）可以定制的數(shù)據(jù)統(tǒng)計(jì)功能［3］

根據(jù)MC模擬的特點(diǎn)，并借鑒相關(guān)研究工作［11］的思路，采用了粒子標(biāo)識(shí)方法記錄粒子及其次級(jí)粒子的輸運(yùn)歷史，并根據(jù)標(biāo)識(shí)變量的值進(jìn)行統(tǒng)計(jì)，以達(dá)到求解特定物理量的目的。除了記錄射線種類、位置、方向、能量、權(quán)重等射線基本參數(shù)外，專門設(shè)置了粒子標(biāo)識(shí)數(shù)組用來(lái)標(biāo)識(shí)射線產(chǎn)生的位置、作用核素、反應(yīng)類型等感興趣參數(shù)。利用粒子標(biāo)識(shí)數(shù)組，可以對(duì)具體特定屬性的粒子進(jìn)行統(tǒng)計(jì)分析，得到某種材料或某種反應(yīng)產(chǎn)生的粒子數(shù)對(duì)總計(jì)數(shù)的貢獻(xiàn)，如特定材料（n，2n）反應(yīng)對(duì)中子通量的貢獻(xiàn)，光子康普頓散射對(duì)光子點(diǎn)通量的貢獻(xiàn)等。通過(guò)粒子標(biāo)識(shí)和用戶定義響應(yīng)函數(shù)，還可以擴(kuò)充數(shù)據(jù)統(tǒng)計(jì)模塊功能。

2 中子與典型核素產(chǎn)生的次級(jí)γ射線特征參數(shù)計(jì)算

獲得不同能量中子與特定核素相互作用過(guò)程中產(chǎn)生的次級(jí)γ射線特征參數(shù)，對(duì)于分析中子－γ聯(lián)合屏蔽問(wèn)題具有重要意義。在模擬過(guò)程中，希望直接記錄中子與原子核相互作用產(chǎn)生的次級(jí)γ信息，而通用的MC模擬程序，如MCNP，在粒子輸運(yùn)模擬的統(tǒng)計(jì)記錄過(guò)程中難以排除介質(zhì)在后續(xù)輸運(yùn)過(guò)程中對(duì)次級(jí)γ射線參數(shù)的影響。為此，在PHEN程序包的中子及γ輸運(yùn)過(guò)程模塊的基礎(chǔ)上進(jìn)行二次開(kāi)發(fā)，能夠計(jì)算不同能量的中子與典型核素相互作用后產(chǎn)生的次級(jí)γ的產(chǎn)額及平均能量。

計(jì)算中不模擬粒子的幾何輸運(yùn)過(guò)程，直接模擬特定能量中子與特定原子核的相互作用過(guò)程，統(tǒng)計(jì)產(chǎn)生的次級(jí)γ射線的產(chǎn)額及能量。設(shè)入射中子總數(shù)為Nn，所產(chǎn)生的次級(jí)γ的總數(shù)為Nγ，次級(jí)γ總能量為Eγ，tot。設(shè)中子能量為En，統(tǒng)計(jì)得到能量為En的中子產(chǎn)生的次級(jí)γ的平均能量Eγ，ave為

定義次級(jí)γ產(chǎn)額Y為

其中，Y表示中子與原子核相互作用過(guò)程產(chǎn)生次級(jí)γ的概率。

針對(duì)典型核素H、C、N、Fe及U，計(jì)算得到了不同能量中子作用下產(chǎn)生的次級(jí)γ產(chǎn)額及平均能量，如圖3所示。

圖3 不同能量中子在典型核素中產(chǎn)生的次級(jí)γ產(chǎn)額及平均能量Fig.3Yield and average energy of secondaryγ－rays generated by neutron interacting with typical nuclei

從圖3可以看到，不同核素與中子相互作用的次級(jí)γ產(chǎn)額和平均能量存在很大差異，通常高Z核素的次級(jí)γ產(chǎn)額很高，但平均能量相對(duì)較小。

3 基于指向概率法的點(diǎn)通量角度譜計(jì)算

在常見(jiàn)粒子輸運(yùn)模擬程序中，MCNP程序具備中子及光子點(diǎn)通量計(jì)數(shù)功能，采用指向概率方法記錄點(diǎn)通量，通過(guò)計(jì)算粒子與介質(zhì)原子核反應(yīng)時(shí)指向記錄點(diǎn)的虛次級(jí)粒子對(duì)記錄點(diǎn)的通量貢獻(xiàn)，得到總的點(diǎn)通量計(jì)數(shù)。在具體應(yīng)用中，用戶有時(shí)還希望得到點(diǎn)通量的角度分布信息，但MCNP程序目前不具備該功能。

利用模塊化程序包PHEN，基于MCNP程序中的粒子點(diǎn)通量計(jì)數(shù)方法，開(kāi)發(fā)得到了包含角度分布信息的點(diǎn)通量計(jì)數(shù)功能。

為了便于記錄分析，需預(yù)先設(shè)定參考方向。進(jìn)行點(diǎn)通量計(jì)算時(shí)，在記錄虛次級(jí)粒子的點(diǎn)通量貢獻(xiàn)的同時(shí)，計(jì)算虛次級(jí)粒子運(yùn)動(dòng)方向與參考方向的夾角，并對(duì)其進(jìn)行統(tǒng)計(jì)記錄。

為了驗(yàn)證開(kāi)發(fā)的點(diǎn)通量角度譜分析功能的正確性，利用MCNP程序球面上面通量角度譜的對(duì)稱性進(jìn)行對(duì)比計(jì)算，計(jì)算模型如圖4所示。模型中材料為鐵。首先利用MCNP程序記錄3個(gè)球面S1，S2，S3上的面通量及其角度譜；其次利用PHEN程序包記錄3個(gè)點(diǎn)P1，P2，P3的點(diǎn)通量及其角度譜，參考方向?yàn)閤軸正方向；最后利用PHEN程序包記錄S1，S2，S3上的面通量及其角度譜。球心能量為1MeV中子產(chǎn)生的點(diǎn)通量和面通量計(jì)算結(jié)果，如圖5所示。

圖4 計(jì)算模型Fig.4Calculation model

圖5 中子通量角度譜結(jié)果Fig.5Angle spectra of neutron fluence

根據(jù)對(duì)稱性分析，上述兩種方法計(jì)算的中子通量角度譜應(yīng)該一致。從圖5可以看出，PHEN程序包計(jì)算給出的面通量和點(diǎn)通量及其角度分布與MCNP程序給出的結(jié)果基本一致，驗(yàn)證了本文開(kāi)發(fā)的點(diǎn)通量角度譜分析功能的正確性。需要說(shuō)明的是，本文開(kāi)發(fā)的點(diǎn)通量角度譜分析功能也適用于非對(duì)稱情況下的角度譜分析，用戶只需根據(jù)研究需求，預(yù)先指定參考方向，即可得到所需的角度分布信息。

4 中子輸運(yùn)時(shí)間特性模擬計(jì)算

在超臨界系統(tǒng)分析中，除了系統(tǒng)裂變中子的有效增殖因子keff，系統(tǒng)的時(shí)間增殖特性也是研究人員關(guān)心的重要方面。MCNP5程序具有成熟的裂變驛站方法，可以方便地給出系統(tǒng)的keff，但難以直接給出時(shí)間增殖特性。

本文利用模塊化程序包PHEN，開(kāi)發(fā)了粒子輸運(yùn)的時(shí)間驛站功能，能夠計(jì)算給出系統(tǒng)的時(shí)間增殖特性。同時(shí)參照MCNP5程序的裂變驛站方法，也能夠計(jì)算超臨界系統(tǒng)裂變中子的有效增殖因子［12］。

針對(duì)圖6所示的超臨界系統(tǒng)模型，利用PHEN程序包和MCNP5程序計(jì)算得到的有效增殖因子，如表1所列?？梢?jiàn)，2種程序的計(jì)算結(jié)果基本一致。

圖6 超臨界系統(tǒng)計(jì)算模型Fig.6Calculation model of a supercritical system

采用時(shí)間驛站功能，對(duì)圖6所示模型計(jì)算了不同參數(shù)條件下，系統(tǒng)內(nèi)中子總數(shù)隨時(shí)間的變化關(guān)系，如圖7所示。

圖7 中子總數(shù)隨時(shí)間的變化Fig.7Sum number of neutron vs.time

從圖7可以看出，PHEN程序包能夠有效處理與時(shí)間相關(guān)的粒子輸運(yùn)問(wèn)題，可用于超臨界系統(tǒng)中子輸運(yùn)模擬，能夠給出系統(tǒng)裂變中子的有效增殖因子及時(shí)間增殖特性。

5 總結(jié)與展望

針對(duì)特定問(wèn)題需要，開(kāi)發(fā)了能夠模擬粒子輸運(yùn)過(guò)程的模塊化程序包PHEN，具備基本的中子、光子、電子和質(zhì)子耦合輸運(yùn)蒙特卡羅模擬能力。PHEN程序包中的模塊化程序能夠較方便地被調(diào)用，既可以根據(jù)用戶需要開(kāi)發(fā)特定功能程序，也能夠較方便地實(shí)現(xiàn)與其他物理過(guò)程的耦合計(jì)算。PHEN程序包的基本功能模塊可用于粒子輸運(yùn)模擬技術(shù)研究和應(yīng)用研究。

在后續(xù)工作中，將不斷改進(jìn)完善PHEN程序包，如增加光電子耦合輸運(yùn)中的厚靶韌致輻射模型、增加切倫科夫輻射模型、改進(jìn)截面數(shù)據(jù)處理算法等；同時(shí)，要不斷提高程序包的接口適應(yīng)性，為粒子輸運(yùn)相關(guān)數(shù)值模擬工作提供更加完備好用的粒子輸運(yùn)模塊化程序包。