朱光昱，張寧娜，王昆鵬，石興偉，左嘉旭,*，劉福東

（1.生態(tài)環(huán)境部核與輻射安全中心，北京 100082；2.國家環(huán)境保護核與輻射安全審評模擬分析與驗證重點實驗室，北京 102488；3.中國核電工程有限公司核電安全研究中心，北京 100840）

壓水堆核電廠發(fā)生嚴重事故后，堆芯由于失去冷卻而熔化，由此形成的高溫熔融物會聚集在壓力容器下封頭內(nèi)形成熔融池[1]。對此，部分電廠設(shè)計了堆腔注水冷卻系統(tǒng)，通過在壓力容器外部注水的方式來冷卻下封頭，以避免壓力容器在嚴重事故后失效。該系統(tǒng)的有效性取決于熔融池施加到下封頭上的熱流密度是否低于冷卻水側(cè)的臨界熱流密度。為了給系統(tǒng)設(shè)計提供必要的依據(jù)，國內(nèi)外設(shè)計了多個實驗臺架并使用多種熔融物模擬物對下封頭內(nèi)熔融池換熱進行了研究，其中包括切片形式的COPO[2]（COrium Pool avec formation de glace aux parois）、BALI[3]（BAin Liquide）、SIMECO[4]（Simulation of In-vessel MElt COolability）和國內(nèi)的 COPRA 等實驗臺架[5]（COrium Pool Research Apparatus）以及下封頭三維縮比形式的 LIVE 實驗臺架（Late In-Vessel Phase Experiments）[6]。

為了準確模擬下封頭熔融池內(nèi)高瑞利數(shù)自然對流工況以及熔融物相變過程，LIVE 實驗臺架和 COPRA 實驗臺架采用了與原型熔融物UO2-ZrO2的相圖存在相似性的20 mol%NaNO3-80 mol%KNO3混合物作為實驗工質(zhì)[5,6]。由于熔融池內(nèi)的物理變化過程十分復雜，LIVE 實驗僅測量了熔融池中的部分區(qū)域的溫度值，出于彌補實驗中可采集數(shù)據(jù)局限性的目的，Zhang等[7]采用基于源項的SIMPLE 算法建立仿真計算模型，通過仿真手段研究了LIVE 實驗的熔融池流場和溫度分層情況。上述方法需要通過自編譯程序進行計算，為了獲得簡便的熔融池數(shù)值模擬方法，本文基于成熟的商用軟件COMSOL Multiphysics 建立了一種熔融池計算模型，對LIVE 實驗熔融池內(nèi)的溫度場和速度場進行了仿真模擬，測試了COMSOL 在熔融池仿真方面的適用性。

1 仿真模型

1.1 相變模型

20 mol%NaNO3-80 mol%KNO3的二元混合物的物性可根據(jù)Gaus-Liu[6]的研究結(jié)果設(shè)置，固相線溫度為497.15 K，液相線溫度為557.15 K，由固相線到液相線的相變潛熱為161 960 J·kg-1，其他主要物性參數(shù)如表1 所示。結(jié)合以往仿真經(jīng)驗等[5,7]研究結(jié)果，由于 20 mol%NaNO3-80 mol% KNO3的二元混合物具有較大固液相線溫差，采用相變模型可以更好地對流場進行模擬。COMSOL 中相變模型原理如圖1 所示，圖中θ為相體積分數(shù)，ΔT為固液相線溫差，固相線溫度和液相線溫度的均值Tpc即為COMSOL 中相變溫度，對于在固液相線之間的混合物，在計算過程中采用固液兩相按體積平均的物性參數(shù)進行計算。

圖1 相變模型原理圖Fig.1 The schematic of phase change

表1 20 mol% NaNO3-80 mol% KNO3 混合物物性Table 1 The property of 20 mol% NaNO3-80 mol% KNO3

根據(jù)Zhang等[5]研究，可采用低雷諾數(shù)k-ε模型對熔融池湍流場進行計算，混合物的粘度設(shè)置為固液兩相的平均值，本文通過設(shè)置較大的粘度來模擬固相不存在流動性的情況，當固相的粘度系數(shù)高于103Pa·s 后對計算的影響很小。

1.2 計算域設(shè)置

本次模擬工作針對LIVE-L5 實驗[8]進行。圖2 為LIVE 實驗臺架示意圖，其整體為典型壓水堆壓力容器下封頭的1/5，內(nèi)徑為1 m。實驗過程中通過熔鹽中的電加熱器模擬熔融物衰變熱，下封頭外部設(shè)有與 IVR（In-Vessel Retention，IVR）技術(shù)類似的冷卻流道。實驗中加注的熔鹽為 210 L，形成的熔融池高度為0.43 m。熔融池上部為熔鹽液面、水冷壁面和頂部蓋板組成的封閉腔體。腔體內(nèi)輻射換熱網(wǎng)格圖如圖3 所示，其中頂部壁面保溫性能很好可認為是重輻射面，則熔融池上表面向水冷壁面輻射熱流Φ可通過式（1）計算：

圖2 LIVE 實驗系統(tǒng)示意圖Fig.2 The schematic of the LIVE test vessel

圖3 輻射傳熱等效網(wǎng)絡(luò)圖Fig.3 The radiation heat transfer network

圖3 和式（1）中：Eb1、Eb2和Eb3以及J1、J2和J3分別為熔鹽液面、水冷壁面和頂部蓋板的黑體輻射熱流密度和有效輻射熱流密度；ε1和ε2為熔鹽液面和臺架內(nèi)不銹鋼壁面的發(fā)射率，參考以往仿真經(jīng)驗分別取0.44 和0.8[5]；A1、A2和A3為分別為熔鹽液面、水冷壁面和頂部蓋板的面積；X1,2、X1,3和X2,3為各表面間的角系數(shù)；Rt為圖3 中所示的各輻射熱阻的總等效熱阻，可以根據(jù)各熱阻串并聯(lián)關(guān)系計算，具體計算方法在文獻［9］中有詳細介紹。

式中：σ——黑體輻射常數(shù)；

T1——熔鹽液面的溫度；

T2——水冷壁面的溫度。

由于熔鹽的液位很高，此時水冷壁面在高度方向上的彎曲程度很小，因此可近似認為封閉腔體是扁平的圓柱體，并據(jù)此計算各輻射面之間的角系數(shù)。通過將黑體輻射公式帶入式（1）可導出式（2），此時將式（2）的總等效熱阻的倒數(shù)作為一個等效的表面發(fā)射率，水冷壁面的溫度作為外界環(huán)境溫度，便可將輻射換熱區(qū)轉(zhuǎn)化為熔融池上表面對外界環(huán)境的輻射換熱過程，這樣熔鹽液面可設(shè)置為COMSOL 軟件中的表面對環(huán)境輻射換熱邊界條件，從而不需要對輻射換熱的區(qū)域進行建模。此外，模型還考慮了頂部通入氮氣帶走的熱量，根據(jù)Zhang等[10]估算這部分提供的冷卻效率最大不超過11.7 W。

綜上所述，最終建立的軸對稱計算模型如圖4 所示。由于LIVE 實驗相關(guān)文獻中并未給出冷卻流道的具體尺寸，因此只能估測水冷壁面對流冷卻換熱系數(shù)。對此，本文根據(jù)實驗中的各穩(wěn)態(tài)工況中[8]下封頭壁面導出熱量占總加熱功率的比例來校核對流換熱邊界的熱通量。

圖4 計算域和網(wǎng)格劃分Fig.4 The calculation domain and mesh generation

1.3 邊界層和網(wǎng)格敏感性

參考以往的研究結(jié)果[7,10]，仿真得到熔融池內(nèi)流速均很低，因此邊界層變化對速度場中最大速度和平均速度影響很小。然而，通過預計算發(fā)現(xiàn)，邊界層的設(shè)置對水冷卻壁面附近相變計算的影響較為明顯，經(jīng)過測試，采用COMSOL 默認的邊界層設(shè)置，在水冷壁面內(nèi)設(shè)置8 層邊界層網(wǎng)格即可滿足一般的計算需要。圖5 所示為加熱功率為18 kW，初始化溫度為615 K 時，不同網(wǎng)格數(shù)下熔融池平均溫度隨時間的變化。隨著網(wǎng)格加密，計算得到的平均溫度逐漸降低，在網(wǎng)格數(shù)加密至4 121 后，進一步加密的計算結(jié)果變化率在1%以內(nèi)，因此最終計算采用的網(wǎng)格數(shù)為4 121。

圖5 不同網(wǎng)格數(shù)下的熔融池平均溫度Fig.5 The average corium pool in different mesh number

2 計算結(jié)果分析

圖6 和圖7 所示為熱功率為18 kW 時熔融池的計算溫度分布和速度場，在穩(wěn)態(tài)工況下，冷卻壁面附近熔鹽受冷后形成了沿壁面向下的自然對流，其他不同高度區(qū)域存在多個渦流。受渦流攪拌的熔融池中、上部溫度分布較為均勻，而底部則的呈現(xiàn)出了明顯的熱分層結(jié)構(gòu)。圖8 所示為熱功率為5 kW 和18 kW 時穩(wěn)態(tài)情況下，距離中心軸365 mm 處的熔融池溫度隨高度的變化。其中底部區(qū)域的計算值與實驗值[7]相比偏差較大，這可能是因為LIVE 實驗的下封頭中底部布置了大量的加熱和測溫設(shè)備，實際實驗條件下該區(qū)域的渦流攪拌作用可能并不明顯，因此溫度場仿真結(jié)果并不準確。對于熔融池中部和頂部區(qū)域，計算得到的溫度分布與實驗值十分接近，誤差不超過4%。

圖6 熱功率18 kW 時熔融池的熱分層結(jié)構(gòu)Fig.6 Thermal stratification of the corium pool under the heating power of 18 kW

圖7 熱功率18 kW 時熔融池內(nèi)速度場Fig.7 The velocity field of the corium pool under the heating power of 18 kW

圖8 熱功率為5 kW 和18 kW 時熔融池內(nèi)溫度穩(wěn)態(tài)分布Fig.8 The steady corium pool temperature profile under the heating power of 5 kW and 18 kW

圖9 和圖10 示出了熱功率為18 kW 時熔融池內(nèi)的結(jié)殼厚度分布以及結(jié)殼厚度的計算值和實驗值[7]。受重力和沿冷卻壁面向下的自然對流影響，凝固的熔鹽會向下封頭底部遷移，因此結(jié)殼厚度總體呈現(xiàn)隨著熔融池圓弧壁面徑向角度增加而減小的趨勢。在嚴重事故緩解過程中，堆內(nèi)熔融物結(jié)殼可以形成一層熱阻，從而有效阻止下封頭的不銹鋼壁面發(fā)生熔化。根據(jù)當前的實驗和仿真得到得結(jié)殼分布情況，即使僅考慮單層熔融池，下封頭的高相位角區(qū)域也屬于薄弱環(huán)節(jié)。當前熔融物堆內(nèi)滯留技術(shù)（Invessel Retention，IVR）普遍采用了從下封頭底部向堆腔注水的方式，由于冷卻流道的截面積隨著相位角增加而增大導致在高相位角區(qū)域冷卻水流速降低，同時冷卻水在底部受熱導致高相位角區(qū)域水溫相對較高，因此高相位角區(qū)域的冷卻效果相對底部較差，更容易發(fā)生熔穿等后果，在后續(xù)IVR 技術(shù)工程設(shè)計過程中，應考慮在下封頭外冷卻流道的高相位角區(qū)域增加強化換熱設(shè)計。

圖9 18 kW 時固相體積分數(shù)Fig.9 The soild volume fraction under the heating power of 18 kW

圖10 18 kW 時冷卻壁面上結(jié)殼厚度Fig.10 The crust thickness along the vessel wall under the heating power of 18 kW

3 結(jié)論

基于COMSOL Multiphysics 軟件建立了核電廠嚴重事故后下封頭熔融池的計算模型，采用LIVE 實驗結(jié)果對模型進行了驗證。仿真結(jié)果表明熔融池中除了沿冷卻壁面向下的自然對流外，熔融池的上、中部還同時存在大量渦流，受其攪拌的區(qū)域分布較為均勻，而熔融池底部黏度較大無法產(chǎn)生渦流因此呈現(xiàn)出了明顯的熱分層結(jié)構(gòu)。在穩(wěn)態(tài)下，仿真得到的中部和頂部熔融池溫度分布以及冷卻壁面附近的結(jié)殼厚度與實驗結(jié)果十分接近，說明當前COMSOL 模型適用于熔融池的仿真計算。