球床高溫堆平衡態(tài)燃耗計算程序的開發(fā)

朱貴鳳，鄒 楊，李明海，嚴(yán) 睿，彭紅花，徐洪杰，＊

（1.中國科學(xué)院 上海應(yīng)用物理研究所，上海 201800；2.中國科學(xué)院 核輻射與核能技術(shù)重點(diǎn)實驗室，上海 201800；3.中國科學(xué)院大學(xué)，北京 100049）

球床高溫堆平衡態(tài)燃耗計算程序的開發(fā)

朱貴鳳1，2，3，鄒 楊1，2，李明海1，2，嚴(yán) 睿1，2，彭紅花1，2，徐洪杰1，2，＊

（1.中國科學(xué)院 上海應(yīng)用物理研究所，上海 201800；2.中國科學(xué)院 核輻射與核能技術(shù)重點(diǎn)實驗室，上海 201800；3.中國科學(xué)院大學(xué)，北京 100049）

基于MCNP5和ORIGEN2耦合方法，開發(fā)了平衡態(tài)下球床高溫堆的燃耗計算程序PBRE，用于堆的性能價值分析。為節(jié)省蒙特卡羅計算時間，對迭代收斂的方法進(jìn)行優(yōu)化，使之可在10個迭代步內(nèi)收斂。使用PBRE對清華大學(xué)HTR-10進(jìn)行建模計算，得到的平均卸料燃耗深度與文獻(xiàn)報道值一致，表明PBRE程序適用于球床堆平衡態(tài)的燃耗分析。

球床堆；平衡態(tài)；燃耗；PBRE

近年來，球床高溫堆因其具有較高的熱轉(zhuǎn)化效率及可實現(xiàn)高溫制氫的潛在產(chǎn)業(yè)優(yōu)勢，受主要核大國的推崇。球床高溫堆的一顯著特點(diǎn)是在線換料，能實現(xiàn)較高的燃耗深度，因而具有較好的燃料利用率，但這種流動性也決定了其燃耗分析與以往的壓水堆相比更困難。

早期的球床堆燃耗程序有VSOP［1］，其將時間項和空間項分開處理，各區(qū)間做完1次點(diǎn)燃耗計算后，再進(jìn)行1次空間平移處理，以此完成區(qū)間滯留時間內(nèi)的流動燃耗計算。由于對堆的性能價值評估只需分析堆平衡態(tài)時的燃耗深度等物理參數(shù)，隨后的確定論球床堆燃耗程序大部分是基于平衡態(tài)而開發(fā)的，如PREC［2］和PEBBED［3］等。PREC只局限于計算1次通過的換料方式；而PEBBED考慮再循環(huán)燃料在入口的邊界條件，使用雙迭代計算，可適用于各種換料方案的設(shè)計。

為避免確定論方法的幾何局限性和等效方法帶來的誤差，出現(xiàn)了基于蒙特卡羅方法開發(fā)的平衡態(tài)燃耗程序。2007年，F(xiàn)ratoni［4］在MOCUP程序的基礎(chǔ)上開發(fā)了球床堆平衡態(tài)燃耗計算程序，根據(jù)全反射邊界單球模型計算提供堆芯燃耗成分。但其關(guān)于滯留時間的修正方法以及程序準(zhǔn)確性的基準(zhǔn)驗算未見相關(guān)報道。

本文根據(jù)PEBBED程序的思路，采用MCNP5和ORIGEN2耦合的方法，開發(fā)1套平衡態(tài)球床高溫堆的燃耗分析程序PBRE。針對蒙特卡羅計算耗時的不足，重點(diǎn)對其收斂方法進(jìn)行分析，以期穩(wěn)定加速收斂，最后采用HTR-10堆芯進(jìn)行程序驗證。

1 程序設(shè)計思路

燃耗耦合程序需聯(lián)合解中子輸運(yùn)、功率和燃耗3個方程。中子輸運(yùn)方程采用MCNP5求解，給出單群截面和相對中子注量率分布；功率方程解出源強(qiáng)并計算出實際的中子注量率；燃耗方程通過中子注量率和單群截面解出堆芯燃耗成分分布，供中子輸運(yùn)計算使用，整個迭代過程設(shè)為內(nèi)迭代計算。為實現(xiàn)keff收斂到預(yù)期值，需不斷修改換料方案，這個迭代過程稱為外迭代。具體耦合流程示于圖1。圖1中，大方框外為輸入?yún)?shù)，其中，下標(biāo)i為區(qū)間號，k為i－1以下的區(qū)間號，j為燃料球種類號；A為矩陣；ˉσi，j為對能譜和空間平均的單群截面；ˉφi為能譜和空間平均的相對中子注量率；ˉΦi為對應(yīng)的實際中子注量率；Ti，j為區(qū)間i燃料球j的滯留時間；Ni，j（t）為新料Nj經(jīng)i－1個區(qū)間輻照后在i區(qū)間t時刻的成分；ˉNi為各類型燃料平均核素密度之和；kn為第n次循環(huán)計算得到的keff；dev和σ為用戶定義的收斂域。每次循環(huán)計算判斷keff是否收斂在σ以內(nèi)，不收斂則更新單群截面和單群中子注量率繼續(xù)計算；收斂則進(jìn)一步判斷是否與預(yù)期值k0相差在dev以內(nèi)，若在dev域值以外則修改添球速率繼續(xù)計算，否則結(jié)束計算。

1.1 平衡態(tài)燃耗方程處理

與一般的燃耗方程不同，球床堆堆芯的燃耗方程［5］為：

圖1 燃耗耦合流程圖Fig.1 Flow chart of coupled burnup code

其中：N（r，t）為堆芯r處t時刻的燃料核素成分組成的向量；v（r，t）為r處t時刻燃料的流動速度矢量，其值由換料方案決定；A（φ（r，t），σ（r，t），λ）為矩陣系數(shù)，矩陣元是r處t時刻的中子注量率、單群截面和衰變常量的函數(shù)。在平衡態(tài)時，忽略局部的隨機(jī)擾動，假定堆芯各物理參數(shù)分布穩(wěn)定，式（1）不含時間變量，可轉(zhuǎn)化為：

式（2）中成分、中子注量率和單群截面與空間相關(guān)，使用MCNP5求解，必須對堆芯劃分區(qū)間，中子注量率和單群截面為區(qū)間內(nèi)的平均值，此即對空間項做離散化處理。為維持動態(tài)平衡，某區(qū)間內(nèi)的核素隨時間的演化與核素在空間上的泄漏等價，故式（2）可轉(zhuǎn)化成：

式中：ni，j（t）為區(qū)間i內(nèi)滯留時刻t的第j種燃料球成分；ˉφi和ˉσi分別為區(qū)間i內(nèi)的平均中子注量率和平均單群截面。區(qū)間滯留時間可通過流動速度和區(qū)間體積獲得，此處假設(shè)所有球在該區(qū)間內(nèi)具有相同的流動軌跡。式（3）可由ORIGEN2求解。至此，球床堆平衡態(tài)堆芯燃耗方程的求解，轉(zhuǎn)化為對各分區(qū)各球在其滯留時間內(nèi)的點(diǎn)燃耗計算。堆芯各區(qū)間成分由流經(jīng)其內(nèi)的各燃料球平均成分求和得到。

1.2 區(qū)域成分平均處理

結(jié)合式（3），區(qū)域內(nèi)單個種類燃料球的平均成分為：

式中，Ti，j為區(qū)間i內(nèi)的j類型燃料球的滯留時間。式中的Ai難以處理，將指數(shù)泰勒展開取二階近似，結(jié)合式（3）簡化為：

式中，ni，j（Ti，j＋Δt）和ni，j（Ti，j）可通過ORIGEN2求出。考慮到ORIGEN2給出的成分保留5位有效數(shù)字，Δt的選取不宜過小。

區(qū)域總的平均成分是流經(jīng)該區(qū)的各類型燃料球平均成分之和。若區(qū)域內(nèi)燃耗分布很大，得到的平均成分與實際成分對keff的影響可能不同，因此，要求堆芯區(qū)間劃分得足夠多。

1.3 MCNP5建模處理

MCNP5模型中燃料球并不流動，模型幾何結(jié)構(gòu)不更改，迭代中，只對燃耗的成分進(jìn)行更新。焚燒的成分既可是燃料，也可是含硼的石墨。VSOP中某區(qū)域的成分是從上一區(qū)域轉(zhuǎn)移下來的，因此，流道上的分段必須是等體積的，所有的區(qū)間必須是等流動時間間隔。而本方法中各區(qū)成分均從新燃料球沿流動軌跡依次解燃耗方程得到，空間和時間已合并，因此，對區(qū)間的劃分無約束。為實現(xiàn)不同燃料類型、富集度、碳／重金屬比（燃料裝載量）的球的混合流動，需堆芯初始按一定比例混合建模，混合比例依賴于換料方案。

MCNP5建模可能帶來的誤差是，與實際的燃料間隔著石墨包殼不同，各區(qū)間的燃料可能在一個球內(nèi)直接接觸，產(chǎn)生額外的空間屏蔽效應(yīng)，尤其是不同燃料類型之間。

1.4 換料方案處理

換料方案包括循環(huán)次數(shù)、循環(huán)流道、添球速率及堆外滯留時間等。為求解單個球的燃耗，需將換料方案轉(zhuǎn)化成各新球從開始加入堆芯到最后卸出堆芯的整個流動軌跡。所以換料方案給定了區(qū)間燃耗計算的次序以及區(qū)間平均成分和平均燃耗深度的計算方法。文中換料方案的修改只針對添球速率，循環(huán)軌跡不變。PBRE可實現(xiàn)在同一流道內(nèi)對不同燃料類型、富集度、燃料裝載量的球分別計算，并能保存其幾何上的異質(zhì)性。

2 燃料球添球速率修改方法

本文提供線性關(guān)系和正相關(guān)兩種方法修改燃料球的添球速率，計算模型則采用HTR-10堆芯。

2.1 線性關(guān)系法

添球速率V（kgU／d）與堆芯總功率P（MW）和堆芯平均卸料燃耗深度U（MW·d／kgU）的關(guān)系為：

平衡態(tài)時，不同堆芯平均卸料燃耗深度與keff在臨界前后存在近似線性關(guān)系（圖2）。因此，假定V與keff存在如下關(guān)系：

V＝B／（－keff＋A）（7）

式中，A和B為正數(shù)。A和B的數(shù)值由前兩次迭代求出，過程如下。

圖2 keff與平均卸料燃耗深度的關(guān)系Fig.2 keffvs average discharged burnup

第1次迭代：初裝載計算認(rèn)為V為無窮大，此時有：

第2次迭代：添球速率為猜測值V2，此時有：

第3次迭代：添球速率采用式（7）計算，其中，keff為預(yù)期設(shè)定的值。實際上，卸料燃耗深度與keff并非嚴(yán)格的線性關(guān)系，為能實現(xiàn)局部線性關(guān)系，A和B需不斷更新，更新方法為：

式（10）取絕對值是為了避免統(tǒng)計擾動帶來負(fù)值。

第n次迭代：當(dāng)keff收斂在設(shè)定域值內(nèi)時，不修改添球速率，因為此時受MCNP5統(tǒng)計擾動的影響大；若接下來的4次循環(huán)仍收斂，則完成計算并結(jié)束任務(wù)，否則修改添球速率繼續(xù)計算。

由于MCNP5存在統(tǒng)計誤差，因此收斂域不能較統(tǒng)計誤差小。

圖3示出了迭代收斂的結(jié)果，兩條虛線之間為1次外迭代，虛線之間的數(shù)據(jù)點(diǎn)為該外迭代下的內(nèi)迭代點(diǎn)。為分析初始設(shè)定添球速率對迭代次數(shù)的影響，分別計算單個通道初始新球添球速率為100、28.37、20、15g（UO2）／d的收斂情況。從圖3可知：不同初始添球速率均可在2～3個外迭代內(nèi)實現(xiàn)收斂，初始添球速率越接近收斂值（28g（UO2）／d附近），收斂越快；初始添球速率越大，第1次內(nèi)迭代數(shù)越少，情況a和b分別用了4次，情況c用了6次，情況d用了10次；情況a、b、c、d下達(dá)到收斂時，MCNP5總的計算次數(shù)分別為15、11、13、20，收斂時間不一致。結(jié)果表明，線性關(guān)系法在收斂計算上相當(dāng)有效，但計算時間與初始添球速率有關(guān)，需根據(jù)預(yù)期燃耗深度推測初始添球速率。

圖3 線性關(guān)系法不同初始添球速率下的keff迭代收斂Fig.3 Convergence of keffwith linear relation method in different initial feedrates

為減少M(fèi)CNP5的計算次數(shù)，將線性關(guān)系法的內(nèi)外迭代合并，計算結(jié)果示于圖4。情況d的震蕩很明顯，收斂很慢，情況a開始時在收斂域徘徊，20個迭代后發(fā)散。原因是內(nèi)迭代不穩(wěn)定的情況下，新燃料球的添球速率與keff的關(guān)系式（式（7））完全破壞，修正不具可靠性。

2.2 正相關(guān)法

內(nèi)外迭代合并的情況下，為保障迭代的收斂性，必須維持keff與添球速率的正相關(guān)性。給出一簡單的關(guān)系式：

式（11）雖收斂穩(wěn)定，但收斂速度很慢。為加速收斂，必須考慮添球速率變化幅度與keff變化幅度的關(guān)系。對式（7）兩邊微分，有：

圖4 內(nèi)外迭代合并下keff和添球速率的收斂Fig.4 Convergence of keffand feedrate in merged iteration process

A等于keff1，選取修正幅度為預(yù)期值keff0附近的幅度，則添球速率可修改為：

式（13）確保了V的數(shù)值為正，同時也考慮了初始裝載不同的影響。

采用正相關(guān)法計算得到的結(jié)果示于圖5，圖5中增加了初始添球速率1 000g（UO2）／d的情況（情況e）。keff在10個迭代步內(nèi)基本收斂，相比線性關(guān)系法，減少了一定的計算時間，同時對不同的添球速率不敏感，可靠性較高。

圖5 正相關(guān)法下的keff迭代收斂Fig.5 Convergence of keffwith positive correction method

3 HTR-10驗算

3.1 模型及區(qū)域劃分

采用清華大學(xué)HTR-10作為PBRE程序的驗算基準(zhǔn)題。堆芯模型示于圖6，具體模型參數(shù)見文獻(xiàn)［6］。每個新燃料球含5g鈾，235U富集度為18%，活性區(qū)球填充因子為60%，所有材料均計入硼當(dāng)量。受MCNP5允許的最大存儲空間限制，本文不對堆芯進(jìn)行溫度分區(qū)，根據(jù)文獻(xiàn)［7］提供的溫度分布圖，將溫度進(jìn)行平均近似處理，選取堆芯燃料溫度為880K，堆芯石墨溫度設(shè)定為860K，周圍反射層石墨取中間位置的溫度，設(shè)為800K。堆芯頂空腔高41.3cm，平衡態(tài)時，堆芯裝載2.7萬個燃料球，計算中卸料管內(nèi)無燃料球。堆芯活性區(qū)徑向根據(jù)文獻(xiàn)［8］提供的7∶10∶12∶16∶21體積比例劃分，但劃分曲線與實際流道曲線不同。軸向等高度將每個流道分成7個區(qū)，堆芯共劃分成35個區(qū)。燃料球的流動方式采用5次通過堆芯的方案［8－10］。

圖6 HTR-10建?？v、橫剖面Fig.6 Vertical and horizontal cross sections of HTR-10modeling

3.2 計算結(jié)果分析

控制棒全抽出，平衡態(tài)時keff為1.01左右［8］。設(shè)定每個流道的新燃料球添球速率為28.37g（UO2）／d，此時卸料平均燃耗深度為80MW·d／kgU，計算得到keff收斂為1.007。設(shè)定預(yù)期keff為1.01，收斂域為0.001，則計算得到每流道新燃料球添球速率為28.8g（UO2）／d，此時平均卸料燃耗深度為78.9MW·d／kgU。

表1列出了燃耗深度的分布（徑向區(qū)域標(biāo)號從小到大表示從內(nèi)到外，軸向區(qū)域標(biāo)號從小到大表示從上到下），從徑向上看，外道燃耗深度較內(nèi)道的大，主要因為外道滯留時間長；軸向上因為采用5次通過的方式，所以燃耗深度分布差別不大；堆芯整體平均燃耗深度為43.16MW· d／kgU。因為軸向劃分與VSOP并不相同，堆芯燃耗分布無法與文獻(xiàn)結(jié)果比較，但各流道的卸料燃耗深度與文獻(xiàn)［9］的值僅相差2MW· d／kgU，分布較為一致。

表1 堆芯各區(qū)燃耗深度分布Table 1 Burnup distribution of core

表2所列為各區(qū)域中子注量率分布，最高注量率位于堆芯中央，約1.0×1014cm－2· s－1。中子注量率徑向分布較平坦；在軸向上則因為泄漏的原因，中子注量率差異略大。

表2 堆芯各區(qū)中子注量率分布Table 2 Fluence rate distribution of core

表3列出了堆芯功率密度分布，最大功率密度為2.46MW／m3，堆芯平均功率密度為2MW／m3，功率峰因子很小，僅1.23，較文獻(xiàn)［9］的小8%，可能是堆芯分區(qū)少的緣故。與中子注量率和燃耗深度在徑向分布不同的是，功率密度在徑向上并非單調(diào)減少，而是在外流道略微增大，尤其是堆芯錐體卸出口的規(guī)律與之完全相反。因為HTR-10堆芯小，中子泄漏大，外流道受石墨反射層的影響大，能譜更軟，所以其裂變截面相對增大，導(dǎo)致功率密度增大。235U的單群裂變截面分布列于表4，整個分布呈堆芯中間小、外圍大的規(guī)律，與中子注量率相反。

表3 堆芯各區(qū)功率密度分布Table 3 Power density distribution of core

表4 堆芯各區(qū)235U單群裂變截面分布Table 4 One-group fission cross section distribution of235U in core

4 小結(jié)

本文針對球床高溫堆流動特性導(dǎo)致的燃耗分析難點(diǎn)，采用MCNP5和ORIGEN2耦合的方法，開發(fā)了1套平衡態(tài)球床高溫堆的燃耗分析程序PBRE，并采用該程序?qū)TR-10堆芯進(jìn)行了模擬，計算了燃耗深度、中子注量率、功率分布等數(shù)據(jù)。分析表明，keff＝1.01時，堆芯平均卸料燃耗以及各流道卸料燃耗深度分布與文獻(xiàn)值基本吻合。同時，經(jīng)實測，通過對迭代收斂的加速分析，有效地減少了計算時間，MCNP5需約10次計算即可得到平衡態(tài)參數(shù)。因此，基本驗證了PBRE程序的可靠性，可將該程序用于球床堆平衡態(tài)堆芯的相關(guān)分析計算。堆芯溫度分布的獲取、燃耗計算的溫度分區(qū)以及結(jié)合流動規(guī)律程序更好地給出換料方案，則將是PBRE下一步的研究內(nèi)容。

［1］ RüTTEN H J，HAAS K A，BROCKMAN H，et al.VSOP（99／05）computer code system［M］.Zentralbibliothek：Forschungszentrum Jülich GmbH，2005.

［2］ SEKIMOTO H，OBARA T，YUKINORI S，et al.New method to analyze equilibrium cycle of pebble-bed reactors［J］.Journal of Nuclear Science and Technology，1987，24（10）：765-772.

［3］ GOUGAR H D.Advanced core design and fuel management for pebble-bed reactors［D］.US：The Pennsylvania State University，2004.

［4］ FRATONI M.Development and applications of methodologies for the neutronic design of the Pebble Bed Advanced High Temperature Reactor（PB-AHTR）［D］.Berkeley，US：University of California，2008.

［5］ MASSIMO L.Physics of high-temperature reactors［M］.Holland：Elsevier，1976.

［6］ 經(jīng)滎清，楊永偉，許云林.蒙特卡羅方法用于HTR-10首次臨界燃料裝料預(yù)估的校算［J］.核動力工程，2005，26（1）：28-34.

JING Xingqing，YANG Yongwei，XU Yunlin.Application of Monte Carlo method for verification calculation in fuel loading prediction for first criticality of HTR-10［J］.Nuclear Power Engineering，2005，26（1）：28-34（in Chinese）.

［7］ GAO Z，SHI L.Thermal hydraulic calculation of the HTR-10for the initial and equilibrium core［J］.Nuclear Engineering and Design，2002，218（1）：51-64.

［8］ 吳宗鑫，經(jīng)滎清.HRT-10的燃料管理［J］.清華大學(xué)學(xué)報：自然科學(xué)版，2001，41（4）：120-123.

WU Zongxin，JING Xingqing.Fuel management of HTR-10［J］.J Tsinghua Univ：Sci ＆Tech，2001，41（4）：120-123（in Chinese）.

［9］ 經(jīng)滎清，張旭，羅經(jīng)宇.HTR-10MW高溫氣冷實驗堆換料方式的研究［J］.核科學(xué)與工程，1993，13（2）：119-125.

JING Xingqing，ZHANG Xu，LUO Jingyu.The study on the shuffling scheme in HTR-10MW test module［J］.Chinese Journal of Nuclear Science and Engineering，1993，13（2）：119-125（in Chinese）.

［10］楊永偉，經(jīng)滎清，許云林，等.HTR-10平衡態(tài)運(yùn)行方式研究［J］.清華大學(xué)學(xué)報：自然科學(xué)版，1995，35（6）：12-16.

YANG Yongwei，JING Xingqing，XU Yunlin，et al.Study of operation model at equilibrium state of core for HTR-10［J］.J Tsinghua Univ：Sci ＆Tech，1995，35（6）：12-16（in Chinese）.

Development of Burnup Calculation Code for Pebble-bed High Temperature Reactor at Equilibrium State

ZHU Gui-feng1，2，3，ZOU Yang1，2，LI Ming-hai1，2，YAN Rui1，2，
PENG Hong-hua1，2，XU Hong－jie1，2，＊
（1.Shanghai Institute of Applied Physics，Chinese Academy of Sciences，Shanghai 201800，China；2.Key Laboratory of Nuclear Radiation and Nuclear Energy Technology，Chinese Academy of Sciences，Shanghai 201800，China；3.University of Chinese Academy of Sciences，Beijing100049，China）

The burnup calculation code PBRE coupling MCNP5and ORIGEN2was developed for pebble-bed high temperature reactor at equilibrium state，and it can be used to analyze the neutronic performance of equilibrium core.The iteration method was optimized in order to save Monte Carlo calculation time，and the convergence can be reached in 10iterative steps.The average discharged burnup for HTR-10is consistent with literature，and it indicates that the PBRE is suitable to analyze the burnup for pebble-bed reactor at equilibrium state.

PBR；equilibrium state；burnup；PBRE

TL32

：A

：1000-6931（2015）05-0890-07

10.7538／yzk.2015.49.05.0890

2014-01-22；

2014-04-22

中國科學(xué)院戰(zhàn)略性先導(dǎo)科技專項資助項目（XDA02010200）；上海市科技創(chuàng)新重點(diǎn)項目資助（11JC1414900）；國家自然科學(xué)基金青年基金資助項目（11005148）；973計劃資助項目（2010CB934501）

朱貴鳳（1987—），男，江西九江人，博士研究生，核科學(xué)與工程專業(yè)

＊通信作者：徐洪杰，E-mail：xuhongjie＠sinap.ac.cn