嚴重事故熔融池實驗熱工水力特性仿真分析

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

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

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

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

1 仿真模型

1.1 相變模型

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

圖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

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

1.2 計算域設置

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

圖2 LIVE 實驗系統示意圖Fig.2 The schematic of the LIVE test vessel

圖3 輻射傳熱等效網絡圖Fig.3 The radiation heat transfer network

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

式中：σ——黑體輻射常數；

T1——熔鹽液面的溫度；

T2——水冷壁面的溫度。

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

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

圖4 計算域和網格劃分Fig.4 The calculation domain and mesh generation

1.3 邊界層和網格敏感性

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

圖5 不同網格數下的熔融池平均溫度Fig.5 The average corium pool in different mesh number

2 計算結果分析

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

圖6 熱功率18 kW 時熔融池的熱分層結構Fig.6 Thermal stratification of the corium pool under the heating power of 18 kW

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

圖8 熱功率為5 kW 和18 kW 時熔融池內溫度穩態分布Fig.8 The steady corium pool temperature profile under the heating power of 5 kW and 18 kW

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

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

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

3 結論

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