基于蒙特卡羅方法的IVR熔融池內(nèi)熱源時(shí)序模型構(gòu)建及分析

陳俊逸，黃善仿，郝以昇，劉國(guó)棟，胡鈺文，黎 陽(yáng),*，宮厚軍，昝元鋒，郭嘯宇，駱 浩

(1.清華大學(xué) 工程物理系，北京 100084；2.中國(guó)核動(dòng)力研究設(shè)計(jì)院 中核核反應(yīng)堆熱工水力技術(shù)重點(diǎn)實(shí)驗(yàn)室，四川 成都 610041)

堆內(nèi)熔融物滯留(IVR)作為反應(yīng)堆嚴(yán)重事故的關(guān)鍵緩解策略，其具體內(nèi)容為：嚴(yán)重事故發(fā)生時(shí)，堆芯熔融物向壓力容器下封頭內(nèi)遷移形成熔融池。熔融池內(nèi)殘余核燃料及裂變產(chǎn)物釋熱對(duì)壁面持續(xù)加熱。冷卻水完全淹沒(méi)下封頭外表面，通過(guò)氣液兩相自然循環(huán)保證充分帶走壓力容器外表面熱量。

IVR由加州大學(xué)圣塔芭芭拉分校(UCSB)的Theofanous教授提出，并基于其自主開發(fā)的面向事故的風(fēng)險(xiǎn)分析方法(ROAAM)對(duì)反應(yīng)堆嚴(yán)重事故下形成的熔融池進(jìn)行了分析。結(jié)果表明，壁面熱流密度是否超過(guò)臨界熱流密度(CHF)，是評(píng)價(jià)IVR結(jié)構(gòu)性能的關(guān)鍵指標(biāo)[1]。

數(shù)年來(lái)國(guó)內(nèi)外開展了大量試驗(yàn)，以基于CHF探索嚴(yán)重事故下IVR的安全特性。Theofanous教授開展了COPO及ACOPO實(shí)驗(yàn)[2]，對(duì)下封頭二維熔融池模型進(jìn)行了傳熱分析，并得到了用于預(yù)測(cè)下封頭CHF隨角度變化的表達(dá)式。Sehgal教授開展了SEMICO實(shí)驗(yàn)[3]，對(duì)池內(nèi)雷諾數(shù)及壁面熱流密度進(jìn)行了探索。西安交通大學(xué)基于不同熔融池結(jié)構(gòu)開展了COPRA實(shí)驗(yàn)[4-5]，針對(duì)熔融池內(nèi)對(duì)流進(jìn)行CFD模擬，探索了速度場(chǎng)與溫度場(chǎng)的分布。俄羅斯Korchatov研究所于RASPLAV與MASCA實(shí)驗(yàn)中發(fā)現(xiàn)了重金屬層分布于底部的熔融池三層構(gòu)型[6]，進(jìn)而分析了IVR的安全性。

上述實(shí)驗(yàn)與模擬，均為人工設(shè)置恒定熱源模擬體積熱源，不能精確反映熔融池釋熱隨時(shí)間與空間的變化，增加了熱流密度測(cè)量及模擬的不確定性。為精準(zhǔn)模擬熔融池內(nèi)衰變熱，清華大學(xué)基于反應(yīng)堆蒙特卡羅(RMC)程序及內(nèi)耦合燃耗/衰變熱計(jì)算程序DEPTH提出了面向熔融池內(nèi)物理過(guò)程的新衰變熱計(jì)算方法，并驗(yàn)證了其可行性與準(zhǔn)確性[7]。然而，堆芯運(yùn)行不同時(shí)刻的核素組分不同，熔融池釋熱規(guī)律及物理過(guò)程不盡相同。前述研究未能進(jìn)一步探索堆芯運(yùn)行不同時(shí)間形成熔融池的釋熱特點(diǎn)及物理規(guī)律，因此需要建立熔融池內(nèi)熱源時(shí)序模型，并分析相應(yīng)的物理機(jī)理。

本文基于RMC程序和BEAVRS基準(zhǔn)題，通過(guò)RMC程序的構(gòu)建實(shí)體幾何(CSG)方法對(duì)熔融池結(jié)構(gòu)進(jìn)行精細(xì)化建模與燃耗區(qū)劃分。通過(guò)蒙特卡羅方法模擬熔融池內(nèi)中子輸運(yùn)及相關(guān)核反應(yīng)，使用線性子鏈方法求解燃耗方程進(jìn)行IVR熔融池內(nèi)燃耗計(jì)算及裂變產(chǎn)物衰變鏈模擬，分析熔融池內(nèi)物理過(guò)程，以建立精準(zhǔn)的熔融池內(nèi)熱源時(shí)序模型，并從物理分析角度對(duì)變化趨勢(shì)進(jìn)行評(píng)價(jià)，為IVR策略安全及設(shè)計(jì)提供參照。

1 計(jì)算方法

1.1 反應(yīng)堆蒙特卡羅程序

RMC程序是清華大學(xué)工程物理系REAL實(shí)驗(yàn)室自主研發(fā)的蒙特卡羅中子輸運(yùn)程序，兼具臨界計(jì)算、燃耗計(jì)算、中子/光子/電子耦合計(jì)算、全堆換料計(jì)算、衰變熱計(jì)算等功能，具有良好的并行計(jì)算效率。通過(guò)面向BEAVRS等基準(zhǔn)題的驗(yàn)證計(jì)算，證明RMC程序具有較高的準(zhǔn)確性[8]。

DEPTH模塊為RMC程序內(nèi)置燃耗計(jì)算模塊，實(shí)現(xiàn)了內(nèi)耦合的燃耗計(jì)算及衰變熱計(jì)算等功能，具備大規(guī)模的燃耗計(jì)算能力，可用于計(jì)算壓水堆全堆燃耗[9-10]及模擬熔融池中核素組分變化。燃耗計(jì)算與衰變熱計(jì)算的準(zhǔn)確性已在文獻(xiàn)[7,10]中得到了驗(yàn)證。

1.2 計(jì)算流程

針對(duì)RMC程序及DEPTH程序，本研究結(jié)合全堆計(jì)算與燃耗/衰變熱計(jì)算兩個(gè)模塊，進(jìn)行迭代求解，計(jì)算流程圖如圖1所示。具體計(jì)算步驟如下：1)BEAVRS堆芯建模及全堆輸運(yùn)計(jì)算，在指定燃耗步輸出堆芯中核素的絕對(duì)原子密度；2)基于BEAVRS基準(zhǔn)題對(duì)堆芯組分質(zhì)量進(jìn)行計(jì)算，得到三層熔融池各層體積，構(gòu)建熔融池模型；3)細(xì)網(wǎng)格劃分熔融池燃耗區(qū)，熔融池燃耗計(jì)算得到各燃耗區(qū)的中子通量及核反應(yīng)功率，更新熔融池材料；4)對(duì)熔融池進(jìn)行衰變熱計(jì)算[7]，輸出燃耗區(qū)內(nèi)釋熱功率，更新材料并進(jìn)入下一個(gè)全堆燃耗步長(zhǎng)；5)循環(huán)運(yùn)行步驟1～4，直到全堆計(jì)算進(jìn)入最終燃耗步，結(jié)束計(jì)算。

圖1 熔融池分時(shí)刻釋熱功率計(jì)算流程圖

1.3 衰變熱計(jì)算

熔融池釋熱過(guò)程中，不能忽視衰變熱的貢獻(xiàn)。RMC程序結(jié)合了衰變熱計(jì)算與燃耗計(jì)算以描述核素的變化規(guī)律。燃耗計(jì)算的核心是求解經(jīng)典點(diǎn)燃耗方程，描述核反應(yīng)系統(tǒng)中核素隨時(shí)間的變化規(guī)律。燃耗鏈中核素與時(shí)間相關(guān)的點(diǎn)燃耗方程如下：

(1)

(2)

(3)

式中：λi為第i個(gè)核素的衰變常量；φ為中子通量；σi,j為第i個(gè)核素反應(yīng)生成第j個(gè)核素的反應(yīng)截面。

常用的燃耗方程求解方法主要有兩類，分別為線性子鏈法與矩陣指數(shù)法。本文采用線性子鏈法(TTA)將復(fù)雜的燃耗鏈分解為獨(dú)立的線性燃耗子鏈，然后對(duì)每條線性子鏈解析求解。

對(duì)1條線性子鏈，假設(shè)頂端核素初始核密度為n1(0)，則線性子鏈第k個(gè)核素經(jīng)時(shí)間t后的核素密度nk(t)為：

(4)

(5)

(6)

針對(duì)反應(yīng)堆停堆以及熔融池等場(chǎng)景中常見的衰變熱計(jì)算需求，DEPTH模塊內(nèi)置了衰變熱計(jì)算功能[5]。基于中子注量與絕對(duì)原子密度直接計(jì)算燃耗區(qū)內(nèi)的衰變熱釋熱功率Q，即：

Q=Rq=Nλq

(7)

式中：R為核反應(yīng)率；q為核素每次衰變放出的熱量。

2 計(jì)算模型構(gòu)建

2.1 BEAVRS全堆模型構(gòu)建

RMC程序具有幾何模型構(gòu)建功能，能基于層級(jí)空間的幾何描述系統(tǒng)以及CSG方法進(jìn)行三維模型的構(gòu)建。

基于層級(jí)空間的幾何描述系統(tǒng)通過(guò)空間結(jié)構(gòu)的層級(jí)嵌套，不僅能分層表示反應(yīng)堆各結(jié)構(gòu)，如燃料棒、導(dǎo)向管、儀表管、燃料組件、全堆堆芯等，還可用于下封頭熔融池的建模。RMC程序構(gòu)建的BEAVRS全堆幾何模型示意圖如圖2[8]所示。全堆計(jì)算設(shè)置為：每輪迭代采用100萬(wàn)粒子數(shù)，100個(gè)非活躍代和500個(gè)活躍代，收斂誤差為0.000 01，進(jìn)行全堆燃耗計(jì)算，獲取各燃耗步的核素信息作為熔融池建模邊界條件。

圖2 BEAVRS全堆幾何模型示意圖[8]

2.2 熔融池模型構(gòu)建

國(guó)內(nèi)外對(duì)熔融池幾何結(jié)構(gòu)的建模已有較多研究，對(duì)于三層熔融池，較有代表性的有Esmaili等[11]基于MASCA理論的三層模型、Seiler等[12]的保守三層模型以及Salay等[13]的基于熔融池?zé)峄瘜W(xué)性質(zhì)的三層模型，后者已應(yīng)用于MAAP5。

本文基于全堆熔化工況，作出下述假設(shè)以進(jìn)行熔融池各層參數(shù)計(jì)算：1)輕金屬層的成分為金屬Zr以及不銹鋼(SS)；2)氧化物層的主要成分為UO2及ZrO2，設(shè)置金屬Zr的氧化率為0.5；3)重金屬層的主要成分為金屬U和金屬Zr；4)裂變產(chǎn)物基于全堆燃耗計(jì)算結(jié)果，均勻分布在氧化物層與重金屬層中。

計(jì)算得到了三層熔融池結(jié)構(gòu)各層建模參數(shù)。本工作建模方法與文獻(xiàn)模型參數(shù)計(jì)算結(jié)果列于表1。

表1 本工作與文獻(xiàn)建模參數(shù)計(jì)算結(jié)果

2.3 熔融池燃耗區(qū)劃分

熔融池形成后，池內(nèi)可能存在活化裂變產(chǎn)物等核素放出中子并引發(fā)次級(jí)核反應(yīng)，對(duì)熔融池中不同區(qū)域的材料密度造成影響，進(jìn)而影響釋熱功率。需要對(duì)熔融池進(jìn)行細(xì)網(wǎng)格區(qū)域劃分，從而獲得更準(zhǔn)確的釋熱功率分布。

本文將下封頭熔融池沿z軸(高度)方向進(jìn)行分層，其中，裂變產(chǎn)物主要聚集于氧化物層。將輕金屬層分為5層、氧化物層分為12層、重金屬層分為4層。每層沿半徑方向分為11層，沿圓周方向分為18層。輕金屬層分989個(gè)燃耗區(qū)域，氧化物層分2 015個(gè)燃耗區(qū)域，重金屬層分287個(gè)燃耗區(qū)域，共3 291個(gè)燃耗區(qū)實(shí)現(xiàn)對(duì)下封頭熔融池三層結(jié)構(gòu)的完全分區(qū)和填充[5]。

嚴(yán)重事故下，壓力容器被兩相流體淹沒(méi)并發(fā)生沸騰換熱。為探索兩相流邊界條件中產(chǎn)生的氣泡可能對(duì)中子輸運(yùn)造成的影響，本文采用弦長(zhǎng)抽樣方法，將氣泡視為外部冷卻水基體中的彌散介質(zhì)，對(duì)下封頭外側(cè)兩相流體邊界條件進(jìn)行建模[14]。三層熔融池的下封頭及外部氣液兩相流模型如圖3所示。

圖3 熔融池燃耗區(qū)劃分策略

3 計(jì)算結(jié)果與分析

3.1 反應(yīng)堆運(yùn)行時(shí)刻選取

反應(yīng)堆運(yùn)行期間，中子變化復(fù)雜，核素分布多樣，因而在運(yùn)行任一時(shí)刻發(fā)生嚴(yán)重事故，形成的熔融池皆具有不同的釋熱特性。本文選取BEAVRS百萬(wàn)千瓦級(jí)壓水堆滿功率運(yùn)行中的數(shù)個(gè)時(shí)間點(diǎn)作為嚴(yán)重事故熔融池示例，分別為裝料后運(yùn)行5、30、40、70、150、200、250、300 d。通過(guò)DEPTH程序進(jìn)行BEAVRS全堆燃耗計(jì)算?；谙鄳?yīng)運(yùn)行時(shí)刻的核素信息，選取重要核素作為熔融池初態(tài)核素組成邊界條件，計(jì)算熔融池釋熱功率[15]。

3.2 熔融池時(shí)序模型構(gòu)建及分析

對(duì)熔融池體積熱源的模擬，一方面要考慮由于核素衰變釋放的衰變熱，另一方面要考慮熔融池中裂變碎片及自發(fā)裂變核素等釋放中子、殘余核燃料裂變等的核反應(yīng)釋熱。

通過(guò)RMC程序在熔融池內(nèi)初始化點(diǎn)中子源，能有效模擬熔融池中的中子輸運(yùn)。RMC程序設(shè)置每輪迭代采用10萬(wàn)粒子數(shù)，20個(gè)非活躍代和100個(gè)活躍代，收斂誤差為0.000 2，時(shí)間步長(zhǎng)為0.05 d。

反應(yīng)堆運(yùn)行5 d和40 d后熔融池釋熱功率變化規(guī)律示于圖4。由圖4a可看出，運(yùn)行5 d后，熔融池釋熱功率呈指數(shù)形式衰減，變化趨勢(shì)與文獻(xiàn)[7]計(jì)算結(jié)果相符。釋熱功率最大值出現(xiàn)在最初燃耗步，為2.5 MW，然后迅速下降，在第3燃耗步處下降到約50 kW且逐漸趨于平緩(圖4b)，于第5 d衰變熱功率降于初始的0.05%，可忽略。這體現(xiàn)了熔融池形成初期，熔融池內(nèi)中子對(duì)熔融池釋熱功率存在影響。釋熱功率較低是由于可裂變核素尚未大量裂變并產(chǎn)生裂變產(chǎn)物，抑制了衰變熱的產(chǎn)生。熔融池遷移初期的峰值對(duì)IVR包容失效的威脅需要得到關(guān)注。

圖4 熔融池釋熱功率變化規(guī)律

堆芯運(yùn)行40 d后，熔融池5 d內(nèi)以及1～2 d內(nèi)的衰變熱總功率變化趨勢(shì)與曲線極值處的變化規(guī)律分別示于圖4c、d。由圖4d可看出，最大釋熱功率出現(xiàn)在約18 MW處，隨后在第2燃耗步陡降于16 MW以下，又于第3燃耗步重新上升，在第0.2 d達(dá)到極值，再以指數(shù)衰減形式持續(xù)下降。極值的產(chǎn)生使釋熱衰減趨勢(shì)減緩，這反映熔融池中裂變產(chǎn)物的衰變釋熱主導(dǎo)了釋熱功率變化。

堆芯運(yùn)行30、70、130、200、250、300 d的熔融池釋熱功率與時(shí)間關(guān)系曲線示于圖5。該模型燃耗步設(shè)置為30步，前10步時(shí)間步長(zhǎng)為0.1 d，后20步時(shí)間步長(zhǎng)為0.2 d，模擬了5 d內(nèi)的衰變熱功率變化情況。圖5顯示，隨著運(yùn)行時(shí)間的增加，熔融池釋熱功率增加。裝料初期由于核燃料未完全消耗，導(dǎo)致所形成熔融池中裂變產(chǎn)物較少，衰變熱釋熱功率較低。另一方面，運(yùn)行250 d與運(yùn)行300 d的熔融池釋熱功率曲線完全重合，說(shuō)明熔融池內(nèi)核素組分在運(yùn)行250 d后已趨于穩(wěn)定，計(jì)算最大釋熱功率為22 MW。

圖5 運(yùn)行不同時(shí)間熔融池釋熱功率變化曲線

運(yùn)行300 d后熔融池截面內(nèi)部柵元釋熱功率密度分布隨時(shí)間的變化示于圖6。圖6顯示，熔融池形成初期(0.1 d)由于熔融池內(nèi)中子的影響，中心部分釋熱功率相對(duì)較高，有較大的功率密度；于0.3 d時(shí)各柵元功率密度尚有微弱差異，隨后趨于均勻并穩(wěn)定下降。這反映出隨著裂變產(chǎn)物的逐漸衰變，燃料裂變的釋熱影響可被忽略，釋熱逐漸被裂變產(chǎn)物衰變釋熱主導(dǎo)。輕金屬層中因存在部分輕裂變產(chǎn)物有少量衰變釋熱，功率密度相對(duì)于氧化物層與輕金屬層可忽略。

圖6 運(yùn)行300 d后熔融池釋熱功率密度變化規(guī)律

4 結(jié)論

本研究采用反應(yīng)堆蒙特卡羅軟件RMC及燃耗/衰變熱計(jì)算模塊DEPTH，對(duì)嚴(yán)重事故下熔融池內(nèi)熱源的變化規(guī)律進(jìn)行了計(jì)算，并基于得到的數(shù)據(jù)進(jìn)行曲線擬合，得到了時(shí)序模型。通過(guò)研究，得到如下結(jié)論。

1)考慮熔融池內(nèi)由于活化裂變碎片、自發(fā)裂變核素等釋放出中子，初期熱源分布呈中心高于四周的特性，體現(xiàn)了熔融池內(nèi)的中子輸運(yùn)影響熔融池核素組成，并主導(dǎo)熔融池形成初期的釋熱變化。

2)反應(yīng)堆運(yùn)行不同時(shí)刻形成的熔融池，釋熱功率變化為相似的指數(shù)衰減趨勢(shì)，其值則隨堆芯中裂變反應(yīng)程度不同而發(fā)生改變。全堆運(yùn)行250 d時(shí)，熔融池材料組成達(dá)到穩(wěn)定，最大釋熱功率達(dá)22 MW，裂變產(chǎn)物衰變主導(dǎo)了熔融池形成后期的釋熱變化。

3)以往的研究認(rèn)為，熔融池內(nèi)熱源僅存在于下封頭氧化物層中，而反應(yīng)堆運(yùn)行期間可能產(chǎn)生輕核素，如Zr的同位素存在于輕金屬層且放出衰變熱。輕金屬層作為熱源對(duì)熔融池傳熱、聚集效應(yīng)、下封頭壁面蠕變等產(chǎn)生的影響需進(jìn)一步驗(yàn)證。

本研究結(jié)果及時(shí)序模型可用于更精確地反映堆嚴(yán)重事故診斷及IVR性能模擬，為IVR的性能測(cè)試及設(shè)計(jì)提供數(shù)據(jù)支撐。