流道結構對超臨界二氧化碳傳熱特性的影響研究

白一凡,呂海財,王 漢

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

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

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

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

1 數值模擬設置

1.1 計算參數選取依據

首先,不同幾何流道的SCO2入口溫度T0相同,保證入口的熱物性參數相同;其次,入口位置的流動應滿足相似性準則的要求,即3個幾何流道的雷諾數相等:

(1)

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

(2)

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

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

(3)

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

(4)

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

根據以上分析,不同幾何流道的U0和qw應選擇相同的參數,同時選取適當幾何尺寸以保證D相同,這樣能夠使數值模擬結果的差異只來源于通道幾何結構。

1.2 控制方程

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

(5)

(6)

(7)

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

(8)

式中,PrSGS為SGS普朗特數,計算中采用OpenFOAM程序的默認常數1.0。

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

(9)

(10)

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

(11)

(12)

(13)

(14)

(15)

1.3 幾何結構和邊界條件

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

a——圓管;b——環形通道;c——子通道

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

1.4 計算工況

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

1.5 網格劃分

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

a——圓管;b——環形通道;c——子通道

1.6 計算方法

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

1.7 驗證

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

圖4 近壁速度與DNS結果對比

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

圖5 LES結果與實驗數據對比

2 計算結果分析

2.1 壁面溫度分布

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

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

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

2.2 平均流動特性

a——圓管;b——環形通道;c——子通道

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

2.3 湍流統計量

圖10 x=40D橫截面上沿徑向湍動能分布

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

2.4 傳熱關系式評估

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

圖12 Nu沿軸向分布

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

表1 選擇的典型傳熱關系式

表2 傳熱關系式預測偏差統計

a——圓管;b——環形通道;c——子通道

(16)

(17)

式中:n為LES計算的數據點個數;δi為第i個相對誤差,其計算公式為:

(18)

式中:Nucorr,i為傳熱關系式計算得到的Nu;NuLES,i為LES計算得到的Nu。

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

3 結論

本文采用LES方法對SCO2在圓管、環形通道和子通道中的流動與傳熱特性進行了研究,對比了不同流道結構的傳熱差異,并根據計算結果對典型傳熱關系式進行了評估,主要結論如下。

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

2) 子通道的傳熱性能優于圓管和環形通道,但其近壁面湍動能和湍流熱流密度最低,說明不同幾何通道湍流變量的強弱與傳熱性能優劣并無直接對應關系。

3) D-B公式不能準確預測SCO2的傳熱強弱,而Kim關系式對LES計算得到的圓管、環形通道與子通道的努塞爾數都具有較高的預測精度。

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