運行狀態下超大型冷卻塔內表面風荷載的數值模擬研究

董國朝 ，張建仁，蔡春聲,2，韓 艷,李春光

(1.長沙理工大學 土木與建筑學院，湖南 長沙 410114； 2. 路易斯安娜州立大學 土木與環境工程系，美國路易斯安娜州 巴吞魯日 70803)

運行狀態下超大型冷卻塔內表面風荷載的數值模擬研究

董國朝1?，張建仁1，蔡春聲1,2，韓 艷1,李春光1

(1.長沙理工大學 土木與建筑學院，湖南 長沙 410114； 2. 路易斯安娜州立大學 土木與環境工程系，美國路易斯安娜州 巴吞魯日 70803)

對運行狀態下的某超大型雙曲冷卻塔的內表面平均風壓進行了CFD數值模擬.在計算流體動力學軟件基礎上進行二次開發，采用DPM模型結合UDF函數方法加入源項來研究某超大型冷卻塔內表面平均風壓分布；塔中水相采用了拉格朗日方法模擬，而空氣相采用歐拉方法模擬，較好地實現了冷卻塔運行狀態下的內外流場計算及其與傳熱傳質的耦合計算，分析了運行狀態下冷卻塔橫風向來流時的內壓分布規律.無側風工況下計算結果顯示，塔運行過程中的內表面壓力對稱性良好，出水溫度與實測結果相符，驗證了本文提出的塔運行過程中的傳熱傳質計算方法的正確性.側風工況下得到塔內壓力系數沿高度方向相應變大，而沿緯向變化不明顯.同時討論了中國規范對內表面壓力系數的取值不完善之處，給出了建議取值，為超大型冷卻塔設計過程中的內壓計算提供方法和依據.

冷卻塔；平均風荷載；數值模擬；傳熱傳質；離散相模型(DPM)

超大型雙曲自然通風冷卻塔是一種大型空間薄壁開口結構，普遍用于發電廠中循環水冷卻，風荷載對其安全性的影響起決定性作用.隨著電站裝機容量的增加，冷卻塔也向高大化發展.如湖南、湖北擬建和在建的塔高分別達到200 m和220 m.

隨著英國渡橋電廠冷卻塔的倒塌，工程界對于冷卻塔的風荷載給予了極大關注[1-5]，并制定了相關設計規范.而英國規范和德國規范則相對完善，并被他國參考.由于塔體包括內、外兩個表面，因此涉及到外風壓載荷以及內風壓載荷，以往的許多研究主要是針對塔的外表面風壓，而對于內表面的風壓研究較少.孫天風等[2]對中國的茂名塔進行了內壓實測，結果顯示，迎風面受到的內外壓的作用力方向一致，兩者的疊加作用將導致子午向的應力顯著增大.且有資料[3]表明，英國渡橋電廠冷卻塔倒塌，主要是由于迎風面子午向鋼筋受拉屈服斷裂造成的.因此，在冷卻塔風荷載設計過程中必須充分考慮內壓的影響.

李鵬飛等[1]認為塔內表面壓力系數沿環向、高度均勻分布，假定其為某個數值.Kawarabata等[6]認為內壓壓力系數為-0.4～-0.50，實際設計中采用內壓壓力系數為-0.45.德國規范[7]將內表面的風壓系數取定值為-0.50.中國規范只對165 m高度以下的雙曲型冷卻塔外表面的風壓進行了規定[8-9]，而對于內壓，規范中并未提及其具體設計標準.內壓研究正被越來越多的學者所重視，但是冷卻塔在運行狀態下的內壓系數研究尚少見文獻涉及.不同的學者和文獻的各種取值差異較大，而如何借鑒各國的規范來指導中國的設計，有必要做充分的研究.然而，內壓研究存在其固有困難：現場實測由于受到外部天氣因素及塔內運行的影響，難于得到最不利工況結果；實驗室測試卻無法重現其真實運行過程中的填料壓降以及上升抽力等條件.因此，數值模擬由于其條件易于實現的優越性成為其中重要的研究手段.

本文采用CFD數值模擬方法，對商業軟件FLUENT進行二次開發，采用其中的離散相(Discrete Phase Model，DPM)模型，結合用戶自定義函數(User Defined Functions，UDF)，在冷卻塔填料區域中加入熱水程序，采用歐拉方法對空氣相進行模擬，采用拉格朗日方法對水相進行模擬，在填料區域中水流流動的水膜性質通過給定速度的液滴流來近似，通過熱水與空氣之間的熱交換來模擬冷卻塔在實際運行過程中塔內發生的傳熱傳質過程.得到一套反映真實運行狀態下冷卻塔的內壓確定方法，并且研究其在50年重現期風速下的內表面風壓系數分布規律及流場特性，為冷卻塔的內壓確定提供依據和方法，并指導其設計.

1 自然通風逆流濕式冷卻塔簡介

自然通風逆流濕式冷卻塔是一種空氣和熱循環水混合接觸式換熱設備，結構如圖1所示.從冷凝器中出來的熱水進入塔內淋水系統，淋水系統將熱水噴灑在下面的多孔介質填料中，熱水經過填料的時候與空氣發生傳熱傳質以后成雨狀落入水池中，其主要傳熱傳質區由配水區、填料區和雨區組成.熱循環水經配水噴嘴噴出后分別以水滴、水膜和水滴的形式依次經配水區、填料區和雨區進行對流傳熱傳質，然后落入底部蓄水池并由循環水泵輸送回凝汽器內循環再用.作為冷源的環境空氣經塔入口依次進入雨區、填料區和配水區，變為溫度較高、密度較小的近乎飽和的熱濕空氣而沿著塔往上運動，在塔內形成上升的抽力，最后攜帶著熱量被排放到大氣中，整過過程中不需要任何機械通風設備，所以稱為自然通風冷卻塔.

以往文獻[10]關于冷卻塔內的傳熱傳質研究主要是針對塔內工藝方面，而通過模擬塔內傳熱傳質來研究塔內壓力分布規律的相對較少.

2 計算理論及方法

在FLUENT中，描述熱、質量和動量傳輸的空氣(連續相)流動方程可以寫成以下統一形式：

▽ (ρmauφ-Γφ▽φ)=Sφ+Spφ.

(1)

式中：ρma為潮濕空氣密度；u為速度矢量；φ為標量如U,V,W,T,YV,k和ω；Γφ為耗散系數；Sφ為空氣相的源項，Spφ為空氣和水滴之間相互作用的附加源項.

圖1 自然通風逆流濕式冷卻塔

離散相水滴的計算包括水滴軌跡和速度的計算，水滴的熱和質量傳輸計算以及解決對流動狀態有影響的兩相之間的耦合計算.其中水滴軌跡計算根據拉格朗日參考體系，相對于水滴速度的軌跡運動方程為：

(2)

式中：rP為水滴的軌跡；up為水滴的瞬時速度.水滴軌跡的計算通過對水滴所受到的力的積分得到.水滴所受到的慣性力與作用在水滴上的力保持平衡，這其中包括阻力、浮力和其他相互作用力.方程(3)舉例給出了水滴在笛卡爾坐標下y方向受到的力的平衡方程.

(3)

式中：ρp為水滴粒子密度；vP為水滴下落速度；v為空氣速度；FD(v-vP)為單位水滴粒子質量所受到的阻力，阻力的定義如方程(4)所示；其他的相互作用力(Fother)是作為源項加在數值程序里面來表述其在特定區域的影響.

(4)

式中：CD為氣流作用于球形液滴的阻力系數，其為水滴雷諾數Re的函數；而雷諾數是液滴和空氣之間相對速度的函數，其分別定義為：

(5)

(6)

式中：DP為水滴粒子直徑.軌跡方程和其他描述液滴中的熱和質量傳輸的輔助方程可以通過在離散時間步上逐步積分得到.通過聯解方程(2)和方程(3)可得到給定時間步下液滴的速度和位置.

自然通風冷卻塔內部空氣與液滴之間的熱傳導包括對流和蒸發，而不考慮因為輻射而產生的熱交換，因為輻射產生的熱交換相對來說并不明顯. 通過求解一個水滴與周圍空氣的對流熱傳導方程來確定TP的值：

(7)

式中：cP為比熱容.將方程(7)在線性時間步Δt上進行積分，假設在一個時間步到下一個時間步的過程中液滴的溫度是緩慢改變的，因此得到：

(8)

式中：MP為液滴質量；AP為液滴質點面積；TP為液滴溫度；Tadb為空氣干球溫度；熱交換系數h可以通過Ranz和Marshal[11-12]提出的相互關系式來估算：

(9)

式中：Prma為濕空氣普朗特數；Nu為無量綱準則數.由于空氣流動而產生的水的蒸發率與水滴表面與空氣之間的蒸發壓力梯度成比例，則

(10)

式中：Psat為飽和空氣壓力；Pop為工作壓力.從熱和質量傳輸分析得到質量傳輸系數hm：

(11)

式中：Sh為無量綱準則數.液滴質量的減少可以通過以下方程得到：

(12)

最后通過將液滴與空氣相之間的對流和蒸發熱傳導與液滴中的熱改變之間的熱平衡關系式來校正液滴的溫度：

(13)

除了水滴軌跡、熱、質量和動量的獲取與損失需要計算外，這些量同樣需要作為附加源項(Spφ)耦合在相應空氣相的計算中.

(14)

式(14)中熱交換項將在接下來的空氣相流場計算之間的能量平衡中作為能量源出現.此外，質量交換將作為一個空氣相連續方程中的質量源項以及守恒方程中的水蒸汽的源項而出現在兩個方程中，通過方程(15)來計算.

(15)

動量交換在空氣相動量方程中作為一個源項或者匯項(Spu)出現，通過方程(16)來計算.

(16)

3 計算參數

幾何參數：受某電力設計院的委托進行本文的相關研究并將此研究作為技術儲備.本文所研究的冷卻塔塔高220m，為當前第一高塔.其主要參數為：淋水面積20 000m2，塔筒底部標高13.45m，塔筒底部直徑169.878m，喉部標高169.4m，喉部直徑103.545m，塔頂出口直徑109m.

計算域及邊界條件：計算域及邊界條件如圖2所示.由圖2可知，滿足阻塞比小于5%的要求.入口邊界條件為速度入口；出口采用壓力出口邊界條件，相對壓力選為零；流域側壁以及頂部采用自由滑移壁面條件；地面采用無滑移壁面邊界條件；冷卻塔表面為具有粗糙度的無滑移壁面邊界條件，粗糙度的施加參考文獻[13]中的方法；塔底只考慮填料區與雨區的影響，初始階段先在填料區頂部建立一個面施加熱水水滴，熱水經過填料區以及雨區過程中進行傳熱傳質交換，通過UDF編寫源項以及阻力項實現.

冷卻塔所在地為B類地貌，地面粗糙度系數α為0.16.入口處速度剖面和湍流度采用規范[8]中的指數率形式，該地區10m高度處50年重現期最大平均風速為25m/s.地面高度30m處湍流強度為16%.

湍流模型：本計算按風剖面輸入風速計算，喉部位置的雷諾數超過107，屬于完全發展的湍流階段；本研究雖然涉及到傳熱傳質計算，而對于風場的模擬依然是一個典型的由于逆壓梯度而引起分離的鈍體繞流問題；同時，本文討論的結果主要是內表面的平均壓力系數，對脈動信息要求低，因而采用模擬逆壓梯度引起流體分離具有較高精度的剪切應力輸運k-ω模型.

(a) 計算區域俯視圖

(b) 計算區域正視圖

計算工況及主要參數：工況參數由設計院通過實測提供.干球溫度308.47K，濕球溫度299.73K，大氣壓力為98.1kPa，進塔水溫為316.5K，出塔水溫為303.61～305.06K，進塔水質量流量為70 328.6t/h，水滴直徑為5mm，相對濕度為80%.

網格劃分：由于計算既包含傳熱傳質計算，同時也包含外流場側風計算，兩種計算的混合非常復雜，存在計算收斂較難的風險，因此全流域采用六面體網格，并對網格進行嚴格的細化處理，網格延伸率都為1.05，Yplus控制在2以下，本計算中的Yplus值符合SST湍流模型的要求，同時網格延伸率也比相關文獻[14]建議值小，遵循了嚴格的網格劃分方法，因此計算中忽略了網格無關性測試.本計算中在網格劃分中忽略人字柱，其引起的壓力損失通過施加源項來補償.網格如圖3所示，網格單元總數為5 738 890個.

4 計算結果及分析

計算包括兩個工況：工況1為無風時塔的傳熱傳質計算，通過與實測結果比較，以驗證本文研究方法的正確性；工況2為B類側風下塔在運行過程中的內壓計算研究.先計算純風場下的穩態結果，然后噴入熱水進行穩態計算收斂，將前面穩態結果作為瞬態計算的初場，最后進行瞬態計算.瞬態計算時間步為0.001s.計算在長沙理工大學的數值仿真研究中心的雙CPU、12核工作站上進行，總耗時約為18d.

本文主要討論的平均風壓系數定義為：

(17)

(a)塔局部網格圖

(b) 填料、雨區網格圖

(c) 整體網格圖

4.1 方法驗證工況計算

對無側風下運行塔進行計算，通過比較運行塔的出水溫度、塔內的流場特性以及內壓對稱性來驗證本文方法在運行工況下內壓計算的準確性.在該工況下計算得到的出塔平均水溫為305.8K，相比于實測出塔水溫303.61～305.06K吻合良好.圖4描述了靜風時刻塔內溫度等值線圖以及水池位置的溫度等值線圖.由圖4可知，冷卻塔左右兩側環境空氣呈現對稱形式流入塔內.環境空氣進入塔內后，由于塔自身形成的抽力作用和雨區水滴對空氣運動的阻力作用，沿徑向雨區空氣流速不斷減少，在塔雨區中心位置出現最小值.而由于雨區中心位置風速最小，塔內的中心位置的空氣溫度最高.

塔內溫度等值線基本呈對稱分布，在塔出口位置對稱性略差，究其主要原因是本文編寫的UDF程序并未考慮一些工藝上的細節，比如只考慮了雨區以及填料區，而在填料區中水流流動的水膜性質通過給定速度的液滴流來近似，從而造成一定的誤差.

(a) 中面溫度等值線圖

(b) 水池出塔水溫等值線圖

然而，由于本文主要研究對象是運行工況下塔內的壓力問題，圖5(a)中塔內壓力分布情況也顯示，塔中面壓力等值線圖對稱性良好，隨著高度越高壓力越大，而在填料區上下呈現出了良好的壓降情況；而圖5(b)中塔的內壓等值線圖也表明，沿著塔的高度方向，壓力均勻分布，且對稱性良好，高度越高壓力越大.塔出水的平均溫度為305.8K也與實測結果相符，塔內壓問題的分析相對于塔工藝的計算分析來說對于塔內的溫度以及速度的對稱性要求相對較小，所以，綜合各種指標表明本文計算方法滿足內壓分析的要求.

4.2 側風工況下運行塔內壓計算結果

圖6為B類風場下塔運行時的內壓系數特征圖，h為冷卻塔高度，z為水平截面高度.圖中顯示，各高度內壓曲線基本上沿緯向角分布均勻，而且壓力系數隨著高度的增加而緩慢變大，表明在側風的影響下，施加在雨區和填料區阻力和壓降對內壓起到明顯的整流作用.不同高度內壓系數平均值為-0.48～-0.41，峰值為-0.51～-0.403.整塔內壓系數平均值為-0.449，喉部位置內壓系數平均值為-0.43.

(a) 中面壓力等值線圖

(b) 整塔內壓等值線圖

角度/(°)

4.3 側風下塔內流場特性分析

圖7描述了B類風場下運行時塔內流動特性.圖7(a)表明塔內溫度場保持相對均勻，由于側風的影響，流場的對稱性被打破，上游入口位置溫度明顯比下游入口溫度低，塔底位置背風區附近出現了一個高溫區，而塔外的背風區也有一個相對于外圍空氣溫度的高溫區.塔底的速度(圖7(b))也由于側風的影響而不均勻，上游側風直接穿過雨區以及填料區而在塔內迎風面形成一個漩渦，而雨區的背風位置出現了一個低速區，導致此處的溫度升高而嚴重影響塔的冷卻效果.塔內壓力圖(圖7(c))顯示，由于塔內抽力的原因，塔內壓力沿著高度方向相應變大，塔內整體壓力分布雖然由于側風的影響而出現不均勻，但是經過填料區和雨區的整流，塔內的壓力分布從數值上沿著緯向角變化并不明顯.

(a) 中面溫度等值線圖

(b) 中面速度等值線圖

(c)中面壓力等值線圖

4.4 規范內壓系數取值的探討

目前最高的冷卻塔在德國，高度為200 m，與文中計算塔高220 m基本上處于同一量級，德國規范規定內壓系數取常值-0.5，本文計算結果最小峰值為-0.51，各截面高度平均值最小為-0.48，兩者取值比較接近.從壓力損失的角度分析，由于塔為開口結構，并且塔內為自然對流運動，內部抽力應該在越接近開口的地方抽力越小，越靠近底部抽力相應變大，中間應該是緩慢過渡，本文計算結果符合此規律，德國規范針對塔所有高度取同一值相對偏安全.文獻[4]建議內壓系數取值為-0.4～-0.5，而且設計時采用-0.45，略大于本文結果，對于本文220 m的超高塔來說，局部內壓取值稍微偏不安全.因此，從安全的角度考慮，建議對于220 m高度級別的冷卻塔風荷載的內壓系數選取參考德國規范.

5 結 論

1)本文以CFD軟件為基礎進行二次開發，得到了冷卻塔真實運行過程中的傳熱傳質過程，還原了塔內真實流動特征，提出了一種計算冷卻塔運行過程中的內壓系數計算方法，并通過無側風時計算結果與實測結果的比較，驗證了本文方法的準確性.

2)側風下，施加在雨區和填料區的阻力和壓降對內壓起到明顯的整流作用，各截面高度內壓系數曲線沿緯向角分布均勻，而且由于塔內自然對流抽力的原因，壓力系數隨著高度的增加而緩慢變大.

3)在B類風場下塔運行時的整塔內壓系數平均值計算結果為-0.449，喉部位置內壓系數平均值為-0.43，不同高度內壓系數平均值為-0.41～-0.48.因此，從安全的角度出發，建議對于220 m高度的級別的冷卻塔風荷載的內壓系數選取參考德國規范.

4)本文結論是基于高度為220 m量級的逆流式自然對流冷卻塔計算分析得到，塔的內壓系數取值大小與塔高度之間的關系還需要做進一步深入研究.

[1] 李鵬飛,趙林,葛耀君,等.超大型冷卻塔風荷載特性風洞試驗研究[J].工程力學,2008,25(6):60-67.

LI Peng-fei, ZHAO Lin, GE Yao-jun,etal. Investigation on wind load characteristics for super large cooling tower in wind tunnel[J]. Engineering Mechanics,2008,25 (6):60-67. (In Chinese)

[2] 孫天風,周良茂.無肋雙曲線型冷卻塔風壓分布的全尺寸測量和風洞研究[J].空氣動力學報,1983,4:68-76.

SUN Tian-feng, ZHOU Liang-mao. A full scale and wind tunnel study of wind pressure distribution around a ribless hyperbolic cooling tower[J]. Acta Aerodynamica Sinica,1983,4:68-76. (In Chinese)

[3] 鄒云峰，牛華偉，陳政清.基于完全氣動彈性模型的冷卻塔干擾效應風洞試驗研究[J].湖南大學學報:自然科學版,2013,40(12):1-7.

ZOU Yun-feng, NIU Hua-wei, CHEN Zheng-qing. Wind tunnel test on wind-induced interference effect of cooling tower based on full aero-elastic model[J]. Jouranl of Hunan University：Natural Sciences, 2013,40(12):1-7. (In Chinese)

[4] 柯世堂，侯憲安，姚友成，等. 強風作用下大型雙曲冷卻塔風致振動參數分析[J].湖南大學學報:自然科學版,2013,40(10):32-37.

KE Shi-tang, HOU Xian-an, YAO You-cheng,etal. Parameter analysis of wind-induced vibration for large hyperbolic cooling towers under strong wind loads[J]. Jouranl of Hunan University：Natural Sciences, 2013,40(10):32-37. (In Chinese)

[5] DAMJAKOB H, TUMMERS N.Back to future of the hyperbolic concrete tower[C]//Proceedings of the 5th International Symposium on Natural Draught Cooling Towers.London:A.A.Balema Publishers,2004:3-21.

[6] KAWARABATA Y,NAKAE S,HARADA N. Some aspects of the wind design of cooling towers[J]. Journal of Wind Engineering and Industrial Aerodynamics, 1983,14:167-180.

[7] VGB-Guideline:Structural design of cooling tower-technical guideline for the structural design, computation and execution of cooling tower(VGB-R610Ue)[S]．Essen:BTR Bautechnik Bei Kühltürmen, 2005.

[8] 中華人民共和國建設部. GB/T 50102-2003 工業循環水冷卻設計規范[S].北京:中國電力出版社,2003.

Ministry of Construction of the People's Republic of China. GB/T 50102-2003 Code for design of cooling for industrial recirculating water[S].Beijing: China Electric Power Press, 2003. (In Chinese)

[9] 中華人民共和國國家發展和改革委員會. NDGJ5-88 火力發電廠水工設計技術規定[S].北京:中國電力出版社,1989.

The People's Republic of China National Development and Reform Commission. NDGJ5-88 Technical regulation for hydraulic design of fossil fuel power plant[S]. Beijing: China Electric Power Press, 1989.(In Chinese)

[10]AL-WAKED R, BEHNIA M. CFD simulation of wet cooling towers[J]. Applied Thermal Engineering,2006,26:382-395.

[11]RANZ W,MARSHALL W. Evaporation from drops:part I[J]. Chemical Engineering Progress,1952,48: 141-146.

[12]RANZ W,MARSHALL W. Evaporation from drops:part II[J]. Chemical Engineering Progress,1952,48:173-180.

[13]董國朝,陳政清,羅建輝,等.冷卻塔混凝土粗糙度對平均風壓系數的影響[J].湖南大學學報:自然科學版,2011,38(7):6-12.

DONG Guo-chao, CHEN Zheng-qing, LUO Jian-hui,etal. Effect of surface roughness on the mesan pressure coefficient of concrete cooling tower[J]. Jouranl of Hunan University：Natural Sciences, 2011,38(7):6-12. (In Chinese)

[14]HUANG Sheng-hong,LI Q S, XU Sheng-li. Numerical evolauation of wind effects on a tall steel building by CFD[J]. Journal of Constructional Steel Research,2007, 63:612-627.

Numerical Simulation of the Internal Surface Wind Load of Super Large Cooling Tower under Operating Conditions

DONG Guo-chao1?, ZHANG Jian-ren1, CAI Chun-sheng1,2, HAN Yan1,LI Chun-guang1

(1. School of Civil Engineering and Architecture,Changsha Univ of Science and Technology, Changsha, Hunan 410114,China；2. Dept of Civil and Environmental Engineering, Louisiana State Univ,Baton Rouge, USA LA 70803)

The average wind pressure on the internal surface of a super large hyperbolic cooling tower under the operating conditions was simulated in CFD method. Based on computational fluid dynamics software for secondary development, DPM model combined with a UDF function method was used to study the average wind pressure on the internal surface of a super large cooling tower. The Lagrangian method was used to simulate the water phase of the tower and Euler method was adopted to simulate the air phase, the coupling calculation between inner flow field and transformation of heat and mass in cooling tower under the operating conditions was well realized, and the internal pressure distribution law of cooling tower under crosswind direction with running state was analyzed. The results of the tower with no cross wind show that the symmetry of the average pressure on the internal surface goes well, and the water temperature is consistent with the test result, which verifies the effectiveness of the method proposed. The value of the average pressure coefficient of the internal surface of cooling tower under cross wind becomes bigger along the height direction, and the value does not change significantly along the latitudinal direction. The shortcomings of the current code about internal surface pressure coefficient were discussed, and the recommended values were given, which provides methods and basis for the calculation of the internal pressure of super large cooling tower design.

cooling tower; average wind load; numerical simulation; heat and mass transformation; DPM(Discrete Phase Model)

1674-2974(2015)01-0017-07

2014-06-11

國家自然科學基金青年基金資助項目(51408061), Youth Fund of National Natural Science Foundation Projects (51408061); 國家重點基礎研究計劃(973計劃)資助項目(2015CB057700)

董國朝(1981-)，男，廣東陽春人，長沙理工大學副教授, 博士?通訊聯系人，E-mail: dgccpu@163.com

TU375.4

A