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

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

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

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

4（清華大學(xué)工程物理系北京100084）

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

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

1 仿真模型

1.1 計算域設(shè)置

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

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

計算域中，氧化物層高度為1.58 m，金屬層高度為0.56 m，金屬層頂部和下封頭延伸部分的內(nèi)壁面均設(shè)置為表面對環(huán)境輻射條件，發(fā)射率分別設(shè)置為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è)設(shè)置為強(qiáng)制對流換熱邊界條件，冷卻水側(cè)流速取0.5 m·s-1，過冷度設(shè)置為50 K。

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

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

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

1.2 物性設(shè)置

熔融池內(nèi)通過引入相變模型模擬凝固過程，材料密度采用Boussinesq近似［4］處理。其中，氧化物層參考溫度為2 650 K，參考密度為8 196.41 kg·m-3，完全凝固后的密度和導(dǎo)熱系數(shù)設(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í)驗(yàn)相關(guān)數(shù)值模擬經(jīng)驗(yàn)［12］，可以采用式（2）計算各層熔融池中糊狀區(qū)的黏度。

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

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

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

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

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

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

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

經(jīng)預(yù)計算測試，網(wǎng)格寬度大于1.7 cm后計算會直接發(fā)散，而網(wǎng)格數(shù)小于1.5 cm時，僅當(dāng)計算時長較短時可以獲得收斂結(jié)果。因此僅在網(wǎng)格寬度1.5～1.7 cm內(nèi)進(jìn)行網(wǎng)格無關(guān)性測試。圖2為不同網(wǎng)格數(shù)下計算時長為0.5 h時，模型中軸線上的溫度分布情況。網(wǎng)格數(shù)對氧化物層的溫度分布影響很小，但對金屬層溫度分布影響較大。隨著網(wǎng)格加密，金屬層處的溫度差異逐漸減小，本文最終采用的網(wǎng)格數(shù)為19 603的模型進(jìn)行后續(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)外，金屬層中還同時存在多個不斷變動的渦流，導(dǎo)致流場不斷發(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 展示了當(dāng)計算時長為2 h時，下封頭和熔融池內(nèi)的溫度場。在自然循環(huán)流動的攪拌作用下，氧化物層和金屬層的溫度近似一致且分布較為均勻，沒有出現(xiàn)明顯的熱分層現(xiàn)象。在金屬層熱聚集效應(yīng)作用下，與其接觸的下封頭溫度明顯升高，而與氧化物層接觸的下封頭溫度則相對較低。

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

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

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

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

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

3 結(jié)語

采用COMSOL軟件建立了一個熔融池仿真模型，對核電廠嚴(yán)重事故后壓力容器下封頭內(nèi)雙層熔融池的演變進(jìn)行了仿真研究，結(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ū)域下封頭熔化。由于金屬層的熱聚集效應(yīng)，與其接觸的下封頭會發(fā)生明顯的熔化，但在良好的冷卻條件下，不會發(fā)生熔穿現(xiàn)象。

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