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

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

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

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

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

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

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

1 計算方法

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

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

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

1.2 計算流程

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

圖1 熔融池分時刻釋熱功率計算流程圖

1.3 衰變熱計算

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

(1)

(2)

(3)

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

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

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

(4)

(5)

(6)

針對反應(yīng)堆停堆以及熔融池等場景中常見的衰變熱計算需求，DEPTH模塊內(nèi)置了衰變熱計算功能[5]?；谥凶幼⒘颗c絕對原子密度直接計算燃耗區(qū)內(nèi)的衰變熱釋熱功率Q，即：

Q=Rq=Nλq

(7)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.1 反應(yīng)堆運行時刻選取

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

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

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

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

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

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

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

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

圖5 運行不同時間熔融池釋熱功率變化曲線

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

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

4 結(jié)論

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

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

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

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

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