朱光昱閔金坤靖劍平王昆鵬劉福東

1（生態(tài)環(huán)境部核與輻射安全中心北京100082）

2（國家環(huán)境保護核與輻射安全審評模擬分析與驗證重點實驗室北京102488）

3（中國核電工程有限公司北京100840）

4（清華大學工程物理系北京100084）

核電廠發(fā)生嚴重事故后，堆芯各類構件由于失去冷卻后熔化，逐漸在壓力容器下封頭內(nèi)形成了Fe、Zr和UO2-ZrO2組成的高溫熔融池。為防止事故進一步發(fā)展造成大規(guī)模放射性釋放后果，第三代核電廠設計過程中均進行了堆芯熔融物滯留設計。主要包括：將熔融物滯留在壓力容器內(nèi)的堆內(nèi)滯留（In-vessel Retention，IVR）技術［1］，以及在壓力容器破損后將熔融物重新收集冷卻的堆外滯留（Exvessel Retention，EVR）技術［2］。我國的先進壓水堆設計更傾向于采用IVR技術，在第三代反應堆（如HPR1000）以及隨后的新型設計中均使用了該策略［3］。

在堆芯熔融池形成后，由于輕金屬與氧化物并不相溶，金屬層會浮于氧化物層之上形成雙層熔融池結(jié)構［4］。在此基礎上，由于部分析出的U與Zr混合可能會形成重金屬層，而形成三層熔融池結(jié)構。然而，根據(jù)Bechta等［5］基于熔融池中U-Zr-Fe-O的相平衡機理分析結(jié)果，三層熔融池僅屬于過渡段形態(tài)，雙層結(jié)構是熔融池發(fā)展的最終形態(tài)。由于輕金屬層的熱聚集效應會使反應堆壓力容器（Reactor Pressure Vessel，RPV）下封頭更容易發(fā)生熔穿，在IVR設計過程中必須考慮雙層熔融池對壓力容器完整性的挑戰(zhàn)。由于真實熔融物的大尺度熔融池實驗存在一定困難，國內(nèi)外一些學者采用數(shù)值模擬手段對雙層［6-7］或三層［8-9］熔融池的熱工水力特性、金屬層的聚焦效應和下封頭機械性能進行了研究。本文基于當前大型三代壓水堆的設計參數(shù)，采用COMSOL Multiphysics多物理場耦合軟件建立了一個雙層熔融池仿真模型，對RPV下封頭內(nèi)雙層熔融池的換熱特性進行了數(shù)值模擬研究。

1 仿真模型

1.1 計算域設置

圖1為采用四邊形網(wǎng)格建立的二維軸對稱計算模型，其中熔融池半徑為2.2 m，下封頭厚度為0.15 m，下封頭頂部區(qū)域延伸了長度為0.4 m的不銹鋼壁面。仿真過程中，溫度場采用固體、流體傳熱模型計算，流場采用低雷諾數(shù)k-ε模型［4，10］計算，通過非等溫流動模塊進行多物理場耦合。

圖1 計算域和網(wǎng)格劃分示意圖Fig.1 Diagram of calculation domain and mesh generation

計算域中，氧化物層高度為1.58 m，金屬層高度為0.56 m，金屬層頂部和下封頭延伸部分的內(nèi)壁面均設置為表面對環(huán)境輻射條件，發(fā)射率分別設置為0.45和0.8［1］。目前，大型壓水堆采用的堆腔注水冷卻系統(tǒng)均是從下封頭底部注入冷卻水，所以沿著下封頭弧形面由下至上的冷卻水流速是逐漸降低，因此循環(huán)流速一般采用通道內(nèi)平均流速表征。在AP1000堆型以及國內(nèi)自主研發(fā)的國和一號堆型中，堆內(nèi)熔融物滯留系統(tǒng)的冷卻水循環(huán)流量分別為3 000 t·h-1和3 800 t·h-1，結(jié)合流道尺寸計算平均流速分別為0.5 m·s-1和0.6 m·s-1。結(jié)合上述參考，將下封頭外側(cè)設置為強制對流換熱邊界條件，冷卻水側(cè)流速取0.5 m·s-1，過冷度設置為50 K。

在熔融物穩(wěn)定分層后，可以認為衰變熱全部集中在氧化物層中，其內(nèi)熱源可采用式（1）計算：

式中：Q為堆芯衰變熱功率，MW；fd為衰變熱修正系數(shù)，可設置為0.75［11］；P0可認為是反應堆正常運行時的熱功率，本文取3 050 MW；t為計算時長，s。

考慮到事故初期并未形成雙層熔融池，仿真過程中，由事故發(fā)生4 h后開始進行計算，氧化物層的初始體積釋熱率約為1.49 MW·m-3。

1.2 物性設置

熔融池內(nèi)通過引入相變模型模擬凝固過程，材料密度采用Boussinesq近似［4］處理。其中，氧化物層參考溫度為2 650 K，參考密度為8 196.41 kg·m-3，完全凝固后的密度和導熱系數(shù)設置為9 512.8 kg·m-3和2.8 W·m-1·K-1。金屬層參考溫度為1 750 K，參考密度為6 818.50 kg·m-3。根據(jù)VULCANO熔融物鋪展實驗相關數(shù)值模擬經(jīng)驗［12］，可以采用式（2）計算各層熔融池中糊狀區(qū)的黏度。

式中：μl為液相黏度，氧化物層和金屬層分別取0.008 12 Pa·s和0.003 01 Pa·s；系數(shù)C一般取3.5～8［13］；θs為固相體積分數(shù)，可直接在計算過程中實時調(diào)用參數(shù)。

熔融池內(nèi)其他計算所需物性參照文獻［4］和［14］設置，具體如表1所示。冷卻水參數(shù)直接引用COMSOL數(shù)據(jù)庫中的物性參數(shù)，不考慮實際工況下水中含有硼酸帶來的影響。

表1 物性參數(shù)Table 1 Physical property

計算過程中依然采用相變材料來模擬下封頭熔化過程，熔化前下封頭的材料物性設置為A508-3鋼的物性，熔點為1 600 K，導熱系數(shù)為32 W·m-1·K-1，相變潛熱為269.55 kJ·kg-1，并設置了20 K的相變溫度區(qū)間。由于實際計算中下封頭與氧化物層接觸區(qū)域沒有發(fā)生熔化，與金屬層接觸區(qū)域熔化后的鋼相對于金屬層總質(zhì)量而言很小，可以認為下封頭熔化后的材料物性參數(shù)與金屬層一致。在當前COMSOL計算模型中，下封頭部分設置為固體域，無法模擬材料熔化后混入臨近熔融池的過程，為了考慮下封頭熔化區(qū)域內(nèi)流動導致的換熱，在材料本身的導熱系數(shù)之外將湍流導熱添加至熔化后區(qū)域的導熱系數(shù)當中，湍流導熱可采用式（3）計算。

式中：kT為湍流流動產(chǎn)生的導熱系數(shù)，W·m-1·K-1；μ為湍流動力黏度，Pa·s，可在用戶定義項中調(diào)用鄰近金屬層內(nèi)的參數(shù)作為參考。PrT為湍流普朗特數(shù)，對于一般的湍流邊界層問題取0.9。Cp為比熱容，參照熔化后的金屬層參數(shù)設置。

初始化過程中，氧化物層、金屬層以及下封頭的初始化溫度分別設置為2 650 K、1 750 K和1 000 K。

1.3 網(wǎng)格無關性分析

經(jīng)預計算測試，網(wǎng)格寬度大于1.7 cm后計算會直接發(fā)散，而網(wǎng)格數(shù)小于1.5 cm時，僅當計算時長較短時可以獲得收斂結(jié)果。因此僅在網(wǎng)格寬度1.5～1.7 cm內(nèi)進行網(wǎng)格無關性測試。圖2為不同網(wǎng)格數(shù)下計算時長為0.5 h時，模型中軸線上的溫度分布情況。網(wǎng)格數(shù)對氧化物層的溫度分布影響很小，但對金屬層溫度分布影響較大。隨著網(wǎng)格加密，金屬層處的溫度差異逐漸減小，本文最終采用的網(wǎng)格數(shù)為19 603的模型進行后續(xù)的仿真工作。

圖2 不同網(wǎng)格數(shù)下的模型中軸溫度分布Fig.2 Temperature distribution on the axis of calculation domain under different mesh number

2 計算結(jié)果分析

圖3給出了計算時長為1 h、1.5 h和3 h的熔融池速度場。兩層流體由于密度差異較大形成了穩(wěn)定的分層狀態(tài)，氧化物層在衰變熱的作用下很快形成了自然循環(huán)，熔融物在靠近冷卻壁面一側(cè)受到冷卻向下流動，在主流中受熱后轉(zhuǎn)向上流動。初始階段金屬層不存在穩(wěn)定的循環(huán)流動狀態(tài)，在約1.5 h時逐漸建立了自冷卻壁面一側(cè)向下、在吸收金屬層熱量后向上的自然循環(huán)流動。然而，與氧化物層始終存在一個穩(wěn)定的自然循環(huán)流動不同，在頂部輻射換熱冷卻和底部氧化層流動剪切作用的共同影響下，除了冷卻壁面附近存在一個較穩(wěn)定的自然循環(huán)外，金屬層中還同時存在多個不斷變動的渦流，導致流場不斷發(fā)生變動。

圖3熔融池內(nèi)的速度場(a)1 h，(b)1.5 h，(c)3 hFig.3 Velocity field of RPV corium pool(a)1 h,(b)1.5 h,(c)3 h

圖4 展示了當計算時長為2 h時，下封頭和熔融池內(nèi)的溫度場。在自然循環(huán)流動的攪拌作用下，氧化物層和金屬層的溫度近似一致且分布較為均勻，沒有出現(xiàn)明顯的熱分層現(xiàn)象。在金屬層熱聚集效應作用下，與其接觸的下封頭溫度明顯升高，而與氧化物層接觸的下封頭溫度則相對較低。

圖4下封頭和熔融池溫度分布Fig.4 Temperature distribution of RPV lower head and corium pool

圖5 展示了當計算時長為1 h時，各計算域內(nèi)的固相體積分數(shù)，此時氧化物層已經(jīng)進入到UO2-ZrO2糊狀區(qū)狀態(tài)而金屬層則依然保持純液態(tài)狀態(tài)。此外，氧化物層中的熔融物會迅速在下封頭處結(jié)殼，在低導熱系數(shù)硬殼的良好保護和外側(cè)冷卻水的共同作用下，該區(qū)域的下封頭并未發(fā)生明顯的熔化。與金屬層中接觸的下封頭區(qū)域由于溫度超過熔點發(fā)生了明顯的熔化，但最薄的區(qū)域的厚度依然大于10 cm。

圖5 下封頭和熔融池內(nèi)固相體積分數(shù)Fig.5 Soild volume fraction of RPV lower head and corium pool

圖6下封頭外壁面溫度分布Fig.6 Temperature distribution on outer wall of RPV lower head

圖6為下封頭計算時間1～3 h時下封頭外壁面的溫度分布，可見對于雙層熔融池結(jié)構，將下封頭外壁面考慮為強制對流冷卻邊界條件時下封頭外壁面溫度分布并不均勻，與金屬層接觸的下封頭區(qū)域溫度會顯著升高。圖7為計算時間1 h時下封頭計算域內(nèi)的熱流密度分布，金屬層導入的熱量除了向冷卻水側(cè)徑向傳遞之外，還會沿著下封頭不銹鋼壁面進行周向傳遞。對比張盧騰等［4］將冷卻壁面設置為400 K恒溫壁面得到的仿真結(jié)果，在徑向角度50°以下的區(qū)域，外壁面處的熱流密度均在0.15～0.20 MW·m-1。然而，在與金屬層接觸的區(qū)域，當前模型計算得到的熱流密度顯著低于張盧騰［4］的計算結(jié)果。

圖7 下封頭計算域內(nèi)熱流密度分布Fig.7 Distribution of heat flux in RPV lower head calculation domain

根據(jù)圖6所示溫度分布，與金屬層接觸的下封頭外側(cè)應處于沸騰狀態(tài)，因此采用強制對流邊界條件計算得到的換熱量會偏低。綜合上述分析，對于下封頭外部的冷卻壁面而言，采用單一的恒溫壁面模型或強制對流換熱模型進行仿真仍具有很大的局限性，相關仿真模型尚有改進空間。

3 結(jié)語

采用COMSOL軟件建立了一個熔融池仿真模型，對核電廠嚴重事故后壓力容器下封頭內(nèi)雙層熔融池的演變進行了仿真研究，結(jié)論如下：

1）在穩(wěn)定分層情況下，氧化物層會迅速形成存在穩(wěn)定的自然循環(huán)。金屬層中除靠近冷卻壁面的區(qū)域存在穩(wěn)定自然循環(huán)外，在其他區(qū)域還存在一些不穩(wěn)定的渦流。

2）在自然循環(huán)和渦流的攪渾作用下，氧化物層和金屬溫度分布較為均勻，不會發(fā)生熱分層現(xiàn)象。

3）氧化物層中，熔融物會迅速在冷卻壁面附近結(jié)殼，從而防止該區(qū)域下封頭熔化。由于金屬層的熱聚集效應，與其接觸的下封頭會發(fā)生明顯的熔化，但在良好的冷卻條件下，不會發(fā)生熔穿現(xiàn)象。

作者貢獻聲明朱光昱：采集數(shù)據(jù)和文章撰寫；閔金坤、靖劍平、王昆鵬：文章撰寫；劉福東：工作支持。