應(yīng)用于管殼式換熱器熱工水力數(shù)值模擬的多孔介質(zhì)算法

陸道綱，王 雨,*，袁 博，隋丹婷

(1.華北電力大學(xué) 核科學(xué)與工程學(xué)院，北京 102206；2.非能動(dòng)核能安全技術(shù)北京市重點(diǎn)實(shí)驗(yàn)室，北京 102206；3.長(zhǎng)江勘測(cè)規(guī)劃設(shè)計(jì)研究院，湖北 武漢 430010)

管殼式換熱器在工程中應(yīng)用十分廣泛，其類型多樣，內(nèi)部結(jié)構(gòu)各不相同，因此與實(shí)驗(yàn)研究相比，數(shù)值計(jì)算方法省時(shí)省力且具有可重復(fù)性。1974年，Patankar和Spalding[1]首次提出利用多孔介質(zhì)的概念對(duì)殼側(cè)流場(chǎng)進(jìn)行數(shù)值模擬，引入體積孔隙率并考慮分布式阻力對(duì)流場(chǎng)的影響。在此基礎(chǔ)上，Sha[2]引入了表面滲透率解釋說明了液態(tài)金屬換熱器中支撐板、傳熱管等對(duì)流場(chǎng)的阻礙作用。Prithiviraj和Andrews[3]開發(fā)了三維計(jì)算程序HEATX，該程序物理模型采用了Patankar提出的分布式阻力概念，結(jié)合表面滲透率和體積孔隙率對(duì)換熱器管側(cè)進(jìn)行建模。Ferng等[4]提出了一個(gè)三維計(jì)算流體力學(xué)沸騰模型研究蒸汽發(fā)生器二次側(cè)不同工況下的熱工水力特性。蒸汽發(fā)生器二次側(cè)傳熱管被看作多孔介質(zhì)，二次側(cè)兩相混合流被視為擬單相流，并假設(shè)其為物性參數(shù)可變化的可壓縮流。可見，許多學(xué)者在采用數(shù)值計(jì)算方法對(duì)管殼式換熱器相關(guān)問題進(jìn)行研究時(shí)，均利用了多孔介質(zhì)模型。尤其對(duì)于核電廠蒸汽發(fā)生器，其傳熱管數(shù)目眾多且結(jié)構(gòu)復(fù)雜，可運(yùn)用多孔介質(zhì)方法進(jìn)行適宜的簡(jiǎn)化模擬。但不同的網(wǎng)格劃分方法，每個(gè)控制體積的表面滲透率和體積孔隙率也會(huì)不同，所以多孔介質(zhì)系數(shù)的計(jì)算方案十分重要。

KAERI等[5]在CUPID的基礎(chǔ)上，開發(fā)了蒸汽發(fā)生器熱工水力分析程序CUPID-SG，該程序在模型的入口、中間和出口處分別取多孔介質(zhì)系數(shù)為0.0、0.5、0.9，3處的流體溫度也被設(shè)定為固定常數(shù)。程序與ATHOS3計(jì)算結(jié)果和瑞士FRIGG試驗(yàn)的結(jié)果進(jìn)行了對(duì)比，結(jié)果符合良好。但多孔系數(shù)并不符合真實(shí)情況下的數(shù)值，該程序仍需進(jìn)一步細(xì)化。鄧斌等[6]設(shè)計(jì)了一種計(jì)算多孔介質(zhì)特性參數(shù)的算法，將傳熱管截取控制體積單元邊界線的情況總結(jié)為16種位置關(guān)系。計(jì)算時(shí)對(duì)每個(gè)控制體積進(jìn)行掃描判斷其與傳熱管的位置關(guān)系，再選取相應(yīng)的計(jì)算公式進(jìn)行系數(shù)計(jì)算，但該方法相對(duì)較為繁瑣。

為能快速且真實(shí)地反映管殼式換熱器中兩相流的熱工水力特性，本文介紹一種適用于直角坐標(biāo)和柱坐標(biāo)系的大型管殼式換熱器各向異性的多孔介質(zhì)系數(shù)自動(dòng)生成程序。

1 基于多孔介質(zhì)模型的氣液兩相流控制方程

管殼式換熱器的殼側(cè)多為氣液兩相流動(dòng)，流動(dòng)情況復(fù)雜，為更真實(shí)反映流場(chǎng)流動(dòng)情況，本文在多孔介質(zhì)模型的基礎(chǔ)上結(jié)合兩相流模型給出了統(tǒng)一的控制方程(式(1))，包括氣相和液相的連續(xù)方程、動(dòng)量方程和能量方程。為使方程封閉，方程源項(xiàng)需輔助方程進(jìn)行計(jì)算。控制方程中各源項(xiàng)表達(dá)式列于表1[7]。

(1)

式中：α為空泡份額；ρ為密度，kg/m3；u、v、w分別為周向、徑向和軸向的流速，m/s；下標(biāo)k表示氣相或液相；fv、fθ、fr、fz為多孔系數(shù)，分別表示體積孔隙率及周向、徑向和軸向的表面滲透率；θ、r、z分別為周向、徑向和軸向的長(zhǎng)度，m；Sφ為方程源項(xiàng)；Γ為擴(kuò)散系數(shù)。

表1中：μ為流體動(dòng)力黏度，Pa·s；λ為導(dǎo)熱系數(shù)，W/(K·m)；cp為比定壓熱容，J/(kg·K)；T為溫度，K；p為壓力，Pa；M為相間質(zhì)量轉(zhuǎn)移速率，kg/(m3·s)；MX為由相間質(zhì)量轉(zhuǎn)移引起的動(dòng)量變化；R為分布式阻力；Q為二次側(cè)流體熱源。

表1 控制方程源項(xiàng)表達(dá)式Table 1 Governing equation source term expression

2 多孔系數(shù)的計(jì)算方案

多孔系數(shù)的計(jì)算參數(shù)主要是體積孔隙率和表面滲透率。以大亞灣核電站蒸汽發(fā)生器為原型，解釋說明多孔系數(shù)的計(jì)算方案。大亞灣蒸汽發(fā)生器內(nèi)部共有4 474根倒U型傳熱管[8]，以正方形排列成管束。傳熱管外徑為19.05 mm，管間距為27.43 mm，下筒體直徑為3 446 mm，管束套筒半徑為1 543 mm。

首先根據(jù)模型建立坐標(biāo)系網(wǎng)格，確定網(wǎng)格中心坐標(biāo)，并記錄網(wǎng)格尺寸、網(wǎng)格體積和柱坐標(biāo)系中3個(gè)方向的投影面積。網(wǎng)格與傳熱管的位置關(guān)系(橫截面)如圖1所示。圖中陰影區(qū)域即為氣液兩相流流通區(qū)域。

其次循環(huán)計(jì)算體積孔隙率和表面滲透率。對(duì)全局傳熱管中心坐標(biāo)進(jìn)行掃描，判斷是否落入有效區(qū)域。有效區(qū)域表示若傳熱管中心坐標(biāo)在此區(qū)域內(nèi)，則此控制體網(wǎng)格被該傳熱管截到。圖2中虛線內(nèi)部為有效區(qū)域，有效區(qū)域?yàn)榭刂企w網(wǎng)格四周各擴(kuò)展1個(gè)傳熱管半徑距離的區(qū)域，Δr為網(wǎng)格徑向尺寸，Δθ為網(wǎng)格周向角度，d為傳熱管外徑。

圖1 網(wǎng)格內(nèi)多孔介質(zhì)示意圖Fig.1 Schematic diagram of porous media in grid

圖2 網(wǎng)格有效區(qū)域示意圖Fig.2 Schematic diagram of grid effective area

2.1 計(jì)算體積孔隙率

從第1個(gè)傳熱管開始掃描，若循環(huán)中檢測(cè)到第n個(gè)傳熱管中心坐標(biāo)位于有效區(qū)域內(nèi)，則調(diào)用散點(diǎn)生成程序，該程序以此傳熱管中心坐標(biāo)為中心，在傳熱管直徑范圍內(nèi)均勻分布N個(gè)點(diǎn)，并記錄點(diǎn)坐標(biāo)。按順序掃描散點(diǎn)并判斷散點(diǎn)坐標(biāo)是否位于控制體內(nèi)。假設(shè)有M個(gè)坐標(biāo)點(diǎn)位于控制體內(nèi)，則此傳熱管占據(jù)控制體的體積V1(式中Δz為網(wǎng)格軸向尺寸)為：

(2)

逐次掃描所有傳熱管，判斷控制體網(wǎng)格是否被傳熱管截到。將在此網(wǎng)格內(nèi)的所有傳熱管體積加和，即為網(wǎng)格內(nèi)的所有傳熱管體積Vs。根據(jù)體積孔隙率定義計(jì)算fv。式(3)中下標(biāo)f、 s分別代表兩相流流體和傳熱管。

(3)

計(jì)算時(shí)，若要提高體積孔隙率的精確程度，可適當(dāng)增加傳熱管直徑范圍內(nèi)的散點(diǎn)數(shù)量。

2.2 計(jì)算表面滲透率

判斷位于網(wǎng)格有效區(qū)域內(nèi)的傳熱管數(shù)量不為0，那么則記錄位于控制體網(wǎng)格有效區(qū)域內(nèi)的傳熱管數(shù)目及徑向坐標(biāo)。

a——周向投影； b——徑向投影圖3 傳熱管投影示意圖Fig.3 Schematic diagram of heat transfer tube projection

根據(jù)表面滲透率定義，則周向表面滲透率為：

(4)

(5)

其中，θ、r分別為周向和徑向的方向矢量。

2.3 計(jì)算方案流程

圖4 計(jì)算方案流程圖Fig.4 Flow chart of calculation scheme

多孔系數(shù)的計(jì)算方案完整流程控制圖如圖4所示。該程序可用于直角坐標(biāo)和柱坐標(biāo)的大型管殼式換熱器各向異性的多孔介質(zhì)系數(shù)自動(dòng)計(jì)算，并且多孔系數(shù)會(huì)隨著網(wǎng)格生成方案的改變而更新。

3 計(jì)算方案驗(yàn)證

3.1 算例實(shí)施

選擇網(wǎng)格(3,7,1)作為算例，如圖5所示。選取1/4的網(wǎng)格劃分橫截面圖，(3,7,1)網(wǎng)格內(nèi)陰影區(qū)域?yàn)闅庖簝上嗔髁魍▍^(qū)域，其余為傳熱管。根據(jù)以上計(jì)算方案可得到該網(wǎng)格中心坐標(biāo)(θ,r,z)為(0.436 3 rad,0.733 7 m,0.025 0 m)，控制體尺寸為Δθ=0.174 5 rad，Δr=0.081 6 m，Δz=0.050 0 m，網(wǎng)格體積5.223 7×10-4m3。周向、徑向面積分別為Aθ=4.080 0×10-3m2，Ar=6.401 5×10-3m2。網(wǎng)格有效區(qū)域?yàn)?.428 7 rad<θ<0.443 9 rad、0.724 2 m

圖5 算例示意圖Fig.5 Schematic diagram of example

根據(jù)傳熱管中心坐標(biāo)位置判斷出有15根傳熱管位于有效區(qū)域內(nèi)，即(3,7,1)網(wǎng)格被15根傳熱管截到。以這15根傳熱管坐標(biāo)中心為中心均勻建立40個(gè)離散點(diǎn)，經(jīng)判斷有556個(gè)點(diǎn)落于控制體內(nèi)，那么實(shí)際固體所占體積為15根傳熱管的556/600，可求該控制體內(nèi)傳熱管體積為1.980 9×10-4m3，則該網(wǎng)格的體積孔隙率為0.620 3。15根傳熱管在周向和徑向上的投影完全覆蓋網(wǎng)格在該方向的面積，則其周向、徑向表面滲透率為0。

3.2 方案驗(yàn)證

1) 多孔系數(shù)驗(yàn)證

為驗(yàn)證多孔系數(shù)生成方案的正確性，本文以大亞灣蒸汽發(fā)生器為原型，應(yīng)用以上多孔系數(shù)計(jì)算程序計(jì)算在本文網(wǎng)格劃分方案下的體積孔隙率。共選取14個(gè)網(wǎng)格，是從(2,7,1)至(8,7,1)周向上連續(xù)7個(gè)網(wǎng)格和(2,8,1)至(8,8,1)周向上連續(xù)7個(gè)網(wǎng)格，計(jì)算孔隙率分別設(shè)為fv2、fv4；同時(shí)在CAD軟件中利用測(cè)量面積工具分別測(cè)量落在網(wǎng)格內(nèi)所有傳熱管的橫截面積以及該網(wǎng)格的軸向面積，進(jìn)而求得體積孔隙率，相應(yīng)網(wǎng)格的體積孔隙率設(shè)為fv3、fv5；且根據(jù)Yoon等[9-10]提到的算法計(jì)算平均孔隙率fv1。

根據(jù)文獻(xiàn)[9-10]，大亞灣蒸汽發(fā)生器的平均體積孔隙率為：

(6)

將3種結(jié)果進(jìn)行對(duì)比，對(duì)比結(jié)果如圖6所示。從圖6可看出，不同控制體的體積孔隙率互不相同，因?yàn)榭刂企w內(nèi)截到的傳熱管數(shù)目和位置關(guān)系不盡相同，導(dǎo)致體積孔隙率有所差異。但總體來說，多孔系數(shù)計(jì)算程序得到的體積孔隙率的值與CAD測(cè)量的真實(shí)值符合良好，最大誤差僅為4.5%。而且計(jì)算值和測(cè)量值基本均分布在平均孔隙率周圍。如若要提高計(jì)算精準(zhǔn)度，可通過增加傳熱管中心坐標(biāo)直徑范圍內(nèi)的散點(diǎn)數(shù)量來減小誤差。

圖6 體積孔隙率計(jì)算結(jié)果Fig.6 Calculation result of volume porosity

2) 數(shù)值模擬結(jié)果及驗(yàn)證

為進(jìn)一步研究GTG多孔系數(shù)計(jì)算方法在管殼式換熱器熱工水力數(shù)值模擬中的有效性，本文運(yùn)用此方法，對(duì)蒸汽發(fā)生器二次側(cè)管板至汽水分離器之間的流場(chǎng)進(jìn)行了三維瞬態(tài)數(shù)值模擬。3個(gè)方向入口流速分別為0.01、0.01和1.33 m/s，入口溫度為543.15 K，入口處空泡份額為0。設(shè)置壓力出口邊界，無滑移壁面邊界條件且壁面溫度保持不變，選取t=5 s時(shí)刻的流場(chǎng)進(jìn)行分析。

圖7為二次側(cè)流場(chǎng)液相和氣相流體的平均速度沿軸向變化示意圖。丁訓(xùn)慎[11]在其研究中提到管束出口處的汽水混合物流速為6.71 m/s， Li等[12]研究中給出汽體出口速度為6.5 m/s左右，液體出口速度為6 m/s左右。本文數(shù)值模擬結(jié)果顯示汽體出口速度為6.9 m/s，液體出口速度為6.2 m/s。參考以上出口速度，該模擬結(jié)果合理。在流體沸騰初期，產(chǎn)生的汽體受浮升力影響，且自身所受流阻較小，因此汽體速度上升較快。汽體沿軸向上升流動(dòng)時(shí)拖曳液體向上運(yùn)動(dòng)，所以氣液兩相流體速率變化基本相同。圖8示出了二次側(cè)平均溫度軸向變化圖以及對(duì)稱面上的溫度分布云圖。與Cong[13]利用漂移流模型研究蒸汽發(fā)生器二次側(cè)流場(chǎng)得到的溫度分布規(guī)律相似，由圖8可看出，二次側(cè)給水溫度在管殼兩側(cè)之間劇烈的傳熱影響下迅速上升，并且很快達(dá)到飽和溫度556 K，在沸騰達(dá)到穩(wěn)定后，流場(chǎng)溫度保持飽和溫度不變。從溫度分布云圖可看出，熱側(cè)流場(chǎng)溫度上升先于冷側(cè)。圖9示出了二次側(cè)壓力軸向變化圖以及軸向高度為7.575 m處相對(duì)于出口邊界的壓力分布云圖。入口壓力為6.891 MPa，出口壓力邊界為6.81 MPa,流場(chǎng)壓力在逐漸降低，總壓降為0.08 MPa，符合大亞灣蒸汽發(fā)生器的設(shè)計(jì)參數(shù)。從壓力云圖可看出，熱側(cè)的相對(duì)壓力略大于冷側(cè)，這是由于熱側(cè)換熱劇烈強(qiáng)于冷側(cè)，因此流體流速更快，導(dǎo)致熱側(cè)相對(duì)壓力大。

圖7 二次側(cè)流場(chǎng)速度分布Fig.7 Velocity distribution in secondary side

圖8 二次側(cè)溫度分布Fig.8 Temperature distribution in secondary side

圖9 二次側(cè)壓力分布Fig.9 Pressure distribution in secondary side

4 結(jié)論

在管殼式換熱器熱工水力數(shù)值模擬中，為使多孔介質(zhì)模型更真實(shí)準(zhǔn)確地反映殼側(cè)流場(chǎng)，本文提出了一種適用于管殼式換熱器三維熱工水力數(shù)值模擬的多孔系數(shù)計(jì)算方法GTG。該方法與網(wǎng)格和換熱管幾何參數(shù)有關(guān)，其基于網(wǎng)格參數(shù)和換熱管位置參數(shù)計(jì)算網(wǎng)格內(nèi)換熱管數(shù)目，進(jìn)而計(jì)算體積孔隙率；又使用區(qū)域縮短法計(jì)算控制體的表面滲透率。

運(yùn)用此方法，以大亞灣蒸汽發(fā)生器為原型，選取14個(gè)網(wǎng)格計(jì)算其多孔系數(shù)，并與利用CAD軟件測(cè)量的真實(shí)值進(jìn)行對(duì)比，結(jié)果表明計(jì)算值與真實(shí)值的最大偏差為4.5%。(6,7,1)網(wǎng)格多孔系數(shù)計(jì)算值和真實(shí)值偏差僅為0.4%，幾近相等。并且計(jì)算值基本分布在平均體積孔隙率周圍，計(jì)算值結(jié)果和真實(shí)值符合良好。又基于該方法對(duì)蒸汽發(fā)生器二次側(cè)流場(chǎng)進(jìn)行了數(shù)值模擬，并與同類研究結(jié)果進(jìn)行了對(duì)比驗(yàn)證，驗(yàn)證結(jié)果良好。總體來說，本文提出的計(jì)算多孔系數(shù)的方法是有效的。