白一凡,呂海財(cái),王 漢

(華北電力大學(xué) 核科學(xué)與工程學(xué)院,北京 102206)

超臨界二氧化碳(SCO2)是能源動(dòng)力系統(tǒng)中最有前景的傳熱和能量轉(zhuǎn)換工質(zhì)之一,在工業(yè)系統(tǒng)中得到了越來越廣泛的應(yīng)用。在核能發(fā)電領(lǐng)域,SCO2氣冷快堆是國內(nèi)外的研究熱點(diǎn),這一系統(tǒng)采用SCO2作為一回路冷卻劑吸收燃料元件的裂變熱,流出反應(yīng)堆壓力容器后直接進(jìn)入汽輪機(jī)做功[1-2]。相較于傳統(tǒng)壓水堆,SCO2氣冷快堆體積小、系統(tǒng)簡單、熱效率高,同時(shí)易于實(shí)現(xiàn)模塊化,可大幅降低核電廠建造成本。

堆芯燃料組件內(nèi)SCO2的流動(dòng)與傳熱行為是反應(yīng)堆熱工水力設(shè)計(jì)中十分關(guān)注的問題。燃料組件一般由若干燃料棒按一定的設(shè)計(jì)方式排列,形成大規(guī)模棒束結(jié)構(gòu)。但是在實(shí)際研究中,特別是在實(shí)驗(yàn)研究方面,研究人員經(jīng)常針對圓管、方管、環(huán)形通道等簡單通道開展超臨界流體的傳熱特性實(shí)驗(yàn)測量,獲得壁面溫度和傳熱系數(shù)的變化規(guī)律,驗(yàn)證或建立傳熱關(guān)系式。Mori等[3]對超臨界壓力下氟利昂在圓管、3棒束和7棒束內(nèi)的傳熱特性進(jìn)行了實(shí)驗(yàn)研究,發(fā)現(xiàn)在相同的實(shí)驗(yàn)條件下,圓管出現(xiàn)了傳熱惡化現(xiàn)象,而在棒束通道中則沒有發(fā)生。Kim等[4]對熱力當(dāng)量直徑為4.5 mm的環(huán)形通道內(nèi)SCO2的傳熱特性進(jìn)行了實(shí)驗(yàn)測量,并將實(shí)驗(yàn)結(jié)果與內(nèi)徑為4.4 mm的圓管進(jìn)行了對比,發(fā)現(xiàn)當(dāng)發(fā)生傳熱惡化時(shí),環(huán)形通道的傳熱系數(shù)略微高于圓管。Bae[5]對圓管和環(huán)形通道內(nèi)SCO2的傳熱系數(shù)進(jìn)行了實(shí)驗(yàn)測量,圓管內(nèi)徑為4.57 mm,環(huán)形通道加熱管直徑為8 mm,間隙為2 mm。對比結(jié)果表明,在相同的實(shí)驗(yàn)條件下,環(huán)形通道相比于圓管對傳熱惡化有一定的抑制作用。Li等[6]在上海交通大學(xué)超臨界水多用途實(shí)驗(yàn)回路上對圓管、環(huán)形通道和2×2棒束通道進(jìn)行了超臨界水的傳熱實(shí)驗(yàn)研究,發(fā)現(xiàn)2×2棒束的傳熱性能優(yōu)于圓管和環(huán)形通道。Liu等[7]使用低雷諾數(shù)湍流模型對圓管和環(huán)形通道中超臨界水的傳熱惡化現(xiàn)象進(jìn)行了數(shù)值模擬研究,發(fā)現(xiàn)在高質(zhì)量流速條件下,圓管的傳熱惡化現(xiàn)象比內(nèi)壁加熱的環(huán)形通道嚴(yán)重,但是當(dāng)環(huán)形通道的內(nèi)外壁面同時(shí)加熱時(shí),其傳熱惡化壁面溫度升高的幅度和變化規(guī)律與圓管相似。研究還認(rèn)為,相比于水力當(dāng)量直徑,熱力當(dāng)量直徑對高質(zhì)量流速下傳熱惡化的發(fā)生起重要作用。

目前國內(nèi)外對超臨界流體異常傳熱現(xiàn)象的實(shí)驗(yàn)或數(shù)值模擬研究以圓管和環(huán)形通道為主,但是由于幾何結(jié)構(gòu)的差異,簡單通道難以模擬燃料組件內(nèi)冷卻劑的橫向交混以及燃料元件周向傳熱不均勻等流動(dòng)和傳熱特點(diǎn),現(xiàn)階段缺乏由流道幾何結(jié)構(gòu)不同所導(dǎo)致的傳熱差異的研究。本文采用大渦模擬(LES)方法對SCO2在相同熱力當(dāng)量直徑、相同工況條件下的圓管、環(huán)形通道和子通道中的流動(dòng)傳熱特性進(jìn)行計(jì)算研究,分析流道結(jié)構(gòu)對壁面溫度、平均速度分布、湍流統(tǒng)計(jì)量等結(jié)果的影響,并根據(jù)LES結(jié)果對典型SCO2傳熱關(guān)系式進(jìn)行評估。典型的SCO2直接循環(huán)反應(yīng)堆工作壓力約為20 MPa,堆芯進(jìn)口冷卻劑溫度為400 ℃。在設(shè)計(jì)工況下,冷卻劑的熱物性變化較為平緩,不容易發(fā)生傳熱惡化。但是考慮到實(shí)際運(yùn)行中系統(tǒng)參數(shù)的瞬態(tài)變化,如果發(fā)生破口事故,系統(tǒng)參數(shù)降至CO2臨界點(diǎn)附近時(shí),傳熱特性與設(shè)計(jì)工況有所不同,可能會(huì)發(fā)生傳熱惡化。因此,本文將計(jì)算工況設(shè)置在臨界點(diǎn)附近來分析此時(shí)CO2的傳熱特性和不同流道對傳熱的影響。

1 數(shù)值模擬設(shè)置

1.1 計(jì)算參數(shù)選取依據(jù)

首先,不同幾何流道的SCO2入口溫度T0相同,保證入口的熱物性參數(shù)相同;其次,入口位置的流動(dòng)應(yīng)滿足相似性準(zhǔn)則的要求,即3個(gè)幾何流道的雷諾數(shù)相等:

(1)

式中:Re0為入口雷諾數(shù);ρ0為入口流體密度,kg/m3;U0為入口平均流速,m/s;μ0為入口動(dòng)力黏度,Pa·s;D為熱力當(dāng)量直徑,m,其定義如下:

(2)

式中:A為通道橫截面積,m2;P為通道橫截面加熱壁面周長,m。

此外,不同幾何形狀流道的格拉曉夫數(shù)與無量綱熱流密度對應(yīng)相等,這兩個(gè)無量綱數(shù)的定義如下:

(3)

式中:Gr0為入口格拉曉夫數(shù);g為重力加速度,m/s2;β0為入口體積膨脹系數(shù),1/K;λ0為入口導(dǎo)熱系數(shù),W/(m·K);ν0為入口運(yùn)動(dòng)黏度,m2/s;qw為加熱壁面熱流密度,W/m2。

(4)

式中:q+為無量綱熱流密度;cp0為入口定壓比熱容,J/(kg·K)。

根據(jù)以上分析,不同幾何流道的U0和qw應(yīng)選擇相同的參數(shù),同時(shí)選取適當(dāng)幾何尺寸以保證D相同,這樣能夠使數(shù)值模擬結(jié)果的差異只來源于通道幾何結(jié)構(gòu)。

1.2 控制方程

本文使用LES方法研究SCO2的湍流傳熱特性,控制方程采用低馬赫數(shù)不可壓縮假設(shè),忽略熱物性隨壓力的變化,認(rèn)為SCO2的物性只與溫度有關(guān)。控制方程包括質(zhì)量守恒方程、動(dòng)量守恒方程和能量守恒方程,表示如下:

(5)

(6)

(7)

式中:下標(biāo)i、j為張量的階數(shù);ρ為密度;ui、uj為速度矢量;t為時(shí)間;xi、xj為位移矢量;μ為動(dòng)力黏度;p為壓力;gi為重力加速度矢量;μSGS為亞格子尺度(SGS)黏度,需要用SGS模型進(jìn)行計(jì)算;h為比焓;α為熱擴(kuò)散系數(shù);αSGS為SGS熱擴(kuò)散系數(shù),可由SGS黏度確定:

(8)

式中,PrSGS為SGS普朗特?cái)?shù),計(jì)算中采用OpenFOAM程序的默認(rèn)常數(shù)1.0。

對于LES計(jì)算,最常用的SGS模型是Smagorinsky模型[8],SGS黏度的定義如下:

(9)

(10)

本文還選取了壁面自適應(yīng)局部渦粘(WALE)模型[9]進(jìn)行SCO2的傳熱計(jì)算,該模型的選擇是基于Brockmeyer等[10]以及Busco等[11]的研究。WALE模型的黏度定義如下:

(11)

(12)

(13)

(14)

(15)

1.3 幾何結(jié)構(gòu)和邊界條件

圖1所示為LES選取的3種幾何結(jié)構(gòu),分別為圓管、環(huán)形通道和棒束子通道。圓管直徑d為5.276 mm;環(huán)形通道內(nèi)管外徑為5.276 mm,外管內(nèi)徑為7.461 4 mm;子通道呈四邊形布置,加熱棒直徑為9.5 mm,柵距p為10.5 mm。圓管壁面和子通道棒束壁面均為加熱壁面,而環(huán)形通道只有內(nèi)管外表面加熱,外管內(nèi)壁面絕熱。根據(jù)上述幾何結(jié)構(gòu)尺寸的選擇,3種通道的D均為5.276 mm。3種通道的流體質(zhì)量流速相等,D相等,q+相等。圓管與環(huán)形通道質(zhì)量流量相等,但略低于子通道的質(zhì)量流量。

a——圓管;b——環(huán)形通道;c——子通道

以圓管為例,邊界條件的設(shè)置如圖2所示。流動(dòng)方向設(shè)為x方向,橫向方向設(shè)為y和z方向。3種通道沿流動(dòng)方向的加熱段長度L均為45D。加熱段壁面為恒定熱流密度邊界條件,壁面處速度設(shè)置為無滑移邊界條件。加熱段出口設(shè)置為對流邊界條件,即?φ/?t+Ub·?φ/?x=0,其中:φ為標(biāo)量或矢量;x為位移矢量;Ub為流體平均速度。加熱段入口前設(shè)置長度L為6D的湍流入口發(fā)生器,其壁面絕熱,出入口設(shè)置為周期性邊界條件,以生成充分發(fā)展的湍流。在湍流入口發(fā)生器橫截面采集足夠長時(shí)間的速度分量,形成入口截面的湍流數(shù)據(jù)庫,然后按設(shè)定的時(shí)間步長將其添加到加熱段的入口,模擬真實(shí)且充分發(fā)展的入口湍流。在子通道中,4個(gè)棒間隙被分為兩對,一對在y方向,另一對在z方向,如圖1c所示。每對間隙采用周期性邊界條件,允許流體通過間隙自由進(jìn)出子通道區(qū)域。

1.4 計(jì)算工況

3種通道的LES計(jì)算工況均相同,即壓力p0=8.0 MPa,U0=0.204 m/s,Reb,0=13 028,T0=301.15 K,qw=36 kW/m2。3種幾何通道的流動(dòng)均為混合對流工況,g的方向與流動(dòng)方向相反,大小為9.8 m/s2。

1.5 網(wǎng)格劃分

圓管、環(huán)形通道與子通道橫截面的網(wǎng)格劃分如圖3所示,圖中數(shù)字表示箭頭所示方向的網(wǎng)格劃分?jǐn)?shù)量。圓管徑向劃分87個(gè)網(wǎng)格,周向劃分136個(gè)網(wǎng)格,軸向劃分750個(gè)網(wǎng)格,網(wǎng)格數(shù)為8 007 000。環(huán)形通道徑向劃分70個(gè)網(wǎng)格,周向劃分180個(gè)網(wǎng)格,軸向劃分645個(gè)網(wǎng)格,網(wǎng)格數(shù)為8 127 000。子通道徑向劃分75個(gè)網(wǎng)格,棒束周向劃分376個(gè)網(wǎng)格,軸向劃分720個(gè)網(wǎng)格,網(wǎng)格數(shù)為1 560萬。網(wǎng)格分辨率根據(jù)Pandey等[12]的計(jì)算方法,代入本文LES工況,計(jì)算得到圓管徑向、周向和軸向網(wǎng)格分辨率分別為0.24～4.8、18和48;環(huán)形通道徑向、周向和軸向網(wǎng)格分辨率分別為0.18～3.6、14～20和55;子通道徑向、周向和軸向網(wǎng)格分辨率分別為0.225～4.5、12和50。根據(jù)Choi等[13]的研究結(jié)論,網(wǎng)格分辨率滿足LES的要求。

a——圓管;b——環(huán)形通道;c——子通道

1.6 計(jì)算方法

本文使用開源程序OpenFOAM求解LES的質(zhì)量、動(dòng)量和能量守恒方程。SCO2的熱物性參數(shù)取自NIST數(shù)據(jù)庫[14],制作物性隨溫度變化的一維查詢表并植入到OpenFOAM中。數(shù)值模擬設(shè)置方面,使用PISO算法求解速度與壓力的耦合,壓力場采用GAMG求解器,殘差為1×10-7;速度場采用smoothSolver求解器,殘差為1×10-6。動(dòng)量和能量方程中的對流項(xiàng)采用二階中心差分格式離散,時(shí)間項(xiàng)采用二階向后隱式格式離散。無量綱時(shí)間步長Δt+=4.7×10-5,其定義為Δt+=Δt/t+,其中無量綱時(shí)間t+=D/uτ,并控制時(shí)間步長使計(jì)算過程中最大庫朗數(shù)小于0.4。當(dāng)入口的湍流發(fā)展到出口后開始時(shí)間平均,采樣時(shí)間約為7倍的流體流通時(shí)間,對應(yīng)約300個(gè)湍流特性時(shí)間(D/U0)。

1.7 驗(yàn)證

本文從兩方面對LES計(jì)算結(jié)果的準(zhǔn)確性進(jìn)行驗(yàn)證。圖4所示為子通道入口湍流發(fā)生器充分發(fā)展的近壁速度與Moser等[15]的DNS數(shù)據(jù)對比。近壁速度的采樣線從棒壁面到子通道中心,橫坐標(biāo)y+=uτy/ν為網(wǎng)格中心到壁面的無量綱距離,縱坐標(biāo)U+=U/uτ為無量綱速度。由圖4可見,LES計(jì)算的速度與DNS結(jié)果吻合較好,說明本文使用的OpenFOAM能夠準(zhǔn)確預(yù)測絕熱條件下的流體流動(dòng)特征,同時(shí)網(wǎng)格分辨率也滿足LES的計(jì)算要求。

圖4 近壁速度與DNS結(jié)果對比

另一方面,選取李志輝等[16]的圓管SCO2實(shí)驗(yàn)數(shù)據(jù)驗(yàn)證LES計(jì)算結(jié)果。圓管垂直布置,直徑為2 mm,長度為290 mm。CO2入口溫度為298.15 K,壓力為8.8 MPa,熱流密度為20 263 W/m2。圖5所示為Smagorinsky模型和WALE模型計(jì)算的壁面溫度和流體溫度與實(shí)驗(yàn)數(shù)據(jù)的對比。計(jì)算發(fā)現(xiàn),所選取的兩個(gè)模型均能準(zhǔn)確預(yù)測流體溫度沿軸向的變化,但是對壁面溫度的預(yù)測存在一定偏差。Smagorinsky模型計(jì)算得到的壁面溫度在進(jìn)出口段與實(shí)驗(yàn)數(shù)據(jù)吻合較好,但是無法預(yù)測傳熱惡化時(shí)的壁面溫度峰值。WALE模型計(jì)算得到的壁面溫度雖然整體偏高于實(shí)驗(yàn)數(shù)據(jù),但壁面溫度分布趨勢與實(shí)驗(yàn)完全一致,壁面溫度峰值與谷值位置也相同。因此,本文認(rèn)為OpenFOAM中的WALE模型能夠較為合理地預(yù)測超臨界流體的傳熱規(guī)律,選取WALE模型進(jìn)行后續(xù)計(jì)算研究,其所預(yù)測的壁面溫度偏差可能來源于恒定湍流普朗特?cái)?shù)的假設(shè)。

圖5 LES結(jié)果與實(shí)驗(yàn)數(shù)據(jù)對比

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

2.1 壁面溫度分布

圖6給出了圓管、環(huán)形通道和子通道的壁面溫度沿軸向的分布?？梢钥闯?3種通道壁面溫度分布趨勢大致相似,在入口附近存在溫度峰值,隨著流動(dòng)向下游發(fā)展,峰值逐漸消失,隨后壁面溫度沿流動(dòng)方向逐漸升高。環(huán)形通道壁面溫度峰值略高于圓管,但遠(yuǎn)高于子通道。SCO2在圓管和環(huán)形通道內(nèi)上升流動(dòng)時(shí),由于受到壁面的加熱作用,靠近壁面的流體溫度高于截面中心處的流體溫度,且近壁流體溫度高于擬臨界溫度,中心流體溫度低于擬臨界溫度。此時(shí),由于SCO2的物性劇烈變化,近壁流體密度減小,受到向上的浮升力作用,流速增大,速度梯度減小,從而降低了流體的切應(yīng)力;邊界層厚度逐漸增加,進(jìn)而抑制了湍流的擴(kuò)散,最終導(dǎo)致了傳熱惡化的發(fā)生,即壁面溫度出現(xiàn)峰值。由于子通道存在橫向交混,橫截面冷熱流體混合均勻,對壁面發(fā)熱量的輸運(yùn)能力更強(qiáng),因此壁面溫度峰值低于圓管和環(huán)形通道。在熱充分發(fā)展區(qū),圓管的壁面溫度最高,其次是環(huán)形通道,子通道壁面溫度最低。從壁面溫度沿軸向分布可以看出,子通道對流換熱效果優(yōu)于圓管和環(huán)形通道。

圖6 不同通道壁面溫度沿軸向分布

圖7 x=44D截面不同通道壁面溫度沿周向分布

2.2 平均流動(dòng)特性

a——圓管;b——環(huán)形通道;c——子通道

圖9 x=40D橫截面上沿徑向平均速度分布

2.3 湍流統(tǒng)計(jì)量

圖10 x=40D橫截面上沿徑向湍動(dòng)能分布

圖11 x=40D橫截面上軸向湍流熱流密度沿徑向分布

2.4 傳熱關(guān)系式評估

圖12所示為計(jì)算得到的圓管、環(huán)形通道和子通道努塞爾數(shù)(Nu)沿軸向的分布,橫坐標(biāo)為無量綱化的加熱段長度x/D。從圖12中可以看到,在通道前半段,環(huán)形通道的Nu與圓管相近,而在后半段環(huán)形通道的Nu分布在圓管和子通道之間。子通道Nu始終高于圓管和環(huán)形通道,說明子通道換熱能力強(qiáng)于相同D的圓管和環(huán)形通道。

圖12 Nu沿軸向分布

本文選取6個(gè)SCO2傳熱關(guān)系式[17-22](表1),將圓管、環(huán)形通道和子通道LES計(jì)算得到的Nu與傳熱關(guān)系式進(jìn)行比較,如圖13所示。Dittus-Boelter(D-B)公式預(yù)測的Nu與LES計(jì)算結(jié)果以及其他傳熱關(guān)系式的偏差很大,不能準(zhǔn)確預(yù)測SCO2的傳熱規(guī)律。這是由于D-B公式僅由Re和普朗特?cái)?shù)組成,因?yàn)镾CO2的熱物性在擬臨界溫度附近變化劇烈,需要對傳熱關(guān)系式進(jìn)行額外修正,以反映浮升力、流動(dòng)加速以及傳熱性質(zhì)變化的影響。表2總結(jié)了6種傳熱關(guān)系式的預(yù)測準(zhǔn)確度,其中參數(shù)δ和σ分別為平均相對誤差和標(biāo)準(zhǔn)差,定義為:

表1 選擇的典型傳熱關(guān)系式

表2 傳熱關(guān)系式預(yù)測偏差統(tǒng)計(jì)

a——圓管;b——環(huán)形通道;c——子通道

(16)

(17)

式中:n為LES計(jì)算的數(shù)據(jù)點(diǎn)個(gè)數(shù);δi為第i個(gè)相對誤差,其計(jì)算公式為:

(18)

式中:Nucorr,i為傳熱關(guān)系式計(jì)算得到的Nu;NuLES,i為LES計(jì)算得到的Nu。

從圖13和表2中可以看到,在選擇的6種傳熱關(guān)系式中,Kim等[20]提出的傳熱關(guān)系式對3種不同通道的傳熱預(yù)測都具有良好的準(zhǔn)確度,平均相對誤差僅為5%左右,最大標(biāo)準(zhǔn)差約為11%。Kim等[20]使用q+表征流動(dòng)加速效應(yīng),其關(guān)系式相比其他關(guān)系式中對流動(dòng)加速效應(yīng)做出了修正,且本研究的q+范圍與Kim關(guān)系式范圍接近,這可能是Kim關(guān)系式對LES計(jì)算結(jié)果預(yù)測準(zhǔn)確度高的原因。Kim關(guān)系式是基于豎直圓管實(shí)驗(yàn)數(shù)據(jù)擬合得到的,其對本文3個(gè)幾何結(jié)構(gòu)的LES計(jì)算結(jié)果都具有較好的預(yù)測效果,這可能與3個(gè)幾何結(jié)構(gòu)的D相同有關(guān)。因此,Kim關(guān)系式并不能體現(xiàn)流道幾何結(jié)構(gòu)對傳熱的影響。

3 結(jié)論

本文采用LES方法對SCO2在圓管、環(huán)形通道和子通道中的流動(dòng)與傳熱特性進(jìn)行了研究,對比了不同流道結(jié)構(gòu)的傳熱差異,并根據(jù)計(jì)算結(jié)果對典型傳熱關(guān)系式進(jìn)行了評估,主要結(jié)論如下。

1) 在相同熱力當(dāng)量直徑、相同工況條件下,子通道的傳熱性能最好,其次是環(huán)形通道和圓管;因此,選取圓管作為研究對象開展SCO2的傳熱特性研究具有保守性。

2) 子通道的傳熱性能優(yōu)于圓管和環(huán)形通道,但其近壁面湍動(dòng)能和湍流熱流密度最低,說明不同幾何通道湍流變量的強(qiáng)弱與傳熱性能優(yōu)劣并無直接對應(yīng)關(guān)系。

3) D-B公式不能準(zhǔn)確預(yù)測SCO2的傳熱強(qiáng)弱,而Kim關(guān)系式對LES計(jì)算得到的圓管、環(huán)形通道與子通道的努塞爾數(shù)都具有較高的預(yù)測精度。

需要指出的是,以上結(jié)論是通過對比分析一個(gè)混合對流傳熱工況所得,不一定適用于SCO2的其他傳熱模式,通道幾何結(jié)構(gòu)對SCO2傳熱特性的影響仍需更全面的深入研究。