基于COMSOL的COPRA熔融池實(shí)驗(yàn)仿真研究

朱光昱 閔金坤 張 麗 元一單

1（中國(guó)核電工程有限公司核電安全研究中心 北京 100840）

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

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

熔融物堆內(nèi)滯留技術(shù)（In-vessel Retention，IVR）是目前廣泛采用的核反應(yīng)堆嚴(yán)重事故堆芯熔融物滯留方案，其有效性取決于反應(yīng)堆壓力容器下封頭壁面是否可以有效導(dǎo)出熔融物的衰變熱。為了研究反應(yīng)堆嚴(yán)重事故后壓力容器內(nèi)熔融池與冷卻壁面的換熱特性，國(guó)內(nèi)外設(shè)計(jì)了許多壓力容器下封頭二維切片 實(shí) 驗(yàn) 臺(tái) 架 ，例 如 COPO［1］（COrium Pool avec formation de glace aux parois）、BALI［2］（BAin Liquide）、SIMECO［3］（Simulation of In-vessel MElt COolability）以及由中國(guó)核電工程有限公司和西安交通大學(xué)針對(duì)華龍一號(hào)堆型搭建的COPRA 實(shí)驗(yàn)臺(tái)架［4?5］（COrium Pool Research Apparatus）。 其 中 ，COPRA實(shí)驗(yàn)采用了穆爾比為20%NaNO3-80%KNO3的二元混合物作為實(shí)驗(yàn)工質(zhì)，該二元混合物與反應(yīng)堆內(nèi)真實(shí)熔融物UO2-ZrO2的相圖存在相似性［4］，可以較為準(zhǔn)確地模擬下封頭熔融池內(nèi)的熔融物相變行為和高瑞利數(shù)自然對(duì)流工況，為IVR 相關(guān)系統(tǒng)設(shè)計(jì)提供了必要的技術(shù)儲(chǔ)備。由于熔融池內(nèi)的物理變化過程十分復(fù)雜，限制了COPRA 實(shí)驗(yàn)中測(cè)量手段，需要建立仿真計(jì)算模型來彌補(bǔ)實(shí)驗(yàn)的不足。對(duì)此，本文基于COMSOL Multiphysics多物理場(chǎng)計(jì)算軟件建立了一種熔融池計(jì)算模型，通過耦合湍流場(chǎng)和傳熱相變過程，實(shí)現(xiàn)了對(duì)COPRA 實(shí)驗(yàn)熔融池內(nèi)的流動(dòng)傳熱、相變過程的仿真模擬。

1 仿真模型

1.1 數(shù)學(xué)物理模型

采用COMSOL 軟件中的非等溫傳熱模塊耦合實(shí)驗(yàn)中的流體傳熱、相變和湍流流動(dòng)。穆爾比為20%NaNO3-80%KNO3的二元混合物具有較大固液相線溫差，適宜采用流體傳熱模塊中的相變材料模型進(jìn)行模擬，其原理如圖1所示，溫度低于固相線區(qū)域采用固相物性計(jì)算，溫度高于液相線區(qū)域采用液相物性計(jì)算，在固液相線之間采用固液兩相按體積平均的物性參數(shù)進(jìn)行計(jì)算。圖1中θ為相體積分?jǐn)?shù)；ΔT為固液相線溫差；Tpc為COMSOL中的相變溫度，這里設(shè)置為固相線溫度和液相線溫度的均值。參考Gaus-Liu［6］和 Janz［7］的 研 究 結(jié) 果 ，20%NaNO3-80%KNO3的二元混合物的固相線溫度為497.15 K，液相線溫度為557.15 K，由固相線到液相線的相變潛熱為 161 960 J·kg?1，其他主要物性參數(shù)如表 1所示。

圖1 相變模型原理圖Fig.1 Schematic diagram of phase change model

湍流流動(dòng)模塊采用低雷諾數(shù)k-ε模型進(jìn)行計(jì)算，參考相變材料模型設(shè)置思路，該模塊中的粘度項(xiàng)設(shè)置為固液兩項(xiàng)的體積平均值，可通過設(shè)置較大的粘度來模擬固相不存在流動(dòng)性的情況。通過預(yù)計(jì)算對(duì)比，當(dāng)固相的粘度系數(shù)高于103Pa·s后對(duì)計(jì)算的影響很小。

求解過程采用分離式求解器，時(shí)間步進(jìn)選擇向后差分，計(jì)算容差因子為0.5，計(jì)算加速性和穩(wěn)定性設(shè)置為Anderson加速度。

1.2 計(jì)算域設(shè)置

圖2 為參考COPRA 實(shí)驗(yàn)臺(tái)架的切片結(jié)構(gòu)建立的二維計(jì)算模型。臺(tái)架整體呈半徑為2.2 m的1/4圓形，加注熔鹽后臺(tái)架可分為底部的熔融池區(qū)域和頂部的輻射換熱區(qū)域，熔鹽加注的最大高度為1.9 m。其中1/4圓形半徑上的邊界（c、d和f）為絕熱邊界，圓弧邊界（a和b）外側(cè)為水冷流道，其平均溫度約維持在323.15 K，可以認(rèn)為是等溫邊界條件。

頂部的輻射換熱區(qū)域可認(rèn)為是4個(gè)表面組成的封閉腔輻射傳熱問題，其輻射換熱網(wǎng)格圖如圖3 所示，其中絕熱壁面c、d可認(rèn)為是重輻射面，則熔融池上表面（e）向水冷壁面輻射熱流Φeb可通過式（1）計(jì)算：

表1 20%NaNO3-80%KNO3混合物物性Table 1 Property of 20%NaNO3-80%KNO3

圖2 計(jì)算域和網(wǎng)格劃分示意Fig.2 Schematic diagram of calculation domain and mesh generation

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

圖3和式（1）中：Ebb、Ebc、Ebd、Ebe分別為各表面上的黑體輻射熱流密度；Jb、Jc、Jd、Je分別為各表面上的有效輻射熱流密度；ε1和ε2為熔融池上表面和臺(tái)架內(nèi)不銹鋼壁面的發(fā)射率，參考張盧騰［4］的成果分別取0.44 和 0.8；Ab、Ac、Ad、Ae分別為各表面的面積，對(duì)于二維模型，可通過將切片臺(tái)架的厚度進(jìn)行歸一化計(jì)算各面的面積；Xe，b、Xe，d、Xe，c、Xb，d、Xd，c、Xb，c為表面 e和b、e和d、e和c、b和d、d和c、b和c間的角系數(shù)，可將封閉腔簡(jiǎn)化為長(zhǎng)方形結(jié)構(gòu)通過查表［8］獲得；Rt為圖3中所示的各輻射熱阻，總等效熱阻需根據(jù)各熱阻串并聯(lián)關(guān)系計(jì)算，具體計(jì)算方法在文獻(xiàn)［8］中有詳細(xì)介紹。

通過引入黑體輻射公式可將式（1）化為式（2）的形式，此時(shí)可將總等效熱阻的倒數(shù)作為一個(gè)等效表面發(fā)射率考慮，則輻射換熱區(qū)域可簡(jiǎn)化為熔融池上表面對(duì)外界環(huán)境的輻射換熱過程，因此僅對(duì)由a、e和f 三個(gè)邊界圍成的熔融池區(qū)域的建立計(jì)算域即可。其中，將a 設(shè)置為邊界恒溫條件，e 設(shè)置為COMSOL軟件中的表面對(duì)環(huán)境輻射換熱邊界條件，f設(shè)置為絕熱邊界條件。

式中：σ為黑體輻射常數(shù)；Te為熔融池上表面的溫度；Tb為水冷壁面的溫度，這里設(shè)置為“表面對(duì)環(huán)境輻射換熱”邊界條件中的外界環(huán)境溫度。

1.3 網(wǎng)格敏感性分析

根據(jù)網(wǎng)格敏感性分析，當(dāng)網(wǎng)格數(shù)大于10 924后，網(wǎng)格變化對(duì)計(jì)算得到的溫度場(chǎng)分布和熔融物結(jié)殼厚度變化率降低至1%以內(nèi)，所以本文僅選取流場(chǎng)中的最大計(jì)算流速作為分析依據(jù)。圖4為加熱功率為10 kW時(shí)的網(wǎng)格敏感性分析結(jié)果，隨著網(wǎng)格加密，計(jì)算得到的流場(chǎng)中最大流速逐漸增加，在網(wǎng)格數(shù)加密至15 462后網(wǎng)格寬度對(duì)計(jì)算得到的最大速度變化率約為2%，因此最終計(jì)算采用的網(wǎng)格數(shù)為15 462。計(jì)算采用的初始時(shí)間步長(zhǎng)為1 ms，由于COMSOL軟件會(huì)自動(dòng)根據(jù)局部誤差對(duì)計(jì)算時(shí)間步長(zhǎng)進(jìn)行調(diào)整，因此未對(duì)時(shí)間步長(zhǎng)的敏感性進(jìn)行分析。

圖4 不同網(wǎng)格數(shù)下的最大計(jì)算速度Fig.4 Maximum simulation velocity under different mesh number

2 計(jì)算結(jié)果分析

COPRA 臺(tái)架通過在熔融池中的電加熱器模擬熔融物衰變熱，對(duì)于實(shí)驗(yàn)中的T1 工況［4］，加熱功率按照A 階段（功率15 kW，持續(xù)時(shí)間 23 700 s）、B 階段（功率10 kW，持續(xù)時(shí)間40 578 s）、C 階段（功率15 kW，持續(xù)至穩(wěn)態(tài)）變化。計(jì)算過程中，參照文獻(xiàn)［5］中的實(shí)驗(yàn)臺(tái)架設(shè)計(jì)信息，將不同的加熱功率折合成體積熱源并采用階梯函數(shù)控制，從而實(shí)現(xiàn)對(duì)實(shí)驗(yàn)過程的模擬。圖5為A階段熔融池的計(jì)算溫度穩(wěn)態(tài)分布，與實(shí)驗(yàn)結(jié)果類似，在熔融池內(nèi)的溫度總體呈現(xiàn)隨高度增加逐漸升高的趨勢(shì)，且不同高度出現(xiàn)了明顯的熱分層現(xiàn)象，在底部結(jié)殼的區(qū)域的溫度梯度較大，中部和頂部的溫度梯度較小。

圖5 A階段熔融池的熱分層結(jié)構(gòu)Fig.5 Thermal stratification of corium pool in stage A

圖6為各階段穩(wěn)態(tài)下的熔融池溫度分布的實(shí)驗(yàn)值和計(jì)算值。其中，底部區(qū)域的計(jì)算偏差相對(duì)較大，這可能是由于當(dāng)前計(jì)算模型沒有引入在熔融物凝固時(shí)會(huì)在硬殼和下封頭間形成氣體間隙導(dǎo)致的［4］。對(duì)于熔融池高度0.5 m以上區(qū)域，計(jì)算得到的溫度分布與實(shí)驗(yàn)值十分接近，最大計(jì)算誤差不超過1.5%。

圖6 各階段熔融池內(nèi)溫度穩(wěn)態(tài)分布Fig.6 Steady temperature profile in different stages of corium pool

圖7（a）為B 階段的流場(chǎng)，熔融池內(nèi)形成了由冷卻壁面向下流動(dòng)，在中部受熱后在浮力作用下緩慢向熔融池頂部和冷卻壁面流動(dòng)的自然對(duì)流趨勢(shì)。受自然對(duì)流和頂部輻射換熱的影響，在熔融池內(nèi)存在大量的渦流，這些渦流攪混的區(qū)域溫度分布趨于均勻，因此易于在熔融池內(nèi)形成熱分層結(jié)構(gòu)。圖7（b）為B階段的固態(tài)體積分?jǐn)?shù)，在冷卻壁面附近，受高速向下主流和重力的共同作用，凝固形成的糊狀物會(huì)更容易向熔融池底部轉(zhuǎn)移并最終結(jié)殼，從而使熔融物結(jié)殼的厚度呈現(xiàn)沿冷卻壁面由熔融池底部到頂部逐漸降低的趨勢(shì)。

圖7 B階段熔融池內(nèi)速度場(chǎng)(a)和固相體積分?jǐn)?shù)(b)Fig.7 Velocity field(a)and soild volume fraction(b)of corium pool in stage B

3 結(jié)語(yǔ)

采用COMSOL Multiphysics軟件的非等溫流動(dòng)模塊、相變材料模型建立的熔融池計(jì)算模型對(duì)COPRA實(shí)驗(yàn)進(jìn)行了仿真模擬，結(jié)果表明：

1）穩(wěn)態(tài)熔融池內(nèi)自然對(duì)流形式為：由冷卻壁面向下流動(dòng)，在中部受熱后在浮力作用下緩慢向熔融池頂部和冷卻壁面流動(dòng)。

2）在自然對(duì)流和頂部輻射換熱的共同影響下，熔融池內(nèi)會(huì)形成了大量的渦流。

3）在冷卻壁面向下主流和重力作用下，熔融物凝固硬殼呈現(xiàn)沿冷卻壁面由熔融池底部到頂部逐漸變薄的趨勢(shì)。