运行状态下超大型冷却塔内表面风荷载的数值模拟研究

董国朝等

摘要：对运行状态下的某超大型双曲冷却塔的内表面平均风压进行了CFD数值模拟.在计算流体动力学软件基础上进行二次开发，采用DPM模型结合UDF函数方法加入源项来研究某超大型冷却塔内表面平均风压分布；塔中水相采用了拉格朗日方法模拟，而空气相采用欧拉方法模拟，较好地实现了冷却塔运行状态下的内外流场计算及其与传热传质的耦合计算，分析了运行状态下冷却塔横风向来流时的内压分布规律.无侧风工况下计算结果显示，塔运行过程中的内表面压力对称性良好，出水温度与实测结果相符，验证了本文提出的塔运行过程中的传热传质计算方法的正确性.侧风工况下得到塔内压力系数沿高度方向相应变大，而沿纬向变化不明显.同时讨论了中国规范对内表面压力系数的取值不完善之处，给出了建议取值，为超大型冷却塔设计过程中的内压计算提供方法和依据.

关键词：冷却塔；平均风荷载；数值模拟；传热传质；离散相模型（DPM）

中图分类号：TU375.4 文献标识码：A

超大型双曲自然通风冷却塔是一种大型空间薄壁开口结构，普遍用于发电厂中循环水冷却，风荷载对其安全性的影响起决定性作用.随着电站装机容量的增加，冷却塔也向高大化发展.如湖南、湖北拟建和在建的塔高分别达到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中，描述热、质量和动量传输的空气（连续相）流动方程可以写成以下统一形式：

3计算参数

几何参数：受某电力设计院的委托进行本文的相关研究并将此研究作为技术储备.本文所研究的冷却塔塔高220 m，为当前第一高塔.其主要参数为：淋水面积20 000 m2，塔筒底部标高13.45 m，塔筒底部直径169.878 m，喉部标高169.4 m，喉部直径103.545 m，塔顶出口直径109 m.

计算域及边界条件：计算域及边界条件如图2所示.由图2可知，满足阻塞比小于5%的要求.入口边界条件为速度入口；出口采用压力出口边界条件，相对压力选为零；流域侧壁以及顶部采用自由滑移壁面条件；地面采用无滑移壁面边界条件；冷却塔表面为具有粗糙度的无滑移壁面边界条件，粗糙度的施加参考文献[13]中的方法；塔底只考虑填料区与雨区的影响，初始阶段先在填料区顶部建立一个面施加热水水滴，热水经过填料区以及雨区过程中进行传热传质交换，通过UDF编写源项以及阻力项实现.

冷却塔所在地为B类地貌，地面粗糙度系数α为0.16.入口处速度剖面和湍流度采用规范[8]中的指数率形式，该地区10 m高度处50年重现期最大平均风速为25 m/s.地面高度30 m处湍流强度为16%.

湍流模型：本计算按风剖面输入风速计算，喉部位置的雷诺数超过107，属于完全发展的湍流阶段；本研究虽然涉及到传热传质计算，而对于风场的模拟依然是一个典型的由于逆压梯度而引起分离的钝体绕流问题；同时，本文讨论的结果主要是内表面的平均压力系数，对脉动信息要求低，因而采用模拟逆压梯度引起流体分离具有较高精度的剪切应力输运kω模型.

计算工况及主要参数：工况参数由设计院通过实测提供.干球温度308.47 K，湿球温度299.73 K，大气压力为98.1 kPa，进塔水温为316.5 K，出塔水温为303.61～305.06 K，进塔水质量流量为70 328.6 t/h，水滴直径为5 mm，相对湿度为80%.

网格划分：由于计算既包含传热传质计算，同时也包含外流场侧风计算，两种计算的混合非常复杂，存在计算收敛较难的风险，因此全流域采用六面体网格，并对网格进行严格的细化处理，网格延伸率都为1.05，Yplus控制在2以下，本计算中的Yplus值符合SST湍流模型的要求，同时网格延伸率也比相关文献[14]建议值小，遵循了严格的网格划分方法，因此计算中忽略了网格无关性测试.本计算中在网格划分中忽略人字柱，其引起的压力损失通过施加源项来补偿.网格如图3所示，网格单元总数为5 738 890个.

4计算结果及分析

计算包括两个工况：工况1为无风时塔的传热传质计算，通过与实测结果比较，以验证本文研究方法的正确性；工况2为B类侧风下塔在运行过程中的内压计算研究.先计算纯风场下的稳态结果，然后喷入热水进行稳态计算收敛，将前面稳态结果作为瞬态计算的初场，最后进行瞬态计算.瞬态计算时间步为0.001 s.计算在长沙理工大学的数值仿真研究中心的双CPU、12核工作站上进行，总耗时约为18 d.

4.1方法验证工况计算

对无侧风下运行塔进行计算，通过比较运行塔的出水温度、塔内的流场特性以及内压对称性来验证本文方法在运行工况下内压计算的准确性.在该工况下计算得到的出塔平均水温为305.8 K，相比于实测出塔水温303.61～305.06 K吻合良好.图4描述了静风时刻塔内温度等值线图以及水池位置的温度等值线图.由图4可知，冷却塔左右两侧环境空气呈现对称形式流入塔内.环境空气进入塔内后，由于塔自身形成的抽力作用和雨区水滴对空气运动的阻力作用，沿径向雨区空气流速不断减少，在塔雨区中心位置出现最小值.而由于雨区中心位置风速最小，塔内的中心位置的空气温度最高.

塔内温度等值线基本呈对称分布，在塔出口位置对称性略差，究其主要原因是本文编写的UDF程序并未考虑一些工艺上的细节，比如只考虑了雨区以及填料区，而在填料区中水流流动的水膜性质通过给定速度的液滴流来近似，从而造成一定的误差.

然而，由于本文主要研究对象是运行工况下塔内的压力问题，图5（a）中塔内压力分布情况也显示，塔中面压力等值线图对称性良好，随着高度越高压力越大，而在填料区上下呈现出了良好的压降情况；而图5（b）中塔的内压等值线图也表明，沿着塔的高度方向，压力均匀分布，且对称性良好，高度越高压力越大.塔出水的平均温度为305.8 K也与实测结果相符，塔内压问题的分析相对于塔工艺的计算分析来说对于塔内的温度以及速度的对称性要求相对较小，所以，综合各种指标表明本文计算方法满足内压分析的要求.

4.2侧风工况下运行塔内压计算结果

图6为B类风场下塔运行时的内压系数特征图，h为冷却塔高度，z为水平截面高度.图中显示，各高度内压曲线基本上沿纬向角分布均匀，而且压力系数随着高度的增加而缓慢变大，表明在侧风的影响下，施加在雨区和填料区阻力和压降对内压起到明显的整流作用.不同高度内压系数平均值为-0.48～-0.41，峰值为-0.51～-0.403.整塔内压系数平均值为-0.449，喉部位置内压系数平均值为-0.43.

4.3侧风下塔内流场特性分析

图7描述了B类风场下运行时塔内流动特性.图7（a）表明塔内温度场保持相对均匀，由于侧风的影响，流场的对称性被打破，上游入口位置温度明显比下游入口温度低，塔底位置背风区附近出现了一个高温区，而塔外的背风区也有一个相对于外围空气温度的高温区.塔底的速度（图7（b））也由于侧风的影响而不均匀，上游侧风直接穿过雨区以及填料区而在塔内迎风面形成一个漩涡，而雨区的背风位置出现了一个低速区，导致此处的温度升高而严重影响塔的冷却效果.塔内压力图（图7（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 Pengfei， ZHAO Lin， GE Yaojun，et al. 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 Tianfeng， ZHOU Liangmao. 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 Yunfeng， NIU Huawei， CHEN Zhengqing. Wind tunnel test on windinduced interference effect of cooling tower based on full aeroelastic model[J]. Jouranl of Hunan University：Natural Sciences， 2013，40（12）：1-7. （In Chinese）

[4]柯世堂，侯宪安，姚友成，等. 强风作用下大型双曲冷却塔风致振动参数分析[J].湖南大学学报：自然科学版，2013，40（10）：32-37.

KE Shitang， HOU Xianan， YAO Youcheng，et al. Parameter analysis of windinduced 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]VGBGuideline：Structural design of cooling towertechnical guideline for the structural design， computation and execution of cooling tower（VGBR610Ue）[S].Essen：BTR Bautechnik Bei Kühltürmen， 2005.

[8]中华人民共和国建设部. GB/T 501022003工业循环水冷却设计规范[S].北京：中国电力出版社，2003.

Ministry of Construction of the People's Republic of China. GB/T 501022003Code for design of cooling for industrial recirculating water[S].Beijing： China Electric Power Press， 2003. （In Chinese）

[9]中华人民共和国国家发展和改革委员会. NDGJ588火力发电厂水工设计技术规定[S].北京：中国电力出版社，1989.

The People's Republic of China National Development and Reform Commission. NDGJ588Technical regulation for hydraulic design of fossil fuel power plant[S]. Beijing： China Electric Power Press， 1989.（In Chinese）

[10]ALWAKED 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 Guochao， CHEN Zhengqing， LUO Jianhui，et al. 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 Shenghong，LI Q S， XU Shengli. Numerical evolauation of wind effects on a tall steel building by CFD[J]. Journal of Constructional Steel Research，2007， 63：612-627.