董国朝 ，张建仁，蔡春声,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