基于OpenFOAM的熔融池自然對流傳熱與凝固數(shù)值研究

王 溪，孟召燦，程 旭（國核（北京）科學(xué)技術(shù)研究院有限公司，北京 100029）

基于OpenFOAM的熔融池自然對流傳熱
與凝固數(shù)值研究

王 溪，孟召燦，程 旭
（國核（北京）科學(xué)技術(shù)研究院有限公司，北京 100029）

摘要：熔融物堆內(nèi)滯留是第3代核電技術(shù)重要的嚴重事故緩解措施之一，堆芯熔融池在壓力容器下封頭壁面的熱流密度分布直接影響該策略的有效性。本文基于開源的數(shù)值計算流體力學(xué)軟件平臺Open－FOAM，應(yīng)用相變模型和浮升力模型二次開發(fā)了用于模擬堆芯熔融物由內(nèi)熱源或溫差驅(qū)動的自然對流傳熱與相變求解器。應(yīng)用該求解器模擬了瑞典皇家理工學(xué)院開展的二維氧化池與金屬層耦合傳熱試驗，獲得了氧化池和金屬層硬殼的相場，以及熔融池內(nèi)的溫度分布及沿容器壁面的熱流密度分布。計算結(jié)果表明，該模型可用于熔融物凝固與自然對流的模擬，為深入分析核電廠采用熔融物堆內(nèi)滯留措施后熔融池的行為奠定了基礎(chǔ)。

關(guān)鍵詞：嚴重事故；CFD；熔融物堆內(nèi)滯留；凝固；自然對流

通過壓力容器外部冷卻滯留堆芯熔融物（IVR－ERVC）是第3代核電技術(shù)重要的嚴重事故緩解措施之一，已應(yīng)用于AP1000、CAP1400等堆型。該措施是在嚴重事故下，使用冷卻劑淹沒壓力容器，熔化后的堆芯重定位到壓力容器下封頭內(nèi)，衰變熱通過熔融池自然對流和頂部輻射傳遞給壓力容器，最終通過壓力容器壁的導(dǎo)熱傳遞給外部冷卻劑［1］。

美國核管理委員會和西屋公司在AP600／1000的研究中，結(jié)合ACOPO、UCLA和MELAD等試驗提出的經(jīng)驗關(guān)系式，采用了集總參數(shù)法分析兩層熔池的熱流密度分配，評價了IVR的有效性［2］。德國卡爾斯魯厄理工大學(xué)（KIT）完成了以熔鹽為模擬材料1∶5比例三維LIVE試驗，重點研究熔融物進入下封頭的瞬態(tài)熱沖擊［3］。瑞典皇家工學(xué)院（KTH）完成了以熔鹽為模擬材料1∶8比例二維切片熔融池金屬層和氧化層耦合相變傳熱SIMECO試驗，機理上深入探索了熔融池傳熱行為，并采用了相變等效傳熱模型（PECM）模擬了SIMECO試驗，但該方法不能用于分析熔融池的流場［4］。在韓國的APR1400熔融物堆內(nèi)滯留的研究中，韓國原子能研究院等采用CFD程序LILAC模擬了氧化池和金屬層的傳熱，計算中假設(shè)氧化池熔融物與硬殼的交界面、金屬層與未熔化的壓力容器壁的交界面，以及金屬層上表面均為等溫壁面邊界條件，僅模擬了熔融池內(nèi)液相區(qū)域的流動和傳熱［5］。本文基于開源計算流體力學(xué)平臺OpenFOAM，將用于模擬凝固相變的Enthalpy－porosity模型和浮升力Buoyancy模型應(yīng)用到該軟件中，并應(yīng)用該軟件模擬KTH所開展的基于模擬材料的二維氧化池和金屬層耦合傳熱試驗，以初步驗證該數(shù)值模型在準穩(wěn)態(tài)熔池傳熱研究方面的可用性。

1　凝固相變數(shù)值模型

1．1 基本控制方程

針對熔融物在壓力容器下封頭內(nèi)的流動和相變問題，在假設(shè)堆芯熔融物為不可壓縮流體的情況下，可建立以下方程。

連續(xù)性方程：

Δ·U＝0（1）

其中，U為流體速度矢量。

動量方程：

Sb＝－ρgβmax［T－Tl，0］（3）

其中：t為時間；μeff為有效黏性系數(shù)；p為壓力；ρ為密度；g為重力加速度；β為熱膨脹系數(shù)；T為流體溫度；Tl為液相線溫度。動量方程包括兩個源項：Sb為浮升力項，當溫度大于液相線溫度時，將溫度與液相線溫度的差值乘以材料的熱膨脹系數(shù)、密度和重力加速度，作為材料的浮升力，方向與重力方向相反，當溫度小于液相線溫度時，認為浮升力消失；Sm為描述介于固相與液相之間的固液兩相過渡區(qū)域類似于一種多孔介質(zhì)對流動產(chǎn)生的阻力，稱為Darcy項。

能量方程：

其中：h為焓；λeff為有效導(dǎo)熱系數(shù)；能量方程的源項是體積內(nèi)熱源Q，用于模擬衰變熱。

1．2 Enthalpy－porosity模型

Enthalpy－porosity模型是工程上應(yīng)用最為成熟的凝固模型，該模型由Voller等［6］提出。在該模型中，動量方程Darcy項由Darcy定律［7］推導(dǎo)而來，Darcy定律方程如下：

其中：K為滲透率；μ為流體的黏度；Δp為壓力梯度。Darcy定律描述了流體以一定速度流過多孔介質(zhì)所產(chǎn)生的壓力下降，或說是阻力。

Darcy方程中的K由Carmen－Kozeny方程獲得［8］，假設(shè)多孔介質(zhì)由若干圓孔通道組成，細孔的半徑為r，空隙體積比為φ，多孔介質(zhì)的長度為L，通道平均實際長度為La，則K為：

K＝φr2L／8La（6）

如果定義單位體積多孔介質(zhì)內(nèi)空隙的表面積為Spv，則K為：

K＝φL／2S2pvLa（7）

Spv＝S／Vp（8）

假設(shè)固液兩相過渡區(qū)域中凝固的固相材料晶體為等直徑的圓球形顆粒狀，熔融的液相的材料在這些固相顆粒之間流動，液相的體積份額即為該多孔介質(zhì)的空隙體積比φ。單位體積內(nèi)固液相接觸的表面積為：

將式（9）代入式（7），可得固液兩相過渡區(qū)域的K為：

將式（10）代入式（6）可得液相速度與多孔介質(zhì)壓降的關(guān)系：

式（11）變形可得：

式（12）的表達為壓力梯度，這與動量方程中的壓力項相當，因此可直接得到動量方程的Darcy源項表達式為：

其中：dp為固相顆粒的平均直徑；ε為一極小數(shù)，防止除零；C為系數(shù)的整合。不難看出，Darcy項在動量方程中為一與速度呈正比的阻力項，且與液相體積份額相關(guān)。實際使用中C的取值越大，則計算結(jié)果的固液交界區(qū)域越薄，相反則越彌散。對于C的選取，需通過真實材料的試驗獲得。本文參考文獻［6］，選取C＝1．6×103kg／（s·m－3）進行計算。

其中：erf為誤差函數(shù)；Ts為固相線溫度；Tm為固液相線溫度的平均值。流體的溫度在固液相線溫度之間，隨溫度的升高，φ增大。溫度高于液相線溫度時，φ＝1；溫度低于固相線溫度時，φ＝0。

能量方程中的焓h與材料的溫度和液相體積份額相關(guān)，如式（15）所示。

h＝cpT＋L（15）

將式（15）代入式（4），可得式（16）。

方程右端多出了一個由相變引起的潛熱的釋放和吸收的源項，如式（17）所示。

將式（14）代入式（17），可得：

上式可認為是關(guān)于溫度T的能量方程右端的另外一個由相變引起的源項，由于引入了誤差函數(shù)，液相體積份額對時間和空間上的一階導(dǎo)數(shù)具有連續(xù)性。

以上控制方程考慮了相變、浮升力和內(nèi)熱源的作用，將這些方程應(yīng)用于OpenFOAM平臺，可求解凝固相變條件下的內(nèi)熱源驅(qū)動自然對流傳熱。

2　數(shù)值模型驗證

2．1 SIMECO試驗

SIMECO是KTH所開展的針對熔池傳熱研究的二維半圓形切片試驗裝置（圖1）。該試驗包括了氧化池和金屬層耦合傳熱試驗，可用于研究氧化池和金屬層硬殼界面對傳熱的影響，初步探索金屬層熱聚焦效應(yīng)［10］。

圖1　 SIMECO試驗裝置示意圖［10］Fig．1 Scheme of SIMECO experiment setup［10］

如圖1所示，該試驗裝置包括1個半圓形切片銅制試驗段，其直徑為62．0cm，高度為53．0cm，厚度為9．0cm。該試驗段的尺寸最高可達到的瑞利數(shù)Ra的范圍為1012～1013。試驗容器壁厚為2．3cm，由冷卻回路進行冷卻，該回路分為兩部分：一部分用于冷卻試驗容器的頂蓋，另一部分用于冷卻試驗容器的底部，冷卻系統(tǒng)主要是為試驗段提供定溫邊界條件。模擬試驗段內(nèi)熱源的加熱系統(tǒng)為1根直徑為3mm、長為4m的電阻加熱絲，最高可達4kW的加熱功率。

KTH在該裝置上開展了研究氧化池和金屬層自然對流傳熱?；囼灒ú捎?0%NaNO3－50%KNO3作為模擬材料的熔池試驗、用銅板將模擬材料水分隔的分層傳熱試驗、石蠟和水的分層傳熱試驗、較高熔點的硝酸鹽與Cerrobend合金的分層試驗以及甘油和Cerrobend合金的無硬殼界面耦合試驗等。

選取SIMECO系列試驗中的硝酸鹽和Cerrobend合金耦合相變傳熱試驗，開展了CFD模擬和模型驗證。該試驗采用NaNO3－KNO3混合物作為氧化池的模擬材料，采用低熔點Cerrobend合金（由鉛、鉍、錫和鎘等材料按一定比例混合而成）作為金屬層的模擬材料。模擬材料的物性列于表1。

表1　 模擬材料NaNO3－KNO3和Cerrobend合金物性Table 1 Material properties of NaNO3－KNO3and Cerrobend alloy

圖2　 SIMECO試驗容器內(nèi)壁溫度分布擬合Fig．2 Interpolation of vessel inner wall temperature distribution of SIMECO experiment

該試驗金屬層和氧化池是由1層薄銅板完全分隔開，由于銅板的傳熱性能很好，可忽略其熱阻。上部的金屬層無內(nèi)熱源，氧化池和金屬層由銅板耦合傳熱。試驗段金屬層和氧化池的高度比分為兩種：3∶27和6∶27。試驗中金屬層頂部的邊界條件采用了4種不同的邊界，包括頂部溫度為6、40、86℃定溫邊界條件和頂部絕熱邊界條件，分別研究較厚金屬層硬殼、較薄金屬層硬殼和無頂部硬殼的情況。共開展了14組試驗，本文選取了9號試驗開展數(shù)值模擬和模型驗證，工況和邊界條件為：輸入功率，2kW；頂部溫度，40℃；側(cè)面溫度，6℃；金屬層厚度，6cm。SIMECO試驗容器底部溫度分布擬合如圖2所示，作為CFD模型的輸入底部壁面邊界條件。其中，0°表示容器底部的溫度，90°表示容器橫向赤道位置的溫度。

圖3　 SIMECO 9號試驗數(shù)值模型Fig．3 Numerical model of SIMECO No．9experiment

2．2 SIMECO試驗數(shù)值模擬

結(jié)合上面給出的初始條件和邊界條件，本文應(yīng)用基于OpenFOAM二次開發(fā)的CFD求解器對SIMECO 9號試驗進行了數(shù)值模擬，數(shù)值模型如圖3所示。氧化池的電阻絲加熱等效為體積內(nèi)熱源，自然對流的作用下將熱量傳遞給上方的金屬層和底部的圓弧形容器壁面；金屬層底部獲得氧化池傳遞的熱量等效為定熱流密度邊界條件，在對流的作用下將這一部分熱流傳遞給頂部和側(cè)壁面。本文依據(jù)文獻［5］對于自然對流換熱系數(shù)網(wǎng)格相關(guān)性的分析結(jié)論，選取網(wǎng)格分辨率為800。計算采用了Spallart－Allmaras湍流模型。由于自然對流屬于準穩(wěn)態(tài)流動，首先基于穩(wěn)態(tài)SIMPLE算法獲得一初步的穩(wěn)態(tài)結(jié)果，然后將計算結(jié)果作為初始條件，采用瞬態(tài)PIMPLE算法獲得一滿足能量平衡的準穩(wěn)態(tài)結(jié)果。

9號試驗的數(shù)值模擬速度場如圖4所示。圖中流場分隔為流體區(qū)域和固相區(qū)域，靠近冷卻壁面的是固相區(qū)域，流體區(qū)域內(nèi)的箭頭是速度矢量。金屬層在下方氧化池的傳熱作用下，形成了類似于Bernard對流的流型，無數(shù)個小漩渦將熱量傳遞給頂部和側(cè)面未熔化的區(qū)域。由于金屬層側(cè)面得到了有效冷卻，硬殼形成的形狀為中間區(qū)域薄，靠近金屬層側(cè)壁位置較厚。受金屬層自然對流漩渦的影響，硬殼與熔融物交界面為波浪形起伏。從氧化池的速度場和液相體積份額分布可看出，氧化池上壁面的硬殼很薄，這是由較高的熱流密度引起的，沿著徑向方向，氧化池頂部的硬殼厚度較為均勻，這說明了氧化池向上傳遞的熱流密度分布相對均勻。氧化池半圓弧形壁面的硬殼越靠近底部區(qū)域越厚，這是由氧化池向容器壁傳熱的熱流密度分布引起的，底部的流體溫度較低，同時熱流密度也較低；熔池靠近半圓形弧面上方的位置，由于以溫度差引起的自然對流換熱為主，熱流密度較高，因此硬殼的厚度較薄。從氧化池的流場可看出，在氧化池上半部分空間區(qū)域內(nèi)形成了明顯的大漩渦，這些漩渦不斷地將內(nèi)熱源產(chǎn)生的熱量傳遞給頂部和底部壁面，然后通過硬殼的導(dǎo)熱載出熱量。CFD計算的結(jié)果直觀呈現(xiàn)出熔池內(nèi)硬殼幾何形狀以及熔融物的流動行為。

圖4　 SIMECO 9號試驗數(shù)值模擬速度場Fig．4 Velocity field of numerical simulation for SIMECO No．9experiment

圖5為SIMECO 9號試驗的數(shù)值模擬溫度場。從金屬層溫度場可看出，在兩個漩渦的中間均形成一股熱流，將從底部獲得的熱量傳遞給未熔化的金屬層材料，受到冷卻后下降。金屬層的硬殼區(qū)域的溫度梯度明顯高于下半部分對流區(qū)，因為對流換熱速率遠高于導(dǎo)熱速率。從氧化池的溫度場可看出，上半部分在大漩渦的作用下，溫度分布非常均勻，而下半部分出現(xiàn)了明顯的熱分層現(xiàn)象，越靠近底部溫度越低，底部硬殼區(qū)域的溫度梯度較明顯。

圖5　 SIMECO 9號試驗數(shù)值模擬溫度場Fig．5 Temperature field of numerical simulation for SIMECO No．9experiment

圖6為SIMECO 9號試驗測量獲得的中心線溫度分布與CFD數(shù)值模擬結(jié)果的比較。從圖中可看出，由于氧化池底部硬殼熱阻較大，溫度上升較為迅速，進入熔融氧化池區(qū)域，溫度變化較小，即熔池內(nèi)熔融區(qū)域的溫度分布較為均勻。氧化池到金屬層交界的區(qū)域，溫度在硬殼和兩個流動邊界層的影響下，出現(xiàn)一階躍下降，因這個硬殼分界面和上下兩個邊界層起到了熱阻的作用。金屬層的熔融區(qū)域溫度分布梯度較小，到了金屬層的凝固硬殼區(qū)域，出現(xiàn)了較為明顯的溫度下降，這是由于導(dǎo)熱速率相對于對流換熱速率很低。

圖6　 SIMECO 9號試驗中心線溫度分布與CFD結(jié)果比較Fig．6 Temperature distribution along center line of SIMECO No．9experiment compared with CFD result

圖7為SIMECO 9號試驗測量獲得的熱流密度分布與CFD模擬結(jié)果的比較。可看出，目前CFD獲得的結(jié)果明顯高于試驗結(jié)果，說明了模擬獲得的熱流密度分布相對試驗結(jié)果存在一定誤差。在CFD結(jié)果中，95°位置附近熱流密度分布有一峰值，這是由金屬層底部等效為均勻熱流邊界和金屬層的高熱導(dǎo)率引起的。總體上，CFD技術(shù)在準確預(yù)測兩層熔池?zé)崃髅芏壬洗嬖谝欢ㄕ`差，但趨勢符合較好。

圖7　 SIMECO 9號試驗熱流密度分布與CFD結(jié)果比較Fig．7 Heat flux distribution of SIMECO No．9 experiment compared with CFD result

3　結(jié)論

本文將Enthalpy－porosity模型應(yīng)用于OpenFOAM平臺，分析了KTH開展的二維模擬材料氧化池和金屬層耦合傳熱試驗，研究了內(nèi)熱源驅(qū)動和溫差驅(qū)動湍流自然對流傳熱與凝固現(xiàn)象。計算獲得了液相熔融物和硬殼的分布，以及溫度和熱流密度分布，與試驗數(shù)據(jù)對比表明，該數(shù)值模型可初步適用于熔融池傳熱和相變研究，為研究IVR條件下堆芯熔融池的行為特性提供了參考。

參考文獻：

［1］ SEHGAL B R，THEERTHAN A，GIRI A，et al．Assessment of reactor vessel integrity（ARVI）［J］．Nucl Eng Des，2003，221：23－53．

［2］ ESMAILI H，KHATIB RAHBAR M，BASU S，et al．Analysis of in－vessel retention and exvessel fuel coolant interaction for AP1000［R］．Washington D．C．：U．S．NRC，2004．

［3］ KRETZSCHMAR F，F(xiàn)LUHRER B．Behavior of the melt pool in the lower plenum of the reactor pressure vessel－review of experimental programs and background of the LIVE program［R］．Germany：Karlsruhe Institute of Technology，2008．

［4］ CHI T T．The effective connectivity model for simulation and analysis of melt pool heat transfer in a light water reactor pressure vessel lower head［D］．Sweden：Royal Institute of Technology，2009．

［5］ REMPE J L．In－vessel retention strategy for high power reactors［R］．Idaho：Idaho National Engineering and Environmental Laboratory，2005．

［6］ VOLLER V R，PRAKASH C．A fixed－grid numerical modeling methodology for convection－diffusion mushy region phase－change problems［J］．Heat Mass Transfer，1987，30：1 709－1 720．

［7］ DARCY H．Les fontaines publiques de la ville de dijon［R］．Paris：Victor Dalmont，1856．

［8］ CARMAN P C．Fluid flow through granular beds［R］．London：Institution of Chemical Engineers，1937．

［9］ R?SLER F，BRüGGEMANN D．Shell－and－tube type latent heat thermal energy storage：Numerical analysis and comparison with experiments［J］．Heat and Mass Transfer，2011，47（8）：1 027－1 033．

［10］ASKAR A G．Natural convection heat transfer in two－fluid stratified pools with internal heat sources［D］．Sweden：Royal Institute of Technology，2002．

Numerical Investigation on Natural Convection and
Solidification of Molten Pool with OpenFOAM

WANG Xi，MENG Zhao－can，CHENG Xu
（State Nuclear Power Research Institute，Beijing100029，China）

Abstract：The in－vessel retention is adopted by the third generation nuclear power technology as an important severe accident mitigation strategy．The integrity of reactor pressure vessel depends on the heat flux distribution of molten pool．In present study，the solidification model in open source CFD software OpenFOAM was applied to simulate solidification and natural convection which was driven by internal heat source or temperature difference．The stratified molten pool heat transfer experiment carried out by Royal Institute of Technology was analyzed in the paper，and the solidified crust，temperature and heat flux distributions were obtained．The simulation results were compared with experimental data．It is shown that this numerical method can be used in the simulation of natural convection and solidification of molten pool，and it will probably be used in the analysis of molten corium behavior in reactor lower head．

Key words：severe accident；CFD；in－vessel retention；solidification；natural convection

作者簡介：王 溪（1985—），男，山東成武人，助理工程師，碩士，核能科學(xué)與工程專業(yè)

收稿日期：2014－04－22；修回日期：2014－06－09

doi：10．7538／yzk．2015．49．08．1393

文章編號：1000－6931（2015）08－1393－06

文獻標志碼：A

中圖分類號：TL33