基于MCNP 的反應(yīng)堆建模方法

王武，夏虹，李偉，巴少華

1.哈爾濱工程大學(xué) 核安全與先進(jìn)核能技術(shù)重點(diǎn)實(shí)驗(yàn)室，黑龍江 哈爾濱 150001

2.哈爾濱工程大學(xué) 核安全與仿真技術(shù)國(guó)防重點(diǎn)學(xué)科實(shí)驗(yàn)室，黑龍江 哈爾濱 150001

由于蒙特卡羅(Monte Carlo,MC)方法相比于確定性方法(如Boltzmann 輸運(yùn)方程的數(shù)值解)，能更準(zhǔn)確地表示幾何結(jié)構(gòu)和核數(shù)據(jù)，例如：確定性方法由于使用數(shù)值求解技術(shù)，因此需要合理地簡(jiǎn)化研究對(duì)象的幾何形狀，并且使用多群截面數(shù)據(jù)來(lái)近似連續(xù)能量截面數(shù)據(jù)，而MC 方法既可以處理簡(jiǎn)單的幾何和多群截面數(shù)據(jù)，也可以處理復(fù)雜的幾何形狀和連續(xù)的截面數(shù)據(jù)。因此對(duì)于許多研究對(duì)象而言，由于其內(nèi)部結(jié)構(gòu)十分復(fù)雜且材料不均勻，MC 是一種更好的模擬方法。

MCNP 作為核能領(lǐng)域的MC 計(jì)算程序，發(fā)展最為成熟、最具代表性。近年來(lái)，國(guó)內(nèi)外學(xué)者利用MCNP 在反應(yīng)堆有效增殖因子[1-3]、反應(yīng)性[4-5]、中子通量和功率分布[3,6-8]以及中子能譜[9-10]計(jì)算等方面做了諸多研究。雖然前人利用MCNP 模擬反應(yīng)堆臨界做了一系列工作，但是，在建模過(guò)程中，多數(shù)學(xué)者只關(guān)心反應(yīng)堆堆芯部件的幾何結(jié)構(gòu)、材料以及數(shù)據(jù)庫(kù)選取，沒(méi)有關(guān)注影響模型準(zhǔn)確性的其他因素，如不可分辨共振區(qū)概率表、自由氣體熱散射的溫度修正以及熱中子S(α,β) 模型等。因此，本文以某壓水反應(yīng)堆為研究對(duì)象，綜合分析了上述因素對(duì)MCNP 模擬的反應(yīng)堆有效增殖因子、堆芯能譜以及功率分布的影響，從而提出一種合理的MCNP 反應(yīng)堆建模方法。

1 反應(yīng)堆模型描述

1.1 堆芯幾何描述

本文以某壓水反應(yīng)堆堆芯為研究對(duì)象，選取0 個(gè)有效滿(mǎn)功率運(yùn)行天(0 effective full power operating days，0EFPD)燃耗時(shí)刻堆芯熱態(tài)滿(mǎn)功率運(yùn)行工況進(jìn)行建模。堆芯組件采用7×7 排布方式，組件按控制棒和可燃毒物棒的裝載方式不同可分為7 類(lèi)，分別用A、B、C、D、E、F、G 表示。每個(gè)組件中都留有16 根導(dǎo)向管，用于容納控制棒和可燃毒物棒，每根導(dǎo)向管柵格占據(jù)4 個(gè)燃料棒柵格的位置。A、B、C、D、E 為含有控制棒的組件，其中A、B、C 為中心棒組，各插入12 根控制棒和4 根可燃毒物棒，D、E 為外圍棒組，各插入16 根控制棒；F、G 為不含控制棒的組件，其中F 類(lèi)組件中插入16 根可燃毒物棒，G 類(lèi)組件中插入12 根可燃毒物棒。各類(lèi)組件在堆芯中呈中心對(duì)稱(chēng)布置如圖1 所示。堆芯內(nèi)部結(jié)構(gòu)均按照實(shí)際尺寸進(jìn)行精細(xì)化建模，燃料棒氣隙、可燃毒物棒氣隙和組件間水隙等微小結(jié)構(gòu)也考慮在內(nèi)。

圖1 組件在堆芯中的排布

使用MCNP 中的lattice(lat)卡、fill 卡以及universe(u)卡對(duì)燃料組件幾何進(jìn)行分級(jí)描述。以E 類(lèi)組件為例，組件結(jié)構(gòu)分三級(jí)描述，第一級(jí)為組件，第二級(jí)為填充組件的柵格，第三級(jí)為填充柵格的燃料棒、控制棒以及水。首先，描述燃料棒以及周?chē)乃⒎峙涫澜珞w編號(hào)u=2。然后，描述控制棒，由于一個(gè)控制棒柵格占據(jù)4 個(gè)燃料棒柵格的位置，而使用lat 卡和fill 卡進(jìn)行矩陣式填充時(shí)柵格的大小必須一致，所以需要將控制棒柵格劃分成4 個(gè)大小相同的柵格并分別分配世界體編號(hào)u=3、u=4、u=5 和u=6，每個(gè)柵格中都包含有1/4 控制棒、1/4導(dǎo)向管以及除控制棒和導(dǎo)向管以外的水。接著，定義參考柵格，并使用lat 卡和fill 卡將各個(gè)世界體填充的柵格按矩陣形式進(jìn)行重復(fù)，其中參考柵格作為矩陣的元素。為所有柵格分配世界體編號(hào)u=1，最后使用fill=1 將柵格填入組件中。圖2 為MCNP5 構(gòu)建的E 類(lèi)組件模型的徑向截面圖，其中不同的顏色代表不同的材料。其他E 類(lèi)組件只需要使用like…but…卡進(jìn)行復(fù)制和坐標(biāo)變換。除燃料組件以外其他堆芯結(jié)構(gòu)(如圍板、吊籃、壓力容器等)的幾何描述相對(duì)比較簡(jiǎn)單，因此不再詳述。圖3 為MCNP 構(gòu)建的堆芯整體模型的徑向截面圖。

圖2 E 類(lèi)組件徑向截面

圖3 反應(yīng)堆MCNP 徑向截面

1.2 堆芯核素份額計(jì)算

堆芯各結(jié)構(gòu)部件的材料為：燃料棒和可燃毒物棒包殼、導(dǎo)向管為Zr-4 合金；圍板、吊籃為OCr18Ni10Ti；壓力容器焊層為奧氏體不銹鋼；壓力容器、上下管座為508-Ⅲ鋼。除堆芯結(jié)構(gòu)材料外，準(zhǔn)確描述燃料棒、可燃毒物棒以及控制棒中核素含量對(duì)于臨界計(jì)算至關(guān)重要。由于控制棒由天然鉿(Hf-nat)組成，因此無(wú)需計(jì)算核素份額。對(duì)于燃料棒以及可燃毒物棒，設(shè)計(jì)資料只給出0EFPD燃耗時(shí)刻堆芯主要同位素235U、238U、10B 的含量分別為100 300、2 848 000 和151 g。因此，其他核素含量需要根據(jù)給出的核素含量進(jìn)行推算。

0EFPD 燃耗時(shí)刻燃料棒中所含的核素有235U、238U 和16O，可燃毒物棒中所含核素有10B、11B、12C和Zr-nat，忽略Zr-2 合金中的微量元素?，F(xiàn)根據(jù)設(shè)計(jì)參數(shù)計(jì)算燃料棒中以及可燃毒物棒中各同位素質(zhì)量份額。

堆芯中UO2總含量為

式中：na為堆芯組件數(shù)，nf為組件中燃料棒數(shù)，Vpin為每根燃料棒的體積，ρf為燃料密度設(shè)計(jì)值。

UO2燃料中235U 的富集度為

式中C5為U 中235U 的核子數(shù)份額。

UO2燃料中235U 的質(zhì)量份額為

由此可以計(jì)算出堆芯中235U 的總含量為

由于根據(jù)堆芯物理參數(shù)計(jì)算出的235U 的總含量與設(shè)計(jì)資料給出的值之間的相對(duì)誤差為9.78%，為了保證堆芯235U 總量與設(shè)計(jì)值相符，需要對(duì)燃料密度進(jìn)行如下修正：

式中：ρf′為修正后的燃料密度，M5為設(shè)計(jì)資料中給出的堆芯中235U 的總含量。

通過(guò)式(2)可以計(jì)算出修正后的燃料密度為9.361 1 g/cm3。

UO2燃料中238U 的質(zhì)量份額為

UO2燃料中16O 的質(zhì)量份額為

根據(jù)式(1)、式(3)和式(4)計(jì)算出的235U、238U、16O這3 種同位素的質(zhì)量份額分別為0.029 97、0.851 47和0.118 56。

對(duì)于可燃毒物，首先計(jì)算單位體積可燃毒物中B 的含量：

式中：np為堆芯內(nèi)可燃毒物棒總數(shù)；Vp為每根可燃毒物棒芯體的體積；MB為堆芯中B 的總量，MB=M10/ε10，其中，M10為堆芯中10B 的總量，ε10為B-nat中10B 的豐度。

然后便可以計(jì)算出單位體積可燃毒物中10B、11B、12C 的含量分別為

式中：CC和CB分別為MCNP5 手冊(cè)中給出的核級(jí)B4C 中的C 和B 的核子數(shù)密度，分別取0.217 4 和0.782 6；AC和AB分別為12C、B-nat 的原子質(zhì)量數(shù)，分別取12 和10.8。

最后，單位體積可燃毒物中Zr 的含量為

式中 ρp為可燃毒物芯體的密度。

由此便可以計(jì)算出可燃毒物中10B、11B、12C、Zr-nat 這4 種元素的質(zhì)量份額分別為0.000 900、0.003 957、0.001 509 和0.993 634。

將上述材料份額分別填寫(xiě)到MCNP 中對(duì)應(yīng)柵元的材料卡Mn 中完成材料描述。由于MCNP 使用的是NJOY 處理的ACE 格式截面庫(kù)，NJOY 處理的截面庫(kù)已經(jīng)包括在截面庫(kù)評(píng)價(jià)溫度下的彈性、俘獲、裂變和其他低閾值吸收截面(小于1 eV)的多普勒展寬，因此為了考慮多普勒效應(yīng)，本研究在選擇堆芯不同構(gòu)件中相關(guān)核素的截面庫(kù)時(shí)，選擇評(píng)價(jià)溫度與該構(gòu)件溫度最相近的截面庫(kù)，對(duì)于235U、238U 還要考慮截面庫(kù)是否包含緩發(fā)中子截面數(shù)據(jù)。由于臨界計(jì)算基本不考慮粒子的抽樣問(wèn)題，因此，設(shè)定堆芯及其周?chē)Y(jié)構(gòu)的柵元重要性IMP:N=1；設(shè)定安全殼以外區(qū)域中的柵元重要性IMP:N=0，某個(gè)柵元的IMP:N=0 表示中子一旦進(jìn)入該柵元域便被“殺死”，中子歷史終止。

1.3 不可分辨共振區(qū)概率表

在不可分辨的共振區(qū)內(nèi)(例如，ENDF/B-VI中，235U 為2.25～25 keV、238U 為10～149.03 keV、239Pu 為2.5～30 keV)，連續(xù)能量中子截面顯示為能量的平滑函數(shù)。這是因?yàn)楣舱穹宓木嚯x太近以至于無(wú)法分辨，而并不是共振峰的缺失。平滑變化的截面未能考慮能量自屏蔽效應(yīng)，但它對(duì)于譜峰在或接近不可分辨共振區(qū)系統(tǒng)可能是重要的。

基于分層抽樣技術(shù)，MCNP5 生成了不可分辨共振區(qū)截面的概率表。只要單次碰撞的平均能量損失比共振的平均寬度大得多，也就是說(shuō)，如果窄共振近似有效，那么從這些概率表中對(duì)隨機(jī)游動(dòng)截面抽樣就是一種有效的物理近似。

從概率表中抽樣截面很簡(jiǎn)單。在每一個(gè)入射能量處，都有一個(gè)累積概率表，以及與這些概率相對(duì)應(yīng)的近總、彈性、裂變和輻射俘獲截面值以及熱沉積數(shù)。這些數(shù)據(jù)補(bǔ)充了通常的連續(xù)數(shù)據(jù)。如果概率表關(guān)閉(在MCNP5 輸入文件數(shù)據(jù)卡中添加PHYS:N J J 1 J 卡)，則使用平滑截面；但是如果概率表打開(kāi)(MCNP5 默認(rèn))，并且如果碰撞能量在不可分辨共振區(qū)內(nèi)，則從概率表中抽樣截面。

1.4 自由氣體熱散射模型

中子和原子之間的碰撞受到原子熱運(yùn)動(dòng)的影響，在大多數(shù)情況下，碰撞也受到附近其他原子的存在的影響。MCNP 采用基于自由氣體近似的熱散射來(lái)解釋熱運(yùn)動(dòng)。由于柵元溫度會(huì)影響彈性散射截面，而自由氣體熱散射的默認(rèn)值為室溫，因此，只要柵元處于其他溫度，就必須在TMP 卡上給出這些值。

MCNP 程序在加載截面時(shí)，如果數(shù)據(jù)庫(kù)中的總截面和彈性截面的溫度與模型中TMP 卡指定的柵元溫度不同，則MCNP 將根據(jù)TMP 卡指定的柵元溫度對(duì)數(shù)據(jù)庫(kù)中的總截面和彈性截面進(jìn)行熱調(diào)整。如果不同的柵元有不同的溫度，MCNP 首先將截面調(diào)整到0K 下的截面，然后在運(yùn)輸過(guò)程中將其再次調(diào)整到柵元溫度對(duì)應(yīng)的截面。由于反應(yīng)堆堆芯溫度高于室溫，調(diào)整后的彈性散射截面會(huì)增大。本研究在TMP 卡上指定堆芯各柵元的平均溫度來(lái)模擬溫度對(duì)自由氣體熱散射模型的影響。

1.5 熱中子S(α,β)模型

S(α,β)模型考慮了化學(xué)(分子)結(jié)合和晶體效應(yīng)，當(dāng)中子的能量足夠低時(shí)(通常針對(duì)能量小于4 ev 的中子)，這些效應(yīng)就變得很重要。S(α,β)模型允許2 個(gè)過(guò)程：1)非彈性散射，截面 σin和相互耦合的能量、散射角均源于ENDFS(α,β)散射定律；2)彈性散射，截面σel和散射角由固體晶格參數(shù)導(dǎo)出。

對(duì)于材料卡Mn 中指定的核素，可以使用自由氣體模型直到低于某一能量，此時(shí)若S(α,β)數(shù)據(jù)存在，S(α,β)模型將自動(dòng)取代自由氣體模型。通常，自由氣體模型用于材料的每一種同位素，直到中子能量低于幾個(gè)電子伏特，然后對(duì)MTn 卡上指定物質(zhì)的同位素進(jìn)行S(α,β)處理。本研究將堆芯中冷卻劑輕水設(shè)置為S(α,β)模型，選擇輕水的S(α,β)表lwtr.03t，lwtr.03t 為500 K 時(shí)輕水的熱中子S(α,β)表，該表只提供輕水中1H 的熱中子非彈性散射數(shù)據(jù)，16O 仍然是默認(rèn)的自由氣體熱散射模型。

2 臨界計(jì)算

臨界計(jì)算最重要的是判斷裂變?cè)捶植寂cKeff是否收斂，因此MCNP5 對(duì)于臨界計(jì)算要求如下[11-12]：1)對(duì)所有可裂變物質(zhì)進(jìn)行充分抽樣；2)確保在有效循環(huán)前獲得基本特征值；3)最終的統(tǒng)計(jì)標(biāo)準(zhǔn)差應(yīng)該小于0.001。MCNP5 具有檢查功能以確保滿(mǎn)足前2 個(gè)要求，如果不滿(mǎn)足MCNP5 將在輸出文件中給出警告信息。

為了對(duì)可裂變物質(zhì)進(jìn)行充分地抽樣，并盡量減少模擬粒子數(shù)，本研究在KSRC 卡上給出每個(gè)燃料組件的中心燃料棒的中心點(diǎn)作為初始源位置。為確保在有效循環(huán)前獲得基本特征值，本研究在KCODE 卡中設(shè)置跳過(guò)的無(wú)效循環(huán)次數(shù)為100。計(jì)算結(jié)果顯示，相關(guān)設(shè)置均滿(mǎn)足要求。

對(duì)于統(tǒng)計(jì)標(biāo)準(zhǔn)差，其大小主要由計(jì)算量決定。用MCNP5 進(jìn)行臨界計(jì)算的計(jì)算量由KCODE卡中的總循環(huán)次數(shù)KCT 和每次循環(huán)的源粒子數(shù)NSRCK 決定。中心極限定理指出，在某置信水平下，統(tǒng)計(jì)標(biāo)準(zhǔn)差與模擬粒子數(shù)的平方成反比，因此欲使統(tǒng)計(jì)標(biāo)準(zhǔn)差減半，需將模擬粒子數(shù)增大為原來(lái)的4 倍，計(jì)算時(shí)間也會(huì)增大為原來(lái)的4 倍。本研究綜合考慮計(jì)算時(shí)間與計(jì)算精度，選擇KCT為500，NSRCK 為100 000，計(jì)算結(jié)果顯示標(biāo)準(zhǔn)差約為0.000 1，遠(yuǎn)小于0.001，可以認(rèn)為足夠精確。

在設(shè)置好臨界源后，分別模擬不考慮不可分辨共振區(qū)概率表、自由氣體熱散射的溫度修正、熱中子S(α,β)模型等因素以及考慮上述所有因素下的堆芯臨界，并使用基于網(wǎng)格計(jì)數(shù)的*FMESHn卡記錄堆芯徑向各組件功率，使用基于柵元計(jì)數(shù)的Fn 卡和En 卡記錄堆芯能譜。

本研究中假設(shè)考慮所有因素的模型為case1，關(guān)閉不可分辨共振區(qū)概率表的模型為case2，去掉TMP 卡(忽略柵元溫度對(duì)自由氣體熱散射模型的影響)的模型為case3，去掉MTn 卡的模型為case4。將case2、case3、case4 的計(jì)算結(jié)果分別與case1 的計(jì)算結(jié)果進(jìn)行比較，分析各因素的影響。

表1 給出了case1—case4 的Keff計(jì)算結(jié)果以及case2、case3、case4 的Keff計(jì)算結(jié)果分別與case1的Keff計(jì)算結(jié)果的相對(duì)誤差?？梢钥闯觯篶ase3 和case4 的相對(duì)誤差均大于MCNP5 給出的統(tǒng)計(jì)誤差，而case1 的相對(duì)誤差小于MCNP5 給出的統(tǒng)計(jì)誤差。這表明不可分辨共振區(qū)概率表對(duì)Keff的計(jì)算結(jié)果的影響可以忽略，自由氣體熱散射的溫度修正、熱中子S(α,β)模型對(duì)Keff的計(jì)算結(jié)果的影響不可忽略。其中，case4 的相對(duì)誤差最大，而且得出了非保守的Keff計(jì)算結(jié)果。所以對(duì)于反應(yīng)堆臨界安全分析必須考慮熱中子S(α,β)模型。

表1 Keff 計(jì)算結(jié)果

計(jì)算case1—case4 的堆芯歸一化能譜，并使用全譜絕對(duì)誤差積分以及各能區(qū)相對(duì)誤差的二范數(shù)說(shuō)明各因素對(duì)能譜的影響程度。case2、case3、case4 能譜與case1 能譜的絕對(duì)誤差的全譜積分分別為4.435 7×10-6、4.943 4×10-5和1.894 6×10-4，低能區(qū)(小于1 eV)、中能區(qū)(1 eV～0.1 MeV)和高能區(qū)(0.1 MeV～20 MeV)相對(duì)誤差的均方根值分別為(0.055 0，0.009 9，0.061 9)、(0.057 1，0.012 1，0.066 8)和(0.310 0，0.034 4，0.083 5)。case4的誤差明顯大于case2 和case3 的誤差，而且尤其在低能區(qū)誤差最為顯著，相對(duì)誤差的均方根值高達(dá)0.31，這也是其Keff計(jì)算結(jié)果非保守的原因。所以在計(jì)算堆芯能譜時(shí)必須考慮自由氣體熱散射的溫度修正以及熱中子S(α,β)模型，而其他2 個(gè)因素相對(duì)于熱中子S(α,β)模型可以忽略。

表2 給出了case1—case4 的1/4 堆芯各組件歸一化相對(duì)功率分布計(jì)算結(jié)果，以及case2、case3、case4 的各組件功率分布計(jì)算結(jié)果分別與case1 的各組件功率分布計(jì)算結(jié)果的相對(duì)誤差。計(jì)算得到case2、case3 和case4 堆芯功率分布相對(duì)誤差的二范數(shù)和正向無(wú)窮范數(shù)分別為(0.029 4，0.012 59)、(0.031 3，0.015 89)和(0.041 0，0.018 72)，統(tǒng)計(jì)誤差的二范數(shù)為0.004 5?？梢钥闯鼋^大部分相對(duì)誤差都大于甚至遠(yuǎn)大于統(tǒng)計(jì)誤差，并且相對(duì)誤差的二范數(shù)遠(yuǎn)大于統(tǒng)計(jì)誤差的二范數(shù)。這表明不可分辨共振區(qū)概率表、自由氣體熱散射的溫度修正、熱中子S(α,β)模型對(duì)堆芯功率分布均有著較大的影響，在計(jì)算堆芯功率時(shí)不可忽略。并且，case4 相對(duì)誤差的二范數(shù)和無(wú)窮范數(shù)最大，case3 次之，case2 最小，表明熱中子S(α,β)模型對(duì)堆芯功率分布影響最大，自由氣體熱散射的溫度修正次之，不可分辨共振區(qū)概率表對(duì)堆芯功率分布影響最小。

表2 1/4 堆芯各組件相對(duì)功率分布值

3 結(jié)論

本文利用MCNP5 對(duì)某壓水反應(yīng)堆堆芯進(jìn)行精細(xì)的幾何描述，并計(jì)算了0EFPD 燃耗時(shí)刻堆芯燃料棒及可燃毒物棒中主要核素的質(zhì)量份額，完成了材料描述。在此基礎(chǔ)上，分別計(jì)算不考慮不可分辨共振區(qū)概率表、自由氣體熱散射的溫度修正、熱中子S(α,β)模型等因素下的反應(yīng)堆有效增值因子、堆芯能譜及功率分布，并以考慮上述所有因素下的計(jì)算結(jié)果為基準(zhǔn)進(jìn)行比較，得出以下結(jié)論。

1)不可分辨共振區(qū)概率表對(duì)Keff值計(jì)算結(jié)果的影響可以忽略，自由氣體熱散射的溫度修正、熱中子S(α,β)模型對(duì)Keff值計(jì)算結(jié)果的影響不可忽略，其中熱中子S(α,β)模型對(duì)臨界值的影響最為顯著，若忽略熱中子S(α,β)模型將得到非保守的Keff值。

2)S(α,β)模型對(duì)能譜有著顯著影響，尤其是在低能區(qū)影響最大，其他2 個(gè)因素對(duì)能譜的影響相對(duì)于熱中子S(α,β)模型可以忽略。

3)3 個(gè)因素對(duì)于堆芯功率分布均有較大影響，并且熱中子S(α,β)模型對(duì)堆芯功率分布影響最大，自由氣體熱散射的溫度修正次之，不可分辨共振區(qū)概率表對(duì)堆芯功率分布影響最小。