楊宇鵬，王成龍，張大林，蘇光輝，田文喜，秋穗正

(西安交通大學(xué) 核科學(xué)與技術(shù)學(xué)院，陜西 西安 710049)

蒸汽發(fā)生器是核反應(yīng)堆動(dòng)力系統(tǒng)的重要組成部分，是將反應(yīng)堆一次側(cè)的熱量向二次側(cè)傳遞的關(guān)鍵樞紐。相較于直管式蒸汽發(fā)生器，螺旋管式直流蒸汽發(fā)生器(HCOTSG)有一些突出優(yōu)勢(shì)[1]：1) HCOTSG的特殊結(jié)構(gòu)，使得其相比較直管式蒸汽發(fā)生器，在體積相同的情況換熱效率大幅提升，同樣熱負(fù)載相同的條件下，其將使得動(dòng)力系統(tǒng)布置更加緊湊，也更符合核能未來(lái)發(fā)展的小型化趨勢(shì)；2) 螺旋管側(cè)與殼側(cè)流體的湍流程度提升，起到強(qiáng)化換熱的效果；3) 螺旋管具有自補(bǔ)償特性，對(duì)熱膨脹具有良好的調(diào)節(jié)能力。但相比之下螺旋管式蒸汽發(fā)生器存在流動(dòng)壓降增加、設(shè)備復(fù)雜性增加等缺點(diǎn)。因此，針對(duì)HCOTSG的研究中，應(yīng)綜合其傳熱與流動(dòng)性能進(jìn)行全面分析。目前已有較多針對(duì)螺旋管內(nèi)流動(dòng)過(guò)程的實(shí)驗(yàn)研究基礎(chǔ)，對(duì)其內(nèi)部的流動(dòng)換熱特性以及特殊現(xiàn)象進(jìn)行了探究[2-5]，對(duì)于其殼側(cè)與管側(cè)的耦合流動(dòng)換熱過(guò)程通常利用商用軟件以及自編程序開(kāi)展數(shù)值模擬計(jì)算研究[6-9]。但針對(duì)液態(tài)金屬外掠螺旋管束的流動(dòng)以及液態(tài)金屬HCOTSG殼側(cè)與管側(cè)的耦合流動(dòng)傳熱過(guò)程的研究相對(duì)較少。

本文提出利用計(jì)算流體力學(xué)軟件(FLUENT)對(duì)液態(tài)金屬HCOTSG殼側(cè)流體、管側(cè)流體以及管壁固體區(qū)域耦合流動(dòng)傳熱過(guò)程進(jìn)行數(shù)值模擬的方法，對(duì)其流動(dòng)傳熱特性進(jìn)行分析。

1 液態(tài)金屬HCOTSG概述

HCOTSG是一種蒸汽發(fā)生器常用形式，其被廣泛應(yīng)用于船用壓水堆、小型模塊化反應(yīng)堆[10]，高溫氣冷堆以及鉛鉍快堆[11]等。鉛基材料(鉛、鉛鉍或鉛鋰合金等)具有中子性能優(yōu)良、熱工性能優(yōu)良、安全性好等優(yōu)點(diǎn)，采用其作為反應(yīng)堆冷卻劑，能使反應(yīng)堆物理特性和安全運(yùn)行具有顯著優(yōu)勢(shì)[12]。在液態(tài)金屬反應(yīng)堆中，HCOTSG也得到較為廣泛應(yīng)用。世界第1座液態(tài)金屬反應(yīng)堆——法國(guó)鳳凰堆使用750 MWe螺旋管式蒸汽發(fā)生器，日本的280 MWe鈉冷快堆(SFR)、徳國(guó)的鈉冷快堆(SNR-300)以及美國(guó)的(ALMR)均采用了HCOTSG。

本文研究對(duì)象為鉛鉍-水HCOTSG，其由分層布置的螺旋管束與同心圓柱套筒組成。圓柱套筒與螺旋管外壁之間的外掠管束區(qū)域?yàn)橐淮蝹?cè)流體域(即殼側(cè)區(qū)域)，螺旋管內(nèi)壁所圍區(qū)域?yàn)槎蝹?cè)流體域(即管側(cè)區(qū)域)，其示意圖如圖1a所示。螺旋管內(nèi)為過(guò)冷水到過(guò)熱蒸汽的兩相流蒸發(fā)過(guò)程，殼側(cè)為液態(tài)鉛鉍外掠螺旋管束冷卻過(guò)程。HCOTSG簡(jiǎn)圖如圖1b所示，螺旋管束分三層均勻布置于圓柱套筒中，螺旋方向分別為左旋、右旋、左旋。為保證換熱均勻，各層螺旋管長(zhǎng)度基本保持一致。各層螺旋管的幾何參數(shù)列于表1。對(duì)此幾何結(jié)構(gòu)的液態(tài)金屬HCOTSG內(nèi)部的流動(dòng)傳熱過(guò)程開(kāi)展了數(shù)值模擬計(jì)算，并對(duì)結(jié)果進(jìn)行了分析。

圖1 HCOTSG示意圖Fig.1 Schematic diagram of HCOTSG

表1 HCOTSG幾何參數(shù)Table 1 Gometrical parameter of HCOTSG

2 數(shù)值模擬方法

2.1 控制方程與計(jì)算方法

1) 控制方程

針對(duì)螺旋管內(nèi)的水蒸發(fā)過(guò)程，其涉及到氣液兩相流動(dòng)以及兩相間的傳熱與界面作用過(guò)程。本文選用適合于界面追蹤的VOF(volume of fluid)多相流模型對(duì)氣液兩相流過(guò)程進(jìn)行數(shù)值模擬，連續(xù)性方程、動(dòng)量方程和能量方程如下。

體積分?jǐn)?shù)方程(連續(xù)性方程)，對(duì)第q相，有：

(1)

(2)

式中：α為體積分?jǐn)?shù)；S為源項(xiàng)；ρ為密度；v為速度；t為時(shí)間。

動(dòng)量方程：

(3)

式中：p為壓力；F為源項(xiàng)；μ為黏性系數(shù)。

物性參數(shù)由各相體積分?jǐn)?shù)決定，如對(duì)于密度有：

ρ=∑αqρq

(4)

能量方程：

(5)

(6)

式中：keff為有效熱導(dǎo)率；當(dāng)i取1～3時(shí)，ui代表x、y和z方向的速度。

當(dāng)涉及到固體區(qū)域時(shí)，其能量方程為：

(7)

式中：ks為導(dǎo)熱系數(shù)；T為溫度。

2) 湍流模型

k-ε模型擁有可靠、收斂性好、內(nèi)存需求低等優(yōu)點(diǎn)。在FLUENT中的不同湍流模型中，Realizablek-ε模型對(duì)于二次流型是針對(duì)低雷諾數(shù)以及射流模擬提出的修正模型。考慮到螺旋管內(nèi)的流動(dòng)為湍流工況，且伴有二次流的發(fā)生，因此選擇此模型為數(shù)值模擬中采用的湍流模型?？紤]到本文的工況，選用的模型為Realizablek-ε湍流模型，其假設(shè)流動(dòng)的情況為完全湍流，從以下的輸運(yùn)方程中可獲得湍動(dòng)能k和湍流耗散率ε：

Gk+Gb-ρε-YM+Sk

(8)

(9)

(10)

式中：Gk為速度梯度引起的湍動(dòng)能增量；Gb為浮力引起的湍動(dòng)能增量；Ym為脈動(dòng)膨脹引起的湍流耗散率增量；C2和C1ε為常數(shù)；σk和σε分別為k、ε的湍流普朗特?cái)?shù)；Sk和Sε為源項(xiàng)。

2.2 液態(tài)金屬流動(dòng)傳熱模型

相比于常規(guī)流體，液態(tài)金屬的熱物理性質(zhì)存在較大差異，如液態(tài)金屬的導(dǎo)熱能力遠(yuǎn)大于常規(guī)流體。因此，對(duì)于液態(tài)金屬，其熱擴(kuò)散能力遠(yuǎn)大于動(dòng)量擴(kuò)散能力，導(dǎo)致其普朗特?cái)?shù)(Pr≈0.01～0.06)與常規(guī)流體(Pr≈1)相差巨大。因此，對(duì)于液態(tài)金屬流動(dòng)的數(shù)值模擬，湍流普朗特?cái)?shù)(Prt)的影響較為顯著。因此在采用CFD方法對(duì)液態(tài)金屬流動(dòng)傳熱進(jìn)行數(shù)值模擬時(shí)通常需對(duì)Prt模型進(jìn)行修正。本文選擇了Cheng等[13]提出的液態(tài)金屬Prt模型，利用用戶自定義文件(UDF)導(dǎo)入FLUENT軟件對(duì)湍流模型進(jìn)行修正，公式如下：

(11)

式中：Pe為佩克萊數(shù)；A為常數(shù)。

2.3 數(shù)值模擬方法驗(yàn)證

由于對(duì)液態(tài)金屬HCOTSG的實(shí)驗(yàn)數(shù)據(jù)相對(duì)較少，因此對(duì)其殼側(cè)液態(tài)金屬流動(dòng)、管側(cè)水蒸發(fā)兩相流動(dòng)分別進(jìn)行模型驗(yàn)證，以證明數(shù)值模擬方法的正確性。

1) 螺旋管內(nèi)水蒸發(fā)模擬驗(yàn)證

對(duì)于管內(nèi)的水沸騰兩相流過(guò)程，首先選擇Bartolomei等[14]于1967年設(shè)計(jì)搭建的豎直加熱圓管內(nèi)流動(dòng)沸騰實(shí)驗(yàn)進(jìn)行了對(duì)比驗(yàn)證，圓管長(zhǎng)度為2 m，直徑為15.4 mm，壁面邊界為恒定熱流邊界，所得數(shù)值模擬結(jié)果如圖2所示。得到的沿軸向高度的空泡份額與實(shí)驗(yàn)值誤差在25%以內(nèi)，溫度分布結(jié)果與實(shí)驗(yàn)值誤差在3%以內(nèi)，均符合較好。

在此基礎(chǔ)上，將模擬方法應(yīng)用于螺旋管內(nèi)水沸騰過(guò)程模擬。選取Santini等[15]所做實(shí)驗(yàn)的結(jié)果與數(shù)值模擬結(jié)果進(jìn)行對(duì)比。螺旋管內(nèi)徑為12.53 mm，螺旋直徑為1 000 mm，節(jié)距為800 mm，長(zhǎng)為32 mm的螺旋管豎直布置。所得數(shù)值模擬結(jié)果如圖3所示。由不同工況下的換熱系數(shù)對(duì)比可看出，所得結(jié)果與實(shí)驗(yàn)值符合較好，大部分工況誤差在20%以內(nèi)，因此螺旋管內(nèi)流動(dòng)的數(shù)值模擬方法正確性得到驗(yàn)證。

圖2 Bartolomei實(shí)驗(yàn)值與計(jì)算值對(duì)比Fig.2 Comparison between experimental result of Bartolomei and calculated result

圖3 Santini實(shí)驗(yàn)值和計(jì)算值對(duì)比Fig.3 Comparison between experimental result of Santini and calculated result

2) 螺旋管外液態(tài)鉛鉍流動(dòng)模擬驗(yàn)證

殼側(cè)的鉛鉍流動(dòng)過(guò)程為液態(tài)金屬外掠螺旋管束流動(dòng)，缺乏相關(guān)的實(shí)驗(yàn)數(shù)據(jù)。因此，首先選擇西安交通大學(xué)圓管內(nèi)的鉛鉍流動(dòng)換熱實(shí)驗(yàn)數(shù)據(jù)進(jìn)了對(duì)比驗(yàn)證，圓管直徑為20 mm，長(zhǎng)度為2 m。圖4為西安交通大學(xué)(XJTU)圓管內(nèi)鉛鉍流動(dòng)實(shí)驗(yàn)值和計(jì)算值對(duì)比。由圖4可知，在試驗(yàn)工況下，圓管平均Nu的數(shù)值模擬結(jié)果與實(shí)驗(yàn)值符合較好，大部分工況誤差在15%以內(nèi)。驗(yàn)證了鉛鉍流動(dòng)數(shù)值模擬方法的準(zhǔn)確性。

圖4 XJTU圓管內(nèi)鉛鉍流動(dòng)實(shí)驗(yàn)值和計(jì)算值對(duì)比Fig.4 Comparison between experimental result of XJTU and calculated result

在驗(yàn)證液態(tài)鉛鉍模擬方法的基礎(chǔ)上，考慮到螺旋管式蒸汽發(fā)生器殼側(cè)區(qū)域特殊的幾何結(jié)構(gòu)，進(jìn)一步建立更為相近的液態(tài)金屬橫掠管束模型，并利用已有的經(jīng)驗(yàn)關(guān)系式對(duì)數(shù)值模擬方法進(jìn)行驗(yàn)證。在建立了液態(tài)金屬鈉橫掠管束模型后對(duì)其開(kāi)展數(shù)值計(jì)算，選擇液態(tài)金屬鈉橫掠管束的Kalish-Dwyer關(guān)系式以及Schad關(guān)系式與數(shù)值模擬結(jié)果進(jìn)行驗(yàn)證，結(jié)果如圖5所示，可看出，數(shù)值模擬結(jié)果與關(guān)系式計(jì)算結(jié)果符合較好，大部分工況誤差在15%以內(nèi)，進(jìn)一步驗(yàn)證了殼側(cè)液態(tài)金屬流體域數(shù)值模擬方法的準(zhǔn)確性。

圖5 關(guān)系式計(jì)算值和模擬值對(duì)比Fig.5 Comparison between calculated result by empirical correlation and simultion result

2.4 模型建立與網(wǎng)格無(wú)關(guān)性分析

1)模型建立

依據(jù)表1中的幾何參數(shù)，利用SOLIDWORKS軟件對(duì)螺旋管束幾何結(jié)構(gòu)進(jìn)行建模。螺旋管共布置3層，共9根，殼側(cè)區(qū)域?yàn)樘淄才c螺旋管束外壁形成的外掠管束區(qū)域。模型包含殼側(cè)流體域、管側(cè)流體域以及管壁固體域三部分。在幾何模型建立完成后進(jìn)行網(wǎng)格劃分，管側(cè)流體與管壁區(qū)域采用六面體網(wǎng)格，殼側(cè)流體域采用四面體網(wǎng)格，部分區(qū)域的網(wǎng)格示意圖如圖6所示。網(wǎng)格劃分完成后將建立模型導(dǎo)入FLUENT軟件進(jìn)行數(shù)值模擬計(jì)算。在進(jìn)行數(shù)值模擬時(shí)，需將液態(tài)鉛鉍的物性以及選用的Prt模型利用UDF導(dǎo)入軟件并進(jìn)行相應(yīng)參數(shù)的選取。入口邊界采用速度入口邊界條件，出口邊界采用壓力出口邊界條件，內(nèi)外套筒壁面均為絕熱邊界條件。綜合考慮管側(cè)與殼側(cè)流動(dòng)過(guò)程，湍流模型選擇Realizablek-ε模型。在監(jiān)控計(jì)算殘差的同時(shí)對(duì)出口溫度、含汽率等參數(shù)進(jìn)行監(jiān)控以保證計(jì)算結(jié)果的準(zhǔn)確性。

圖6 部分區(qū)域網(wǎng)格示意圖Fig.6 Schematic diagram of partial area grid

2) 網(wǎng)格無(wú)關(guān)性分析

由于本文的數(shù)值模擬包含HCOTSG殼側(cè)液態(tài)鉛鉍流動(dòng)與管側(cè)水蒸發(fā)兩相流動(dòng)過(guò)程，兩側(cè)流動(dòng)傳熱過(guò)程的差異較大，網(wǎng)格劃分方式對(duì)數(shù)值模擬結(jié)果產(chǎn)生的影響也差異較大，因此對(duì)兩側(cè)的流體流動(dòng)傳熱過(guò)程分別進(jìn)行網(wǎng)格無(wú)關(guān)性分析。對(duì)于螺旋管內(nèi)的兩相流過(guò)程，本文選取了沿軸向高度不同位置的截面平均壓降作為網(wǎng)格無(wú)關(guān)性的判定依據(jù)，所得結(jié)果如圖7所示。對(duì)螺旋管周向節(jié)點(diǎn)數(shù)量與軸向節(jié)點(diǎn)數(shù)量分別進(jìn)行加密處理，發(fā)現(xiàn)當(dāng)周向網(wǎng)格數(shù)量為30、軸向網(wǎng)格數(shù)量為700時(shí)，加密網(wǎng)格對(duì)軸向沿程壓降的影響較小。因此選擇此時(shí)的網(wǎng)格劃分方式對(duì)螺旋管內(nèi)流體域及管壁固體域進(jìn)行網(wǎng)格劃分。

圖7 管側(cè)網(wǎng)格無(wú)關(guān)性分析Fig.7 Grid independence analysis at tube side

對(duì)于HCOTSG殼側(cè)的液態(tài)鉛鉍流動(dòng)過(guò)程，本文選取沿軸向高度的截面平均溫度作為網(wǎng)格無(wú)關(guān)性分析的判定依據(jù)，所得結(jié)果如圖8所示。首先對(duì)網(wǎng)格進(jìn)行加密，發(fā)現(xiàn)當(dāng)網(wǎng)格數(shù)量達(dá)到580萬(wàn)左右時(shí)，繼續(xù)加密網(wǎng)格對(duì)殼側(cè)液態(tài)鉛鉍的沿程溫度影響較小。因此選擇此時(shí)的全局尺寸劃分網(wǎng)格。在確定全局網(wǎng)格尺寸的基礎(chǔ)上，選取了不同的邊界層網(wǎng)格劃分方式進(jìn)行網(wǎng)格劃分，當(dāng)壁面邊界y+縮小至34時(shí)，繼續(xù)加密邊界層網(wǎng)格軸向溫度的影響較小。因此選擇邊界層網(wǎng)格處理方式對(duì)邊界層網(wǎng)格進(jìn)行劃分。

3 流動(dòng)換熱特性分析

在方法正確性驗(yàn)證的基礎(chǔ)上，選取典型工況的計(jì)算結(jié)果進(jìn)行分析，工況參數(shù)列于表2。圖9示出殼側(cè)區(qū)域中截面速度分布。由圖9可看出，在不同層螺旋管之間存在明顯的流體加速區(qū)域，而同層螺旋管的軸向節(jié)距之間存在明顯的低速區(qū)域，這使得液態(tài)金屬域內(nèi)的流體攪混明顯提升，強(qiáng)化了流體與管壁的傳熱。且由于各層螺旋管螺旋方向布置方式不同，使得在截面對(duì)稱區(qū)域的速度分布并非對(duì)稱，因此在后續(xù)的敏感性分析中，螺旋方向的布置情況也應(yīng)作為考慮的敏感性因素。

圖8 殼側(cè)網(wǎng)格無(wú)關(guān)性分析Fig.8 Grid independence analysis at shell side

表2 工況參數(shù)Table 2 Working condition

圖9 殼側(cè)區(qū)域中截面速度分布Fig.9 Velocity distribution of middle section in shell side area

在殼側(cè)液態(tài)鉛鉍域內(nèi)，其沿高度方向的截面平均溫度與壓降分布如圖10所示。進(jìn)出口的直管段區(qū)域，管束對(duì)殼側(cè)流體擾動(dòng)較小，因此換熱器內(nèi)流體的換熱效果明顯弱于螺旋管束區(qū)域。相應(yīng)地，螺旋管束區(qū)域的壓降也明顯高于直管區(qū)域。在利用HCOTSG的優(yōu)點(diǎn)時(shí)，應(yīng)同時(shí)關(guān)注其特殊結(jié)構(gòu)帶來(lái)的壓降損失。

對(duì)于管內(nèi)的水蒸發(fā)兩相流動(dòng)，各層內(nèi)的螺旋管沿程溫度與熱流密度分布如圖11所示?？煽闯?，各層螺旋管在同一軸向高度下的溫度與壁面熱流密度分布一致，因此在此幾何布置方式下，螺旋管束傳熱均勻性較好。在入口段與蒸發(fā)段，較大的溫差以及水蒸發(fā)的汽化潛熱使得壁面的熱流密度維持較高的水平。而在水完全蒸發(fā)后，螺旋管內(nèi)的換熱過(guò)程為水蒸氣的過(guò)熱過(guò)程，此區(qū)域內(nèi)熱流密度處于較低水平。因此，為更好發(fā)揮蒸汽發(fā)生器的換熱能力，針對(duì)此結(jié)構(gòu)的HCOTSG，應(yīng)在分析工況參數(shù)影響的同時(shí)選擇出合適的工況范圍。

圖10 殼側(cè)區(qū)域沿程溫度和壓降分布Fig.10 Distribution of temperature and pressure drop along shell side area

圖11 管側(cè)溫度與熱流密度分布Fig.11 Distribution of temperature and wall heat flux along tube side

圖12 螺旋管內(nèi)溫度與壓降沿程分布Fig.12 Temperature and pressure drop distribution in helical coil tube

螺旋管內(nèi)溫度與壓降沿程分布如圖12所示，由于入口段和出口段與接管處存在較大的彎曲角度，因此在這兩處壓降損失較大，螺旋區(qū)域內(nèi)壓降變化較為均勻。沿程溫度在入口段快速升高，上升趨勢(shì)逐漸平緩。螺旋管內(nèi)部流體除了受到黏性力作用，還受到慣性力作用。最內(nèi)層螺旋管內(nèi)的沿程截面液相分布變化如圖13所示。在螺旋管內(nèi)，液態(tài)水所受離心力較大，因此其位于螺旋管外側(cè)；而在部分區(qū)域，在螺旋管內(nèi)會(huì)出現(xiàn)特殊的液膜倒置現(xiàn)象，即液態(tài)水集中在螺旋管內(nèi)側(cè)的現(xiàn)象，該現(xiàn)象與已有的螺旋管內(nèi)兩相流動(dòng)實(shí)驗(yàn)情況一致。

圖13 管側(cè)區(qū)域截面液相體積份額Fig.13 Liquid volume fraction in tube side area

4 總結(jié)

本文提出了利用計(jì)算流體力學(xué)對(duì)液態(tài)金屬HCOTSG殼側(cè)流體、管側(cè)流體耦合流動(dòng)傳熱過(guò)程進(jìn)行數(shù)值模擬的方法，分析了其流動(dòng)傳熱特性，得出以下結(jié)論。

1) 本文提出的數(shù)值模擬方法可對(duì)液態(tài)金屬螺旋管式蒸汽發(fā)生器內(nèi)部的流動(dòng)換熱情況進(jìn)行模擬，并將所得計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果和經(jīng)驗(yàn)關(guān)系式結(jié)果對(duì)比：對(duì)于管內(nèi)的水蒸發(fā)兩相流模擬，通過(guò)與Bartolomei的直管沸騰實(shí)驗(yàn)以及Santini的螺旋管內(nèi)沸騰實(shí)驗(yàn)進(jìn)行對(duì)比，誤差均在25%以內(nèi)；對(duì)于殼側(cè)液態(tài)金屬模擬，通過(guò)與西安交通大學(xué)圓管鉛鉍流動(dòng)換熱實(shí)驗(yàn)以及Kalish-Dwyer關(guān)系式和Schad關(guān)系式進(jìn)行對(duì)比，誤差均在15%以內(nèi)，數(shù)值模擬方法的正確性得到驗(yàn)證。

2) 在本文提出的螺旋管式直流蒸汽發(fā)生器幾何結(jié)構(gòu)下，針對(duì)典型工況的流動(dòng)換熱情況進(jìn)行了分析。殼側(cè)鉛鉍域存在明顯的加速區(qū)域，且由于各層螺旋管的布置方式存在差異，因此在各層間的流體速度也存在差異，增加了流體域內(nèi)的攪混。對(duì)于螺旋管內(nèi)部的兩相流動(dòng)過(guò)程，過(guò)冷水被加熱最終變?yōu)檫^(guò)熱蒸汽，在過(guò)熱蒸汽區(qū)的熱流密度明顯下降，因此應(yīng)針對(duì)不同幾何結(jié)構(gòu)，確定合適的工況范圍。同時(shí)發(fā)現(xiàn)管內(nèi)存在液膜倒置現(xiàn)象，與已有實(shí)驗(yàn)所得結(jié)果一致。