矩形窄通道再淹沒現(xiàn)象實驗和數(shù)值研究

姬建業(yè) 許 巍 柴 翔 劉曉晶 曾 未

1（上海交通大學 上海 200240）

2（中國核動力研究設計院 成都 610041）

矩形窄通道再淹沒現(xiàn)象實驗和數(shù)值研究

姬建業(yè)1許 巍1柴 翔1劉曉晶1曾 未2

1（上海交通大學 上海 200240）

2（中國核動力研究設計院 成都 610041）

在反應堆發(fā)生大破口事故時，再淹沒階段可以有效地降低燃料元件溫度，防止堆芯熔毀。為了預測再淹沒過程中板狀燃料元件的換熱特性，進行了豎直矩形窄縫通道底部再淹沒過程的實驗研究。針對實驗工況，基于商用軟件CFX，通過耦合分析加熱板和流體的方法研究豎直矩形窄縫通道底部再淹沒過程。通過將數(shù)值模擬結果與實驗結果進行對比，評價了相關模型的適用性，并驗證了計算流體動力學(Computational Fluid Dynamics, CFD)方法在預測再淹沒過程的有效性?；隍炞C后計算模型，對壁面初始溫度、入口流速對再淹沒過程的影響進行了分析，獲得了相關初始條件對壁面溫度變化的影響規(guī)律。

兩流體模型，再淹沒，矩形通道，流固耦合

板式燃料元件與常規(guī)棒式燃料元件的相比具有更大的傳熱面積，其通道的流動沸騰的換熱系數(shù)與常規(guī)通道相比有很大提高[1]。在反應堆設計中，一些研究用的高通量反應堆通常使用板式燃料元件，這些反應堆的燃料組件會包含一系列的窄縫通道。冷卻劑流經(jīng)這些窄縫通道時將熱量帶走。窄縫通道由于幾何形狀與常規(guī)圓形通道不同，其流動和換熱特性與常規(guī)通道有很大差別[2-3]，尤其是氣泡演化特性和氣液兩相速度。當大破口事故發(fā)生時，注入堆芯的冷卻劑對堆芯進行再淹沒。再淹沒的過程中壁面溫度的變化會導致通道中出現(xiàn)不同的流動形態(tài)[4]。不同流動形態(tài)換熱特性之間的巨大差異，增加了數(shù)值模擬的復雜性。

模擬再淹沒過程通常使用系統(tǒng)程序如Relap5[5]、COBRA-TF等[6]。系統(tǒng)程序通常是用來進行安全分析和復雜系統(tǒng)設計，一維塊狀網(wǎng)格和經(jīng)驗關系式的使用導致流動和傳熱的細節(jié)被忽略。目前計算流體力學(Computational Fluid Dynamics, CFD)方法進行再淹沒的研究相對較少。CFD方法通過求解連續(xù)性、動量和能量方程，使用歐拉兩流體模型對溫度、空泡份額等參數(shù)進行求解。通過壁面換熱模型對流體和固體之間的換熱量進行求解。采用CFD方法計算的幾何形狀可以更復雜，局部的流動傳熱信息也更為詳盡。

本文首先通過再淹沒實驗得到固體的溫度，之后使用商用CFD軟件CFX進行模擬。通過再淹沒實驗得到了壁面溫度對驟冷前沿推進速度及再潤濕溫度的影響。通過將數(shù)值模擬結果與實驗值進行比對分析，研究了數(shù)值模擬對再淹沒過程預測的準確性及不同初始條件對該過程的影響。通過數(shù)值模擬得到了平均空泡份額和壁面熱流密度分布，以及壁面熱流密度最大時的空泡份額。該研究可以從機理上分析再淹沒過程，為以后CFD方法開發(fā)再淹沒計算模型奠定基礎。

1 矩形窄通道再淹沒實驗

1.1 矩形通道再淹沒實驗簡介

實驗臺架THERMAL是矩形窄通道再淹沒實驗臺架，實驗循環(huán)水溫度為(24±1) oC，流量為(24±1)L·h-1，壁面初始溫度為[(250-500)±1] oC。實驗回路在常壓下運行，研究不同初始條件對矩形窄通道中再淹沒過程的影響。

其中實驗回路系統(tǒng)圖如圖1所示。實驗系統(tǒng)由主回路系統(tǒng)和實驗回路系統(tǒng)組成。主回路主要用于調(diào)節(jié)包括進口流速和進口溫度等在內(nèi)的實驗初始條件。實驗系統(tǒng)主要包括實驗段、汽水分離裝置、夾帶水箱和蒸汽過熱器。首先使用陶瓷加熱器將實驗段加熱到指定溫度，實驗開始時關閉加熱器，然后控制氣動閥開關改變水的流動方向，使穩(wěn)定流量的水從實驗段底部向頂部流動。實驗回路中蒸汽發(fā)生器用來產(chǎn)生飽和蒸汽，在實驗開始前預熱實驗回路。實驗中再淹沒發(fā)生時，蒸汽會夾帶液體通過汽水分離器。液滴會從汽水分離器中流入夾帶水箱。用來觀測實驗中的夾帶量。壓力控制閥用于調(diào)節(jié)系統(tǒng)壓力，在高壓工況下需要使用該裝置。目前尚未進行夾帶量和高壓工況實驗。

圖1 再淹沒實驗系統(tǒng)回路Fig.1 Schematic diagram of reflooding test facility.

圖2 為實驗段，圖2中標注的尺寸是窄通道的尺寸。其中流道是由不銹鋼板、玻璃和石墨墊片所圍成的矩形通道。實驗段在實驗中由保溫材料包圍。不銹鋼板的厚度為10 mm。流道高度為600 mm，長度為60 mm，寬度為2 mm。

圖2 實驗段流道示意圖Fig.2 Overview of the test-section.

圖3為實驗段剖面圖，實驗中采用陶瓷加熱板將不銹鋼板加熱到指定溫度。采用夾持裝置來固定玻璃和不銹鋼板，以保持流道的密封性。

圖3 實驗段剖面Fig.3 View of cross-section of test-section.

不銹鋼溫度由熱電偶測得，在實驗段中有5排，每排3個共15個熱電偶。熱電偶高度方向上位置分別100 mm、250 mm、300 mm、350 mm、500 mm，圖4為熱電偶位置圖。熱電偶插入的深度為9 mm，即熱電偶的測量點距離加熱表面為1 mm。由于熱電偶測點位置與流固交界面距離較近，數(shù)據(jù)處理中并未考慮1 mm距離所造成的溫差。

圖4 熱電偶位置Fig.4 Locations of different thermocouples.

1.2 實驗工況

實驗中保持實驗段的入口溫度不變。通過調(diào)節(jié)壁面的初始溫度和入口流速，從而分析壁面溫度和入口流速對再淹沒現(xiàn)象的影響。

表1 實驗工況Table1 Experiment condition.

1.3 驟冷前沿推進速度

圖5為工況1中熱電偶5和熱電偶8的溫度變化，實驗中將壁面溫度急劇下降點TAR作為該時刻驟冷前沿位置，根據(jù)熱電偶5和熱電偶8的溫度變化計算驟冷前沿推進速度。

圖5 工況1中不同位置熱電偶溫度變化Fig.5 Temperature variance of different thermocouples in case 1.

圖6 是實驗工況中熱電偶5位置的壁面溫度變化。在工況1-3中，壁面溫度分別在31 s、55 s、110s時迅速下降。當壁面溫度高于再潤濕溫度時，此時流動狀態(tài)為反環(huán)狀膜態(tài)沸騰。之后壁面溫度迅速下降，流動狀態(tài)從過度沸騰轉變?yōu)楹藨B(tài)沸騰，最后轉變?yōu)閱蜗嗔鲃?。由圖6可以看出，隨著壁面初始溫度的增高，驟冷前沿推進速度逐漸降低。工況1-3驟冷前沿推進速度分別為8.06 mm·s-1、4.54 mm·s-1、2.27 mm·s-1。

圖7是工況4-6中熱電偶5溫度分布。從圖7中可以看出，隨著入口流速的增加，驟冷前沿推進速度逐漸增加。這是由于入口流速增加時，液滴的夾帶量增加，先驅冷卻過程換熱增強，從而使壁面溫度先降到驟冷溫度，之后再潤濕發(fā)生。

圖6 不同初始壁面溫度熱電偶溫度變化Fig.6 Temperature variance of thermocouples at different surface temperature.

圖7 不同入口流速熱電偶溫度變化Fig.7 Temperature variance of thermocouples at different inlet velocity.

從以上實驗數(shù)據(jù)可知，再淹沒過程中驟冷前沿推進速度與初始壁面溫度和入口流速有關。初始壁面溫度增高時，驟冷前沿推進速度減小。入口流速增大時，驟冷前沿推進速度增大。

1.4 再潤濕溫度

再潤濕溫度是指液體可以重新接觸干燥的固體表面時刻的溫度[7]。當表面溫度冷卻到兩相流動流型發(fā)生變化，從膜態(tài)沸騰到過渡沸騰或核態(tài)沸騰時，即再潤濕發(fā)生。在圖6中TAR表示工況3中的再潤濕溫度，即驟冷前溫度曲線的切線與溫度曲線斜率絕對值最大點的切線的交叉點。從圖6中可以看出，再潤濕溫度與壁面初始溫度有關，壁面初始溫度增高時，再潤濕溫度增高。而對比圖7中不同入口流速情形，當入口流速改變時，再潤濕溫度幾乎不變。

2 數(shù)學模型

2.1 歐拉兩流體模型

在CFD兩相計算中，歐拉兩流體的模型被廣泛的應用于CFD兩相計算。歐拉兩流體模型有以下優(yōu)點：可以得到全局的參數(shù)、適用的體積分數(shù)范圍比

較大、計算量相對較少。

質(zhì)量守恒方程：

式中：rα為體積份額；SMSα為質(zhì)量源項；Γαβ為從α到β相的質(zhì)量流密度。

動量守恒方程：

式中：SMα為動量源項；Mα為表面力對α相的作用力。

能量守恒方程：

式中：SEα為外部熱源；Qα為其它相通過相間對α相的傳熱。

2.2 相間作用力模型

由于考慮到再淹沒過程中流型的變化比較復雜，相間作用力模型采用Mixture模型。相間曳力的表達式為：

式中：Cd為無量綱的曳力系數(shù)。其中混合項密度為：

相間傳熱面積為：

式中：dαβ為定義的混合尺度。

2.3 相間能量傳遞模型

在氣相和液相的交界面，用一個總的換熱系數(shù)來表示相間界面的換熱系數(shù)hαβ。

相間的換熱表達式如下：

在Mixture模型中，換熱系數(shù)hαβ描述為混合導熱系數(shù)與混合尺度的函數(shù)：

式中：λαβ=rαλα+rβλβ。

2.4 壁面沸騰模型

在近壁面，采用了計算壁面沸騰最常用的美國倫斯勒理工學院模型(Rensselaer Polytechnic Institute, RPI)[8]。在RPI模型中，壁面的換熱量被認為是三個量的積分。其中當流體的溫度超過飽和溫度時，流體就會發(fā)生相變。

式中：Qtot為壁面總換熱量；Qc為單相對流換熱項；

Qq為驟冷換熱項；Qe為蒸發(fā)項。

2.4.1 單相對流項

式(10)中Qc為通過對流通過壁面?zhèn)鬟f的熱量，表達式如下：

式中：Abubble為被氣泡占據(jù)的壁面份額；Tw為壁面溫度；nwT+為近壁面流體的溫度。其中對流換熱系數(shù)hc采用Kader壁面溫度函數(shù)關系式計算得到。其計算關系式如下：

式中：ρL是液相密度；Cp,L是液相等壓比熱容；uτ是壁面剪切速度；nwT+由Kader的溫度壁面函數(shù)關系式求得，即：

2.4.2 驟冷項

驟冷項Qq的表達式如下：

2008年全球金融危機爆發(fā)，導致全球經(jīng)濟衰退，全球貿(mào)易環(huán)境急劇惡化，貿(mào)易保護主義勢頭明顯回升，加上人民幣匯率浮動、勞動力成本和原材料價格等各種因素，國內(nèi)外貿(mào)企業(yè)在殘酷的現(xiàn)實面前逐漸認識到自身轉型升級的必要性。

其中驟冷項換熱系數(shù)hq為：

式中：f為氣泡脫離頻率，表達式為[9-10]：

twait為氣泡等待時間：

氣泡影響面積為：

式中：F2是一個可以調(diào)節(jié)的影響因素，默認值為2。dw為氣泡脫離直徑，根據(jù)Kurul和Podowski提出將氣泡脫離直徑采用如下關系式：

式中：dmax=1.4 mm，dref=0.6 mm，Tref=45 K。

2.4.3 蒸發(fā)項

蒸發(fā)項的Qe表達式如下：

式中：n為氣化核心密度，其表達式為：

式中：nref=0.8×9.922×105m-2；ΔTsup=min[max (Twall-Tsat,0),ΔTmax]；ΔTref=10 K；p=1.8；ΔTmax是為了防止在初始收斂過程中發(fā)生溢出而設置的壁面過熱度最大值，默認值為25 K。

3 CFD二維計算模型

圖8為CFD計算幾何模型。

數(shù)值模擬中入口為流速入口，出口為壓力出口，其他邊界均為絕熱邊界。不銹鋼固體給定初始溫度進行瞬態(tài)模擬。由于實驗段長度較長，進口流速較小，且驟冷前沿的上部固體溫度變化較小。驟冷前沿上部固體對下部的導熱與再淹沒過程中固體與流體的對流換熱相比很小。為了提高計算效率，將CFD計算幾何模型高度簡化為0.2 m。根據(jù)網(wǎng)格敏感性分析，計算中采用20×5的網(wǎng)格。該網(wǎng)格計算得到的溫度變化最接近實驗值。圖9是使用0.6 m與0.2 m的幾何得到的壁面溫度分布對比。圖9中溫度點的位置高0.1 m中間熱電偶位置。后面的計算中也取該點作為參考點。從圖9中可以看出，采用 0.2 m的幾何與采用0.6 m的幾何計算得到的在0.1 m位置中間熱電偶測點的溫度分布幾乎沒有差別，為了提高計算速度，采用0.2 m的幾何進行后續(xù)計算。

圖8 CFD計算幾何模型Fig.8 CFD geometry.

圖9 0.2 m與0.6 m的幾何溫度變化Fig.9 Temperature difference between 0.2-m and 0.6-m geometry.

4 計算結果及分析

4.1 不同初始條件分析

圖10為實驗結果和數(shù)值模擬結果在指定點的溫度變化。通過圖10對比模擬結果和實驗結果可以發(fā)現(xiàn)，數(shù)值模擬對驟冷前沿推進速度和再潤濕溫度的估計與實際值趨勢一致，但細節(jié)上有一定差別。在工況1中，t=13 s時驟冷前沿已經(jīng)推進至熱電偶的位置，而在模擬中直到約23 s后壁面溫度才開始大幅的下降。在工況2結果中，t=25 s時驟冷前沿已經(jīng)推進至熱電偶的位置，而模擬中直到約32 s壁面溫度開始大幅下降。在工況3中，t=47s驟冷前沿推進至熱電偶位置，而數(shù)值模擬中到約53 s時壁面溫度開始迅速下降。

通過實驗與數(shù)值模擬對比可以發(fā)現(xiàn)，數(shù)值計算對再淹沒過程中壁面溫度預測的趨勢與實驗值一致。但在數(shù)值模擬中，壁面溫度的下降趨勢比實驗值快，這是由于數(shù)值模擬對壁面換熱的估計比真實值大。而壁面換熱估計比真實值大的原因可能是窄通道壁面氣泡生成對換熱的阻礙較大，而數(shù)值計算中難以模擬這一現(xiàn)象。在流動初期，數(shù)值模擬與實驗溫度較為接近，而在流動末期，數(shù)值模擬中壁面溫度下降較快，而實驗中壁面溫度下降緩慢。這是由于RPI模型需設置湍流流動，而實際流動由于入口流速較低，雷諾數(shù)較低，流動后期主要為層流流動。數(shù)值模擬對于再潤濕溫度的預測普遍偏低，這可能是由于數(shù)值模擬對壁面換熱估計過大。

圖10 不同初始壁溫指定點溫度比較Fig.10 Temperature comparison of the given point at different surface temperature.

通過圖11對比不同入口流速的數(shù)值模擬結果可以發(fā)現(xiàn)，入口流速對驟冷前沿推進速度存在影響。當入口流速較快時，驟冷時間較短，再潤濕溫度降低。

圖11 不同入口流速指定點溫度比較Fig.11 Temperature comparison of the given point at different inlet velocity.

4.2 空泡份額與壁面熱流密度

圖12是數(shù)值模擬實驗工況1中通道內(nèi)水平方向的平均空泡份額在豎直方向的瞬態(tài)分布曲線。通過對比該圖可以發(fā)現(xiàn)，在不同位置氣體空泡份額有很大變化。進口處空泡份額較小，而出口處空泡份額較大。由于瞬態(tài)計算，初始時刻計算的不穩(wěn)定導致t=1 s時刻曲線的波動。之后隨著時間的增長，該曲線幾乎平行向右移動。

通過圖13對比數(shù)值模擬實驗工況1中不同時刻壁面熱流密度在豎直方向上分布可以發(fā)現(xiàn)，在軸向不同位置壁面熱流密度在(0.4-6.5)×105W·m-2變化。隨著時間的增長，熱流密度的最高點隨著驟冷前沿的推進逐漸向流道頂部移動。

圖14是t=11 s時水平方向的平均空泡份額分布和壁面熱流密度分布，通過對比圖11可以發(fā)現(xiàn)壁面熱流密度最大處的空泡份額約為0.3。

圖12 空泡份額隨時間變化分布Fig.12 Volume fraction variance with time.

圖13 壁面熱流密度隨時間變化分布Fig.13 Heat flux density variance with time.

圖14 11 s時刻空泡份額與壁面熱流密度分布Fig.14 Volume fraction and heat flux density distributions in t=11 s.

5 結語

本文通過矩形窄通道再淹沒現(xiàn)象實驗和數(shù)值研究得到如下結論：CFD方法在預測再淹沒現(xiàn)象中得到的壁面溫度變化趨勢與實驗值一致。但數(shù)值模擬中對壁面換熱和再潤濕溫度的估計要小于實驗值。這可能是矩形窄通道中再淹沒過程中壁面氣泡生成對換熱的阻礙較大，而數(shù)值計算中難以模擬這一現(xiàn)象。實驗過程發(fā)現(xiàn)壁面溫度增高，驟冷前沿推進速度降低，再潤濕溫度升高。入口流速減小，驟冷前沿推進速度降低，再潤濕溫度不變。數(shù)值模擬可知流道空泡份額與壁面熱流密度分布曲線隨時間的增長，不斷向流道頂部移動，在壁面熱流密度最大時空泡份額約為0.3。

通過本研究得到了驟冷前沿推進速度與再潤濕溫度的數(shù)值。這些實驗值將應用于計算板式元件反應堆大破口事故情形下再淹沒時間及包殼溫度，為提高計算準確性打下基礎。

后續(xù)的研究將會嘗試改變RPI模型中一些子模型，從而更加準確的預測再淹沒過程。

1 陳沖, 高璞珍, 譚思超, 等. 豎直窄矩形通道內(nèi)環(huán)狀流的流動傳熱特性[J]. 化工學報, 2015, 66(2): 537-544.

CHEN Chong, GAO Puzhen, TAN Sichao, et al. Flow and heat transfer characteristics of annular flow in vertical rectangular narrow channel[J]. CIESC Journal, 2015, 66(2): 537-544.

2 董化平, 樊文遠, 郭赟. 板狀燃料組件流量分配CFD研究與優(yōu)化[J]. 核技術, 2016, 39(8): 080603. DOI: 10.11889/j.0253-3219.2016.hjs.39.080603.

DONG Huaping, FAN Wenyuan, GUO Yun. CFD investigation and optimization on flow distribution of plate-type assembly[J]. Nuclear Techniques. 2016, 39(8): 080603. DOI: 10.11889/j.0253-3219.2016.hjs.39.080603.

3 Qu X X, Zhang Y J, Ji H J, et al. Experimental research on velocity distribution in narrow slots of plane type reactor fuel[J]. Atomic Energy Science and Technology, 2003, 37(6): 519-522.

4 Lokanathan M, Hibiki T. Flow regime, void fraction and interfacial area transport and characteristics of co-current downward two-phase flow[J]. Nuclear Engineering and Design, 2016, 307(1): 39-63.

5 曾未, 余紅星, 孫玉發(fā), 等. 基于RELAP5的窄縫通道再淹沒模型適應性研究[J]. 核動力工程, 2013, 34(3): 50-57.

ZENG Wei, YU Hongxing, SUN Yufa, et al. Research on reflooding model of narrow channel based on RELAP5[J]. Nuclear Power Engineering, 2013, 34(3): 50-57.

6 Moon S K, Kim J, Kim K, et al. Reflood experiments in rod bundles with flow blockages due to clad ballooning[J]. Kerntechnik, 2016, 81(3): 251-256.

7 Debbarma A, Pandey K M. Influence on rewetting temperature and wetting delay during rewetting rod bundle by various radial jet models[J]. Kerntechnik, 2016, 81(1): 50-59.

8 Mohamedi B, Hanini S, Ararem A, et al. Simulation of nucleate boiling under ANSYS-FLUENT code by using RPI model coupling with artificial neural networks[J]. Nuclear Science and Techniques, 2015, 26(4): 040601. DOI: 10.13538/j.1001-8042/nst.26.040601.

9 Kon?ar B, Matkovi? M. Simulation of turbulent boiling flow in a vertical rectangular channel with one heated wall[J]. Nuclear Engineering and Design, 2012, 245: 131-139.

10 李權, 焦擁軍, 于俊崇. 豎直加熱圓管內(nèi)過冷沸騰及CHF數(shù)值模擬[J]. 核動力工程, 2015, 36(1): 168-172.

LI Quan, JIAO Yongjun, YU Junchong. Simulation of subcooled boiling and critical heat flux in uniformly heated tube[J]. Nuclear Power Engineering, 2015, 36(1): 168-172.

Experimental investigation and CFD simulation of reflooding
in narrow rectangle channel

JI Jianye1XU Wei1CHAI Xiang1LIU Xiaojing1ZENG Wei2
1(Shanghai Jiao Tong University, Shanghai 200240, China)
2(Nuclear Power Institute of China, Chengdu 610041, China)

Background:Reflooding phenomenon plays an important role in nuclear reactor accident. During reflooding phase, the water level will rise from the bottom to the top with quenching the hot fuel rods, which is the most important phase to guarantee the integrity of cladding. Purpose: The aim is to predict the heat transfer in reflooding phenomenon and to investigate the initial conditions effect on the temperature variance. Methods: The experiments and numerical method were carried out to obtain the temperature variance. The THERMAL facility was built and the computational fluid dynamics (CFD) software CFX was used during the study. Results: The experiment data shows that the wall initial temperature has a big influence on the quench velocity and rewet temperature. After the comparison with the numerical simulation, the numerical model was validated. Through the numerical simulation, the relationship between inlet velocity and wall temperature was received. Conclusion: It can be concluded that the CFD method can be used to predict the heat transfer in reflooding phenomenon. And some of the heat transfer models need to be improved to improve the accuracy.

Two fluid model, Reflooding, Rectangle channel, Fluid solid coupled

JI Jianye, male, born in 1991, graduated from Harbin Engineering University in 2014, master student, focusing on two phase flow and heat transfer

LIU Xiaojing, E-mail: xiaojingliu@sjtu.edu.cn

TL3

10.11889/j.0253-3219.2017.hjs.40.010601

No.11275178)資助

姬建業(yè)，男，1991年出生，2014畢業(yè)于哈爾濱工程大學，現(xiàn)為碩士研究生，研究領域為兩相流動與傳熱

劉曉晶，E-mail: xiaojingliu@sjtu.edu.cn

2016-11-01，

2016-11-28

Supported by National Natural Science Foundation of China (No.11275178)

Received date: 2016-11-01, accepted date: 2016-11-28