AP1000核電機(jī)組巨型冷卻塔型體優(yōu)化數(shù)值計(jì)算

周蘭欣， 馬少帥， 弓學(xué)敏， 孫會(huì)亮

（1.華北電力大學(xué) 電站設(shè)備狀態(tài)監(jiān)測與控制教育部重點(diǎn)實(shí)驗(yàn)室，保定071003；2.保定華電電力設(shè)計(jì)研究院，保定071003）

近年來核電建設(shè)由沿海向內(nèi)陸地區(qū)擴(kuò)展，大批內(nèi)陸核電廠的選址、選型工作正在進(jìn)行.受內(nèi)陸地區(qū)水源限制，核電廠多采用閉式循環(huán)的自然通風(fēng)冷卻塔作為冷源設(shè)備.與燃煤電廠相比，核電機(jī)組循環(huán)水量大，所配備的冷卻塔淋水面積多在16 000m2以上，塔高超過165m，此類規(guī)模的巨型冷卻塔在常規(guī)火電機(jī)組中至今仍無應(yīng)用.巨型冷卻塔的結(jié)構(gòu)參數(shù)直接影響循環(huán)水的冷卻溫降，在設(shè)計(jì)選型過程中，以往規(guī)范標(biāo)準(zhǔn)所推薦的經(jīng)驗(yàn)值有待進(jìn)一步論證.為了尋找最優(yōu)結(jié)構(gòu)參數(shù)，需要對塔內(nèi)流場、溫度場進(jìn)行數(shù)值計(jì)算.

作為第三代核電技術(shù)的AP1000機(jī)組具有較高的安全可控性，是未來我國核電建設(shè)的主流機(jī)型之一.根據(jù)內(nèi)陸中部省份地區(qū)的年平均氣象參數(shù)得知，年平均氣溫為17℃，當(dāng)冷卻塔循環(huán)倍率采用50時(shí)，AP1000核電機(jī)組的循環(huán)水量為180 000t／h，所需淋水面積為20 000m2.針對這一淋水面積的冷卻塔，采用CFD三維數(shù)值計(jì)算軟件對雙曲線塔體進(jìn)行了模擬計(jì)算，以確定最優(yōu)結(jié)構(gòu)參數(shù).

在設(shè)計(jì)施工過程中，雙曲線冷卻塔外型結(jié)構(gòu)根據(jù)其塔筒旋轉(zhuǎn)殼體的母線方程進(jìn)行尺寸放樣定位，如圖1所示，將橫坐標(biāo)軸取在喉部，其母線方程為

式中：r為塔筒半徑；r0為喉部半徑；h為坐標(biāo)點(diǎn)處到喉部的距離（向上為正，向下為負(fù)）；α為決定雙曲線曲率的特征值［1］.

圖1 塔體結(jié)構(gòu)尺寸Fig.1 Structural parameters of the cooling tower

冷卻塔的高度由熱力計(jì)算得出，原則是塔筒抽力應(yīng)克服塔內(nèi)阻力，本文中冷卻塔塔高為210m.針對AP1000機(jī)組，當(dāng)淋水面積20 000m2確定時(shí)，決定雙曲線塔筒外型的參數(shù)是母線方程中的曲率特征值α及喉部半徑r0.以往相關(guān)規(guī)范標(biāo)準(zhǔn)建議塔筒高度不高于165m，并給出了塔體相關(guān)結(jié)構(gòu)尺寸的經(jīng)驗(yàn)推薦值［2］.顯然，在巨型塔的設(shè)計(jì)、建造過程中，需要對以往經(jīng)驗(yàn)值進(jìn)一步驗(yàn)證.

以往常規(guī)火電機(jī)組冷卻塔在選型設(shè)計(jì)時(shí)母線方程中的曲率特征值α常取經(jīng)驗(yàn)值0.15～0.17.相關(guān)標(biāo)準(zhǔn)中喉部面積與塔底面積比的推薦值為0.30～0.40［2］.筆者采用Fluent三維數(shù)值模擬軟件對核電站用巨型冷卻塔母線方程中的曲率特征值α、喉部半徑r0的取值重新進(jìn)行計(jì)算驗(yàn)證.α取值為0.14、0.15、0.16、0.17和0.18共5種情況，喉部半徑與塔底半徑比r0／r2取值為0.5、0.6、0.7和0.8共4種情況.

1 計(jì)算方法

計(jì)算流體動(dòng)力學(xué)（CFD）是建立在經(jīng)典流體力學(xué)與數(shù)值計(jì)算方法基礎(chǔ)上的一門新興獨(dú)立學(xué)科，F(xiàn)luent作為大型商業(yè)化CFD軟件，目前被廣泛應(yīng)用于數(shù)值計(jì)算領(lǐng)域.該軟件提供了大量計(jì)算流體相關(guān)模型和算法，并設(shè)置了外接用戶端口（UDF），具有穩(wěn)定性好、適用范圍廣、精度高和拓展性好等優(yōu)點(diǎn)［3］.因此，選用Fluent作為冷卻塔數(shù)值模擬平臺(tái).

在濕式自然通風(fēng)冷卻塔三維數(shù)值計(jì)算中，通常分噴濺區(qū)、填料區(qū)和雨區(qū)3個(gè)區(qū)域?qū)諝馀c冷卻水之間的傳熱傳質(zhì)過程進(jìn)行模擬.由于冷卻塔內(nèi)水和空氣的體積比小于10%，在噴濺區(qū)和雨區(qū)采用離散相模型計(jì)算塔內(nèi)流場［4］，其中水滴作為離散相采用拉格朗日法計(jì)算，空氣作為連續(xù)相采用歐拉法求解.填料區(qū)采用自定義函數(shù)添加源項(xiàng)的方法進(jìn)行模擬［5－11］，過程是根據(jù)填料層的熱力特性及氣水熱質(zhì)傳遞原理計(jì)算出冷卻水的蒸發(fā)量和溫度降低值，相關(guān)計(jì)算過程可參閱文獻(xiàn)［1］，再通過自定義源項(xiàng)的方式加入到連續(xù)相的主控方程中，實(shí)現(xiàn)對填料層氣水熱質(zhì)傳遞的耦合計(jì)算［4］.

Fluent求解器采用分離隱式算法，流場計(jì)算采用典型的Simple算法.計(jì)算中，采用穩(wěn)態(tài)雷諾應(yīng)力平均N－S方程，選用標(biāo)準(zhǔn)k－ε湍流模型，在輸運(yùn)方程中考慮了浮力項(xiàng).

1.1 濕空氣控制方程

冷卻塔內(nèi)濕空氣作為連續(xù)相，當(dāng)機(jī)組在穩(wěn)定工況下運(yùn)行時(shí)，其內(nèi)外流場可以當(dāng)作穩(wěn)態(tài)計(jì)算，通用控制方程如下［4］

式中：ρ為空氣密度；ui為速度矢量；Φ為通用變量，分別表示各向速度分量（u、v、w）、濕空氣組分Y、溫度T、湍動(dòng)能k和湍流耗散率ε；ΓΦ為廣義擴(kuò)散系數(shù)；SΦ為廣義源項(xiàng).

1.2 水滴控制方程

冷卻塔內(nèi)循環(huán)水以水滴形式在噴濺區(qū)和雨區(qū)自由下落，可采用拉格朗日方法計(jì)算其流場.離散相模型可追蹤水滴的運(yùn)動(dòng)軌跡，并耦合計(jì)算水滴與濕空氣間的熱質(zhì)傳遞.水滴溫度變化方程［4］：

式中：Mp、Tp、Tadb、Ap分別為控制單元內(nèi)水滴質(zhì)量、水滴溫度、氣相干球溫度和水滴面積；cp為水滴的比定壓熱容；h、hfg和t分別為傳熱系數(shù)、水的汽化潛熱和經(jīng)歷時(shí)間.

水滴蒸發(fā)速率為：

式中：cs為水滴表面蒸汽的物質(zhì)的量濃度；c∞為濕空氣中水蒸氣的物質(zhì)的量濃度；Mw為水的摩爾質(zhì)量；hD為以濃度差為動(dòng)力的傳質(zhì)系數(shù).

h、hD由 Nusslet關(guān)系式得出［4］：

式中：Re為雷諾數(shù)；Prva為濕空氣的普朗特?cái)?shù)；Sc為施密特?cái)?shù).

當(dāng)Fluent耦合計(jì)算氣水間熱質(zhì)傳遞時(shí)，相間能量、質(zhì)量和動(dòng)量傳遞將作為源項(xiàng)添加到主控方程中.

1.3 阻力計(jì)算

熱態(tài)的冷卻塔阻力包括氣水兩相間作用力和結(jié)構(gòu)阻力，其中相間作用力對連續(xù)相壓降的影響通過離散相模型實(shí)現(xiàn)，結(jié)構(gòu)阻力主要包括進(jìn)風(fēng)口阻力、填料阻力、配水管網(wǎng)阻力和收水器阻力等4部分.由于這些細(xì)密的結(jié)構(gòu)體在Gambit中建模困難，因此它們對連續(xù)相的作用作為附加動(dòng)量匯添加到主控方程中，通過用戶自定義函數(shù)UDF編程實(shí)現(xiàn).

動(dòng)量匯方程一般形式如下［4］：

式中：Vp為空氣通過幾何邊界面的垂直速度分量；K為壓力損失系數(shù).

K由基于試驗(yàn)的經(jīng)驗(yàn)關(guān)系式計(jì)算得到［11］

式中：mw為水滴質(zhì)量流量；ma為濕空氣質(zhì)量流量；L為填料層高度.

Kjfk、Kgw、Kssq分別為進(jìn)風(fēng)口、管網(wǎng)和收水器的壓力損失系數(shù)，根據(jù)有關(guān)試驗(yàn)資料可以取定值［12］，Kjfk＝0.5、Kgw＝0.5和Kssq＝3.5.

1.4 計(jì)算邊界及網(wǎng)格劃分

以AP1000機(jī)組為例，冷卻塔高210m，進(jìn)風(fēng)口高度12.5m，淋水面積20 000m2，喉部半徑與塔底半徑比依次選用0.5、0.6、0.7、0.8等4種情況，雙曲線塔體母線方程中的特征值α依次選用0.14、0.15、0.16、0.17和0.18等5種情況.塔筒外型依據(jù)其母線方程確定，計(jì)算環(huán)境區(qū)域?yàn)楦叨?00m、直徑500m的圓柱體.

圖2給出了計(jì)算邊界條件.靜風(fēng)條件下進(jìn)出口均為壓力邊界條件；塔壁、地面為無滑移壁面邊界條件；離散相邊界在進(jìn)出口、水池面、地面均為逃逸邊界條件.冷卻水在噴淋層采用面射流模型引入冷卻塔內(nèi)，水滴當(dāng)量直徑為3mm，對水滴下落至水池面的運(yùn)動(dòng)軌跡進(jìn)行追蹤.

圖2 計(jì)算邊界條件Fig.2 Computational boundary conditions

網(wǎng)格劃分如圖3所示，利用尺寸函數(shù)工具Sizefunction，采用非結(jié)構(gòu)化四面體網(wǎng)格劃分，并對網(wǎng)格無關(guān)性進(jìn)行驗(yàn)證，計(jì)算了140萬、160萬、180萬網(wǎng)格總數(shù)下的平均出塔水溫.160萬、180萬網(wǎng)格數(shù)量下的塔內(nèi)流場無明顯變化，所得平均出塔水溫相近，考慮到計(jì)算時(shí)間，選用網(wǎng)格總數(shù)為162萬.

圖3 冷卻塔網(wǎng)格劃分Fig.3 Grid division for the cooling tower

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

2.1 冷卻塔內(nèi)流場分析

經(jīng)計(jì)算，得到三維冷卻塔x－y截面上的速度等值線分布圖.喉部半徑與塔底半徑比r0／r2＝0.6、曲率特征值α＝0.16時(shí)的速度等值線分布如圖4所示，靜風(fēng)環(huán)境下塔內(nèi)流場呈軸對稱分布.空氣在塔內(nèi)吸熱后溫度升高，密度小于塔外空氣密度，這種密度差使塔內(nèi)產(chǎn)生抽吸力.冷卻塔外氣流在抽吸力作用下經(jīng)進(jìn)風(fēng)口流入塔內(nèi)，受進(jìn)風(fēng)口阻力、填料阻力和雨區(qū)阻力的作用，風(fēng)速逐漸降低.

圖5給出了不同喉部半徑與塔底半徑比下的填料層風(fēng)速值.由圖5可以看出，隨著喉部半徑與塔底半徑比的增大，塔內(nèi)風(fēng)速也逐漸增大，當(dāng)r0／r2＝0.7時(shí)，風(fēng)速達(dá)到最大值1.53m／s，當(dāng)r0／r2超過0.7后，流場發(fā)生惡化，風(fēng)速降低.過大的喉部半徑與塔底半徑比使氣流通過面積增加，速度下降，動(dòng)壓降低，不利于發(fā)揮雙曲線塔筒對氣流的抽吸作用，故r0／r2＝0.7為一臨界值.

圖4 當(dāng)r0／r2＝0.6、α＝0.16時(shí)，x－y 截面上的速度等值線分布圖Fig.4 Contours of air velocity in the x－y plane at r0／r2＝0.6，α＝0.16

圖5 填料層風(fēng)速隨喉部半徑與塔底半徑比的變化Fig.5 Air velocity in fill region vs.r0／r2

2.2 喉部半徑與塔底半徑比對進(jìn)塔風(fēng)量和蒸發(fā)水量的影響

圖6 進(jìn)塔風(fēng)量隨喉部半徑與塔底半徑比的變化Fig.6 Inlet air flow vs.r0／r2

圖6為不同喉部半徑與塔底半徑比下進(jìn)塔風(fēng)量的統(tǒng)計(jì)值.由圖6可以看出，進(jìn)塔風(fēng)量隨喉部半徑與塔底半徑比的增大而增大，這是因?yàn)楫?dāng)淋水面積確定時(shí)，喉部尺寸的增大使冷卻塔風(fēng)筒的容積增大，空氣通過量相應(yīng)增加.當(dāng)r0／r2＝0.7時(shí)，冷卻塔內(nèi)氣水比為0.82，此時(shí)填料區(qū)的容積質(zhì)散系數(shù)達(dá)到最大值，由水向空氣的傳熱傳質(zhì)最充分，蒸發(fā)水量也達(dá)到最大值.圖7給出了冷卻塔內(nèi)蒸發(fā)水量隨喉部半徑與塔底半徑比的變化，質(zhì)量傳遞的趨勢可從圖中看出.

圖7 蒸發(fā)水量隨喉部半徑與塔底半徑比的變化Fig.7 Water evaporation rate vs.r0／r2

2.3 塔內(nèi)溫度場和平均出塔水溫分析

圖8為不同喉部半徑與塔底半徑比下水池面的溫度分布圖.由圖8可以看出，中心區(qū)域溫度最高，外圍溫度較低.隨著喉部半徑與塔底半徑比的增大，中心區(qū)域及外圍的溫度均有所下降，在r0／r2＝0.7時(shí)溫度達(dá)到最低值，此后繼續(xù)增大喉部半徑與塔底半徑比，水池面溫度反而有所回升.

圖8 水池面溫度分布（單位：K）Fig.8 Contours of water temperature on pool surface（unit：K）

表1為相同運(yùn)行工況下冷卻塔采用不同喉部半徑與塔底半徑比時(shí)計(jì)算所得的平均出塔水溫.由表1可以看出，平均出塔水溫隨喉部半徑與塔底半徑比的增大呈現(xiàn)先降低后升高的趨勢.空氣作為冷源對塔內(nèi)較高溫度的噴淋水進(jìn)行冷卻，由2.2節(jié)分析可知塔內(nèi)溫度的回升是由于進(jìn)塔風(fēng)量的減小造成的.水滴下落時(shí)間統(tǒng)計(jì)結(jié)果如圖9所示，水滴下落時(shí)間越長，其與空氣間的換熱就越充分，平均出塔水溫就越低.

表1 不同喉部半徑與塔底半徑比下的平均出塔水溫Tab.1 Average outlet water temperature vs.r0／r2

圖9 水滴下落時(shí)間隨喉部半徑與塔底半徑比的變化Fig.9 Drop time of droplet vs.r0／r2

2.4 曲率特征值α對平均出塔水溫的影響

表2為r0／r2＝0.6、相同運(yùn)行工況下，塔筒母線方程選用不同的曲率特征值α?xí)r計(jì)算所得的平均出塔水溫，可以看出隨α的增大，平均出塔水溫降低.

表2 r0／r2＝0.6時(shí)，不同α下的平均出塔水溫Tab.2 Outlet water temperature vs.αat r0／r2＝0.6℃

2.5 綜合分析α及r0／r2對平均出塔水溫和進(jìn)塔風(fēng)量的影響

圖10為不同α和r0／r2組合方案下計(jì)算所得的平均出塔水溫.由圖10可以看出，當(dāng)r0／r2一定時(shí)，隨特征值α的增大，4種喉部半徑取值下的平均出塔水溫均呈現(xiàn)下降趨勢，且均在α＝0.18處取得最低值.

圖10 綜合α及r0／r2對平均出塔水溫的影響Fig.10 Comprehensive influence ofαand r0／r2on average outlet water temperature

表3給出了不同α及r0／r2組合下計(jì)算所得的進(jìn)塔風(fēng)量.由表3可知，進(jìn)塔風(fēng)量的變化趨勢與出塔水溫的變化趨勢相一致，在α＝0.18處達(dá)到最大值.同時(shí)可以看出，在計(jì)算范圍內(nèi)的任意特征值α下，r0／r2＝0.7均為最佳喉部半徑尺寸，此時(shí)進(jìn)塔風(fēng)量最大，出塔水溫最低.當(dāng)r0／r2大于0.7時(shí)，出塔水溫開始回升.

表3 綜合α及r0／r2對進(jìn)塔風(fēng)量的影響Tab.3 Comprehensive influence ofαand r0／r2on inlet air flow kg／s

3 結(jié) 論

（1）進(jìn)塔風(fēng)量隨喉部半徑與塔底半徑比的增大而增大，當(dāng)r0／r2＝0.7時(shí)進(jìn)塔風(fēng)量達(dá)到最大值，此時(shí)平均出塔水溫最低，與r0／r2＝0.5時(shí)相比，平均出塔水溫降低了0.53K，當(dāng)r0／r2大于0.7時(shí)，平均出塔水溫又開始回升.

（2）當(dāng)r0／r2確定時(shí)，平均出塔水溫隨曲率特征值α的增大而降低，α＝0.18時(shí)的平均出塔水溫比α＝0.14時(shí)的平均出塔水溫低0.18K.

（3）在選型設(shè)計(jì)時(shí)，AP1000核電機(jī)組用冷卻塔的喉部半徑與塔底半徑比r0／r2應(yīng)取為0.6～0.7，塔筒母線方程中的曲率特征值α應(yīng)取0.17～0.18.

［1］趙振國.冷卻塔［M］.2版.北京：中國水利水電出版社，2001：36－37.

［2］建設(shè)部.GB／T 50102—2003工業(yè)循環(huán)水冷卻設(shè)計(jì)規(guī)范［S］.北京：中國計(jì)劃出版社，2003.

［3］韓占忠.流體工程仿真計(jì)算實(shí)例與應(yīng)用［M］.北京：北京理工大學(xué)出版社，2004.

［4］Fluent Inc.User＇s guide［M］.Lebanon，USA：Fluent Inc，2003.

［5］張志剛，王瑋，曾德良，等.冷卻塔出塔水溫的迭代計(jì)算方法［J］.動(dòng)力工程學(xué)報(bào)，2010，30（5）：372－377.ZHANG Zhigang，WANG Wei，ZENG Deliang，et al.The iterative calculation of cooling tower outlet water temperature［J］.Journal of Chinese Society of Power Engineering，2010，30（5）：372－377.

［6］趙元賓，孫奉仲，王凱，等.十字墻濕式冷卻塔冷卻特性的數(shù)值研究［J］.中國電機(jī)工程學(xué)報(bào)，2009，29（8）：6－13.ZHAO Yuanbin，SUN Fengzhong，WANG Kai，et al.Numerical analysis of cooling performance of wet cooling tower with cross wall［J］.Proceedings of the CSSEE，2009，29（8）：6－13.

［7］POPPE M，ROGENER H.Berechnung von Ruckkuhlwerken［J］.VDI－Warmeatlas，1991（1）：1－15.

［8］WILLIAMSON N，ARMFIELD S，BEHNIA M.Numerical simulation of flow in a natural draft wet cooling tower－the effect of radial thermo fluid fields［J］.Applied Thermal Engineering，2008，28（2）：178－189.

［9］AL－WAKED R，BEHINA M.CFD simulation of wet cooling tower［J］.Applied Thermal Engineering，2006，26（4）：382－395.

［10］KLIMANEK A.Solution of heat and mass transfer in counter flow wet－cooling tower fills［J］.International Communications in Heat and Mass Transfer，2009，36（6）：547－553.

［11］KLOPPERS J C，KROGER D G.Loss coefficient correlation for wet cooling tower fills［J］.Applied Thermal Engineering，2003，23（17）：2201－2211.

［12］周蘭欣，蔣波，陳素敏.自然通風(fēng)濕式冷卻塔熱力特性數(shù)值模擬［J］.水利學(xué)報(bào)，2009，40（2）：208－213.ZHOU Lanxin，JIANG Bo，CHEN Sumin.Numerical simulation for thermal performance of natural draft wet cooling tower［J］.Journal of Hydraulic Engineering，2009，40（2）：208－213.