C形換熱器一二次側(cè)流動傳熱耦合特性數(shù)值模擬研究

褚文俊，朱宏曄，成 翔，劉 洋，段日強(qiáng),*

(1.清華大學(xué) 核能與新能源技術(shù)研究院，先進(jìn)核能技術(shù)協(xié)同創(chuàng)新中心，先進(jìn)反應(yīng)堆工程與安全教育部重點實驗室，北京 100084；2.中國核動力研究設(shè)計院，四川 成都 610213)

C形換熱器是第三代反應(yīng)堆非能動余熱排出系統(tǒng)的主要設(shè)備。在反應(yīng)堆事故工況下，C形換熱器通過強(qiáng)迫或自然循環(huán)，將一次側(cè)介質(zhì)的衰變熱轉(zhuǎn)移到二次側(cè)(儲水箱)，防止堆芯超溫。研究其流動傳熱特性對反應(yīng)堆安全具有重要意義[1-4]。由于換熱器本身結(jié)構(gòu)復(fù)雜，傳熱管數(shù)量多，設(shè)備的換熱、流動阻力和傳熱管內(nèi)流量分配與設(shè)備結(jié)構(gòu)密切相關(guān)，且一、二次側(cè)之間的流動和傳熱狀態(tài)互相耦合，難以采用常規(guī)的一維分析方法進(jìn)行準(zhǔn)確評估[5-6]。薛若軍等[7]、Men等[8]通過CFD模擬了余熱排出系統(tǒng)的熱工水力現(xiàn)象；Zhang等[9]模擬了包含9個支管的C形棒束的瞬態(tài)傳熱，研究了RANS和LES等湍流模型對模擬結(jié)果的影響，并與實驗數(shù)據(jù)進(jìn)行了對比；劉愛瓊等[10]采用一維均相流模型計算了C形管內(nèi)冷凝換熱，并采用CFD程序分析水池空間的自然對流，建立了C形管換熱器的內(nèi)外耦合傳熱分析模型；陳薇等[11]對C形換熱器外流場分布進(jìn)行了研究。上述研究均對C形換熱器的形狀和排列進(jìn)行了重大簡化。此外，部分研究采用多孔介質(zhì)的方法對傳熱管束進(jìn)行模擬。Wright等[12]開發(fā)了CFD代碼PRHRCFD來分析AP1000的余熱排出系統(tǒng)，將含有689個C形管的換熱器以多孔介質(zhì)表示，并通過ROSA實驗數(shù)據(jù)驗證了該代碼的有效性；Ge等[13]對C形管進(jìn)行了數(shù)值模擬，采用多孔介質(zhì)模型和分布阻力法分析了傳熱管中流體速度和溫度的三維分布。在現(xiàn)有研究中，C形換熱器的全尺寸CFD模擬較少，對于一、二次側(cè)流動和傳熱耦合的研究不足，對傳熱管間流量分配的規(guī)律和機(jī)理缺乏認(rèn)識[14-15]。

本文擬建立C形換熱器流動和傳熱問題的全尺寸CFD模型，進(jìn)行不同熱工條件和幾何條件下的數(shù)值模擬研究，獲取流動參數(shù)分布和傳熱管間流量分配規(guī)律，分析流動參數(shù)和結(jié)構(gòu)參數(shù)對流量分配的影響機(jī)理，為C形換熱器的設(shè)計提供參考。

1 幾何模型與網(wǎng)格劃分

本文對C形換熱器進(jìn)口管道、進(jìn)口聯(lián)箱、傳熱管、出口聯(lián)箱和出口管道建立了全尺寸CFD計算模型，如圖1所示。

圖1 C形換熱器一次側(cè)流場幾何模型示意圖Fig.1 Schematic diagram of geometric model of primary side flow field of C-shaped heat exchanger

網(wǎng)格劃分時，將流場分為3部分：1) 進(jìn)口管道+進(jìn)口聯(lián)箱；2) 傳熱管(500根)；3) 出口管道+出口聯(lián)箱。其中，第1、3部分采用四面體非結(jié)構(gòu)化網(wǎng)格，第2部分采用六面體結(jié)構(gòu)化網(wǎng)格。各部分之間通過interface邊界進(jìn)行連接，并在壁面設(shè)置邊界層網(wǎng)格。網(wǎng)格結(jié)構(gòu)如圖2所示。

a——進(jìn)口管道與聯(lián)箱；b——出口管道與聯(lián)箱；c——管束截面

2 物理模型及仿真設(shè)置

2.1 一次側(cè)計算模型

1) 控制方程

本文采用三維穩(wěn)態(tài)k-ε兩方程湍流模型，控制方程如下。

連續(xù)方程：

(1)

動量方程：

(2)

其中，μt=ρCμk2/ε為湍流黏性系數(shù)，k為湍動能，ε為湍動能耗散率，由半經(jīng)驗公式確定。

能量方程：

(3)

2) 邊界條件

進(jìn)口為速度邊界，出口為定壓邊界，進(jìn)口管道、進(jìn)口聯(lián)箱、出口管道和出口聯(lián)箱壁面設(shè)置為無滑移壁面邊界且絕熱。傳熱管道壁面為無滑移壁面+對流換熱邊界，對流換熱系數(shù)為1～9 kW/(m2·K)，來流溫度設(shè)為100 ℃。

3) 物性參數(shù)

流體密度、黏性系數(shù)、導(dǎo)熱系數(shù)、比熱容等物性參數(shù)隨流體溫度變化較大，在計算中必須予以考慮。在本文計算中，一次側(cè)水壓力為15 MPa，一次側(cè)進(jìn)口水溫度為200～300 ℃；二次側(cè)水壓力為常壓，溫度為100 ℃。在此范圍內(nèi)，物性參數(shù)可擬合為以溫度為自變量的多項式函數(shù)：

ρ=701.219 7+1.920 3T-0.003 28T2

(4)

cp=13.489 86-0.044 55T+

5.344 68×10-5T2

(5)

λ=0.824 52-0.002 25T-8.670 81×

10-6T2-9.668 51×10-9T3

(6)

μ=0.001 49-4.827 28×

10-6T+4.154 23×10-9T2

(7)

其中：ρ、cp、λ、μ分別為密度、比熱容、導(dǎo)熱系數(shù)以及動力黏性系數(shù)；T為流體溫度。

金屬物性取100～300 ℃范圍內(nèi)的均值。

2.2 二次側(cè)計算模型

由于氣泡云和底部管束流量增大的影響，在中低熱流密度下，水平管束間池沸騰對流換熱系數(shù)高于單管池沸騰對流換熱系數(shù)。在管束間氣泡云速度分布與單管生成氣泡群速度分布相似的條件下，可通過以下方式計算管束外表面池沸騰對流換熱系數(shù)[16]：

(8)

αu=αone,nb+fαone,conv

(9)

系數(shù)f由管束尺寸和預(yù)期進(jìn)口速度決定，取值范圍為0.5～1。單管αone,nb由參考狀態(tài)下的換熱系數(shù)、熱流密度、壓力和壁面特性計算，忽略壁面影響，取實驗值和經(jīng)驗關(guān)系式，可得：

αone,nb=5.6×0.403×

0.101×q0.766 5=0.228q0.766 5

(10)

式(9)中對流換熱的貢獻(xiàn)αone,conv在湍流工況下由下式計算:

(11)

其中：g為重力加速度；αv為體積熱膨脹系數(shù)；ΔT為表面與環(huán)境溫度差；D為水力直徑；ν為動力黏度。通過代入0.1 MPa下飽和水的物性，得到對流換熱系數(shù)：

(12)

(13)

其中：q為管束中所有管表面的平均熱流密度，kW/m2；φ為管外換熱面積與光滑圓管外表面面積之比，對于外表無翅片的光管，φ=1。將式(10)、(12)、(13)代入式(8)，最終得到平均換熱系數(shù)與熱流密度的關(guān)系：

(1+(2+q)-1)

(14)

當(dāng)系數(shù)f在0.5～1區(qū)間變化時，可得到平均換熱系數(shù)的取值范圍。然而，由于熱流密度q為未知數(shù)，故在計算中需進(jìn)行迭代以獲得逼近真實值的平均換熱系數(shù)。

2.3 網(wǎng)格無關(guān)性檢驗

為檢驗網(wǎng)格尺寸對數(shù)值模擬的影響，采用不同網(wǎng)格(全局網(wǎng)格尺寸因子F分別為0.5、1.0和1.5)在給定邊界條件下進(jìn)行流場參數(shù)計算，獲得傳熱管內(nèi)流量的分布情況，如圖3所示。可觀察到，在全局網(wǎng)格尺寸因子F為0.5～1.5范圍內(nèi)，傳熱管內(nèi)橫向和縱向流量分布的差別不明顯。本研究選擇全局網(wǎng)格尺寸因子為1.0的網(wǎng)格，滿足計算精度需求。

圖3 不同網(wǎng)格流量分配計算結(jié)果比較Fig.3 Comparison of calculation results of different grid flow distributions

3 結(jié)果及分析

本文計算工況列于表1。為考察進(jìn)口管道與聯(lián)箱連接角度對流量分配的影響，本文計算了3種情況：1) 0°，流體自上而下進(jìn)入聯(lián)箱；2) 90°，流體水平進(jìn)入聯(lián)箱；3) 180°，流體自下而上進(jìn)入聯(lián)箱。

表1 計算工況熱工參數(shù)Table 1 Thermal parameter of calculation condition

3.1 流量分配特性

C形換熱器典型速度分布如圖4所示。從圖4a可看到，在中心對稱面(垂直于z軸)上，外層傳熱管道內(nèi)流速較高，內(nèi)層傳熱管道內(nèi)流速較低，差異較明顯。由于C形換熱器關(guān)于中心平面左右對稱，因此速度分布也關(guān)于中心平面左右對稱。以傳熱管進(jìn)口截面為參照，分別選取豎直方向中心列和水平方向中心行(圖4b)的傳熱管道進(jìn)行比較，其流量分布如圖5所示。由圖5可見，在同一行內(nèi)，傳熱管間流量變化幅度很小且左右對稱；而在同一列內(nèi)，傳熱管間流量變化非常明顯，總體呈現(xiàn)內(nèi)低外高的趨勢。在最內(nèi)幾層管道流量出現(xiàn)負(fù)值，說明出現(xiàn)了反向流動。

a——對稱面；b——傳熱器進(jìn)口截面(紅框內(nèi)為分析參考的中心列和中心行)進(jìn)口溫度200 ℃，流量2 kg/s

管束中心列和中心行分別有26根傳熱管，共51根(1根重疊)。可利用各管流量標(biāo)準(zhǔn)差來表示流量分布的不均勻程度，標(biāo)準(zhǔn)差越大，說明流量分配越不均勻。計算公式如下：

(15)

(16)

(17)

進(jìn)口溫度200 ℃，流量8 kg/s，二次側(cè)換熱系數(shù)5 kW/(m2·K)

3.2 流量的影響

進(jìn)口流量對傳熱管內(nèi)流量分配的影響如圖6、7所示。管束流量分配有以下變化規(guī)律：1) 中心行管道內(nèi)的流量逐漸增大，且中間位置的管道流量增長較快，形成中心向上凸起的流量分布(圖6a)；2) 中心列管道內(nèi)的流量整體逐漸增大，外層管道流量增長，內(nèi)層流量降低(絕對值增大，流向與外層相反)，流量分布曲線斜率增大，另外，中心位置的管道流量增大最快，同中心行管道流量分配情況類似，也形成向上的凸起(圖6b)；3) 中心行管道流量不均勻性較小，且與流量相關(guān)性不明顯(圖7)；4) 中心列管道流量不均勻性較大，隨流量增加不均勻性降低(圖7)。

圖6 C形換熱器中心行和中心列傳熱管流量Fig.6 Flow rate of heat transfer tubes in center row and center column of C-shaped heat exchanger

圖7 C形換熱器傳熱管流量標(biāo)準(zhǔn)差隨進(jìn)口流量的變化Fig.7 Flow variation standard deviation of C-shaped heat exchanger heat transfer tube changing with inlet flow

3.3 進(jìn)口溫度的影響

進(jìn)口溫度對流量分配的影響如圖8所示。進(jìn)口溫度對中心行管道流量幾乎沒有影響。但隨著進(jìn)口溫度的增加，中心列流量分布曲線斜率變大，即內(nèi)層管道反向流動流量更大、外層管道正向流動流量更大，中心列管道的流量不均勻性增大。

圖8 不同進(jìn)口溫度C形換熱器中心行和中心列傳熱管流量分布Fig.8 Flow rate distribution of center row and column of C-shaped heat exchanger heat transfer tubes with different inlet temperatures

3.4 進(jìn)口管道位置的影響

不同進(jìn)口位置的典型速度分布、進(jìn)口管道位置對流量分配的影響如圖9、10所示。可看到，盡管改變進(jìn)口管道位置能使聯(lián)箱內(nèi)流場發(fā)生明顯變化，但對傳熱管束的流量分配影響有限。在內(nèi)層管道中，仍存在反向流動現(xiàn)象。此外，通過比較不同進(jìn)口管道方位中間行、列流量標(biāo)準(zhǔn)差，可知在流量較低時，90°和180°設(shè)計的進(jìn)口管道流量差異高于0°的進(jìn)口管道；隨著流量增大，三者逐漸接近。這說明原設(shè)計(0°)的流量均勻性較另外兩種進(jìn)口管道設(shè)計方案的好。

3.5 C形換熱器流量分配機(jī)理分析

以上計算結(jié)果表明，C形換熱器傳熱管束流量分配存在以下特點：1) 同層管道之間流量分配基本一致；2) 外層管道流量較大，內(nèi)層管道流量為負(fù)值，即由出口聯(lián)箱流至進(jìn)口聯(lián)箱，各層管道流量由內(nèi)至外基本呈線性分布；3) 流量增大到一定值時，中間位置管道流量高于周圍管道流量。

造成不同層管道之間流量分配不均勻的主要原因是熱虹吸現(xiàn)象，如圖11所示。當(dāng)流體通過進(jìn)口聯(lián)箱進(jìn)入傳熱管束時，流經(jīng)最外側(cè)和最內(nèi)側(cè)的兩條路徑分別為O-A1-A2-E和O-B1-B2-B3-B4-E。由于管外換熱系數(shù)較高，當(dāng)流體通過傳熱管水平段O-A1和B1-B2時，一次側(cè)流體溫度迅速降至接近二次側(cè)流體溫度(100 ℃)，因此A1-B2’流體密度高于O-B1流體密度。而在B2’-B2-B1連線以下部位一次側(cè)流體溫度變化很小，可忽略不同路徑的密度差。這樣，路徑O-A1-A2-E和O-B1-B2-B3-B4-E之間存在密度差(圖11b)。根據(jù)進(jìn)口溫度的不同，密度差可達(dá)71～240 kg/m3，導(dǎo)致兩條路徑存在較大重力壓差，進(jìn)而出現(xiàn)熱虹吸現(xiàn)象，即流體有在O-A1-A2-E-B4-B3-B2-B1-O回路中流動的趨勢，在內(nèi)層管道出現(xiàn)反向流動。

a——溫度分布和可能流動路徑示意圖；b——密度分布進(jìn)口流量6 kg/s，進(jìn)口溫度250 ℃

進(jìn)口溫度250 ℃，流量6 kg/s

為進(jìn)一步驗證這一現(xiàn)象，將一、二次側(cè)設(shè)置為不換熱，各管道流量標(biāo)準(zhǔn)差大幅降低(Cx=15.87%，Cy=17.13%，Ct=16.66%)，如圖12所示。

圖10 不同進(jìn)口管道位置流量分布Fig.10 Flow distribution in different inlet pipe orientations

圖12 絕熱條件下中心行和中心列傳熱管流量分配Fig.12 Flow distribution of heat transfer tube in center row and center column under adiabatic condition

流量較大時，中間位置傳熱管流量較大，其原因是慣性力增強(qiáng)后主流方向與進(jìn)口管道方向一致，從而形成射流，指向管板中間位置，如圖13a所示。在到達(dá)管板后流體垂直于管板的速度分量轉(zhuǎn)化為靜壓，從而造成局部壓力升高，流量增大。而當(dāng)流速較低時，由于浮升力強(qiáng)于慣性力，主流受力彎曲向上，不能形成射流，因此沒有出現(xiàn)局部流量增大的現(xiàn)象(圖13b)。

a——流量10 kg/s，進(jìn)口溫度200 ℃；b——流量2 kg/s，進(jìn)口溫度300 ℃

3.6 C形換熱器流動阻力分布特性

中間列傳熱管流動阻力(壓差)計算結(jié)果如圖14所示。由于內(nèi)側(cè)傳熱管流量較小且存在反向流動，內(nèi)層傳熱管流動阻力較低，并由內(nèi)層到外層逐漸降低，然后再逐漸增大。當(dāng)進(jìn)口溫度為200 ℃時，阻力極小值管道出現(xiàn)在第5～6層；進(jìn)口溫度為300 ℃時，阻力極小值管道出現(xiàn)在第8～9層。由于中間部位管道流量存在1個峰值，反映在流動阻力分布上，也存在1個峰值，使得流量分布曲線呈W形。隨著溫度的升高，熱虹吸作用更加明顯，導(dǎo)致內(nèi)、外層管道流量差別增大，且流量絕對值增大，因此阻力分布曲線首尾兩端上翹更加明顯。

圖14 中間列傳熱管流動阻力分布Fig.14 Flow resistance distribution of heat transfer tubes in center row

4 結(jié)論

通過C形換熱器一、二次側(cè)耦合流動與傳熱的數(shù)值模擬，發(fā)現(xiàn)原設(shè)計傳熱管間流量分配非常不均勻，內(nèi)層流量低、外層流量高，且內(nèi)層存在反向流動現(xiàn)象。造成這一現(xiàn)象的原因是C形換熱器外層傳熱管在與進(jìn)口聯(lián)箱連接的水平段與二次側(cè)流體進(jìn)行換熱，溫度迅速降低、密度增大，導(dǎo)致在豎直方向上外層傳熱管內(nèi)重力壓差較大；而內(nèi)層重力壓差較低，形成熱虹吸現(xiàn)象。在進(jìn)口流量較低、進(jìn)口溫度較高時這一現(xiàn)象更為明顯。在本文研究的熱工參數(shù)范圍內(nèi)，熱虹吸現(xiàn)象一直占流量分配的主導(dǎo)地位。比較了0°(原設(shè)計)、90°和180° 3種進(jìn)口管道位置，發(fā)現(xiàn)進(jìn)口管道位置變化對流量不均勻性的改善不明顯；流量較低時，現(xiàn)有設(shè)計的流量不均勻性反而最優(yōu)。在進(jìn)一步的研究中，可探索降低流量分配不均勻性的措施，包括改變傳熱管與聯(lián)箱的連接方式，一、二次側(cè)傳熱系數(shù)分布，以及進(jìn)口聯(lián)箱的安裝方式等。