超大型冷却塔风-雨双向耦合作用机理和气动力分布研究

余文林, 柯世堂

(1.南京航空航天大学 土木工程系，南京 210016；2.中国能源建设集团江苏省电力设计院有限公司，南京 211102)

相比常规冷却塔，超大型冷却塔表面气动力分布与流场作用机理更加复杂[1-2]。尤其在极端气候条件下，结构要同时承受强风和暴雨的共同作用，此时雨滴在风力和重力的共同驱动下运动轨迹发生倾斜，以较大速度击打在冷却塔壁面，使其表面气动力分布发生显著改变。同时，暴雨会极大地恶化无雨状态下脉动风的湍流效应，此时塔筒近壁面气流运动形式复杂，进而改变雨滴轨迹、附加作用力和表面压力作用。鉴于此，对于复杂风雨联合作用下的超大型冷却塔，研究其在不同风速和雨强组合下的塔筒表面压力作用机理和气动力分布，具有重要的理论意义和工程价值。

对于风雨荷载共同作用的研究，文献[3-4]针对建筑物在风雨荷载作用下的风驱雨量以及风雨荷载进行了定性和定量化的分析；文献[5-6]采用数值模拟手段分析了稳态风驱动下低矮房屋表面雨荷载效应以及与对应区域风荷载的比值；文献[7-8]基于数值模拟和有限元结合的方法对大型风力机风雨荷载特性及受力性能展开了详细研究；文献[9-10]采用有限元分析与理论推导的方法提出了一种输电塔体系的风雨激励动力分析模型；文献[11-13]基于风洞试验和数值模拟的方法进行了斜拉桥拉索风雨理论模型、激振机理及参数影响的系统研究。已有关于风驱雨的研究多集中在低矮房屋、斜拉索、输电塔及风力机等结构，且仅仅考虑了风单项驱动雨的影响，忽略了雨对于风的反作用。而位于气候条件较差的超大型冷却塔结构，是否需要考虑风雨耦合作用对其塔筒表面气动性能的影响，目前鲜有文献报道。

为解决该问题，以国内某建成210 m世界最高超大型冷却塔为例，基于风-雨双向耦合算法为核心，采用CFD手段首先对不同风速下的冷却塔周围风场进行数值模拟，在稳定后的风场中添加离散相模型(DPM)以输入不同等级的雨量，然后进行风雨场的同步迭代计算。在此基础上，对比研究风雨联合作用下塔筒表面压力作用机理，提炼出不同风速和雨强对超大型冷却塔表面风驱雨量、雨滴附加荷载及等效压力系数的影响规律。最终，提炼出最不利风雨组合工况，并基于非线性最小二乘法原理拟合给出了超大型冷却塔等效压力系数的计算公式和对应的二维空间曲面。

1 风-雨双向耦合算法

1.1 降雨强度

表1给出了降雨强度分类列表，两种雨强分类的采样时间不同，同一场雨的测量结果差别较大，采用小时雨强进行结构验算相对日雨强偏于安全，小时降雨量更能直观反映出工程中最为关心的极端气候条件下瞬时雨强对结构性能的影响。

表1 降雨强度等级划分

1.2 雨滴谱分布

雨滴谱近似服从负指数分布，常用模型[14]有Best谱、Marshall-Palmer谱及Gamma谱。本文选用Marshall-Palmer谱，表达式为

n(D)=N0e-λD

(1)

式中：D为雨滴直径，单位为mm；n(D)为不同直径雨滴个数浓度谱；N0为浓度，取常数值8 000；λ为尺度参数，其表达式为

λ=4.1×R-0.21

(2)

1.3 雨滴末速度

雨滴受重力作用下降速度不断增加，空气阻力随之加大，当重力与阻力平衡时，雨滴以最终速度匀速下落，该速度称为雨滴末速度或极限速度。

文献[15]指出降落高度≥20 m即可使几乎所有粒径雨滴达到最终末速度。文献[16]认为小直径雨滴在降落过程中可忽略受力变形，从而近似为球体下落。

1.4 壁面碰撞方程

雨滴冲击到塔筒壁面过程服从动量守恒定律，计算中忽略雨滴在冲击过程中可能发生的蒸发、飞溅、破裂等现象，认为雨滴与结构间相互作用遵循牛顿第二定律。由动量定理：

(3)

式中：f(t)为单个雨滴冲击力矢量，单位为N；v为雨滴速度矢量。

雨滴在单位时间内对结构的冲击力F(τ)为：

(4)

将下落时雨滴近似看作球体，则：

(5)

雨滴直径一般小于6 mm，且撞击前水平末速度相对较大，故为简化计算，将碰撞时间τ取为：

(6)

则雨滴对结构的冲击力可简化为：

(7)

2 工程简介与工况设置

2.1 工程简介

该建成超大钢筋混凝土双曲线自然通风间接空冷塔，塔高210 m，喉部标高157.5 m，进风口标高32.5 m，喉部内面直径110 m，进风口内面直径159 m，零米直径为180 m。冷却塔塔筒采用指数变厚，最小厚度位于喉部断面，壁厚为0.37 m，最大壁厚在下环梁位置，壁厚为2.0 m。塔筒由52对X型支柱支撑且与环板基础连接，X型柱采用矩形截面，尺寸为1.8 m×1.2 m。环板基础为现浇钢筋混凝土结构，宽为12.0 m，高为2.5 m。表2给出了该工程冷却塔的主要结构尺寸及示意图。

表2 超大型冷却塔主要结构尺寸及示意图

2.2 工况设置

该冷却塔位于B类地貌，考虑到冷却塔百叶窗的常规工作状态，按30%透风率考虑百叶窗开启效应[17]，对比研究风雨联合作用下3种风速和3种雨强组合对冷却塔表面气动性能的影响。其中，小风、中风和大风分别以重现期为10年、50年和100年最大风速为进行划分；雨强均以大暴雨气候条件的弱大暴雨、中大暴雨和强大暴雨为基准进行分类，共9种对比工况(见图1)。

图1 对比工况组合示意图

3 风-雨双向耦合数值模拟

3.1 建立风雨场模型

为同时保证冷却塔处于降雨区和尾流充分发展，计算域尺寸设置为顺风向长3 000 m，横风向宽1 500 m，高度方向为600 m。选取塔底中心为坐标系原点，x轴与顺风向一致。为了兼顾计算效率和精度，划分网格时将其分为局部和外围风雨场，局部风雨场内含冷却塔模型，采用非结构化网格进行划分，外围风雨场形状规整，采用高质量的结构化网格进行划分。核心区最小网格尺寸为0.2 m，整体模型总网格数量超过1 800万，网格最小正交质量大于0.4(要求大于0.1且杜绝出现负体积[18])，网格数目及质量均满足计算要求。图2给出了整体计算域和模型网格划分示意图。

(a)整体网格划分(b)局部网格划分

图2 整体与局部网格划示意图

Fig.2 Grid division of total and local models

设置计算域入口为速度入口(velocity inlet)边界，出口为压力出口(pressure oulet)边界，两侧壁及顶面采用对称(symmetry)边界，冷却塔及地面均设为壁面(wall)边界，局部与外围计算域的重合面设为交界面(interface)。风雨场计算域及其边界条件，如图3所示。

图3 计算区域与边界条件示意图

3.2 风-雨场耦合计算

数值计算采用3D单精度、分离式求解器，流场流速为绝对速度，空气模型等效为理想不可压缩流体，湍流模型选取k-ω剪切应力(SST)控制方程。计算域入口采用幂指数为0.15的风廓线模型，离地高度10 m高度处的风速分别设置为2.2节中三种基准风速。流场求解采用SIMPLEC算法实现速度与压力之间的耦合，对流项求解格式为二阶，计算过程中设置了网格倾斜校正以提高混合网格计算效果，控制方程的计算残差设置为1×10-6，最后初始化风场进行迭代计算。图4给出了平均风速、湍流度剖面模拟结果与理论值的对比曲线，结果表明平均风速和湍流度剖面均与理论值吻合良好，风场模拟标准满足工程要求。

风场求解稳定后插入离散相继续进行风雨场耦合迭代运算。采用1.0～6.0 mm范围内6种直径的雨滴来模拟连续直径分布的降雨(见表3)，每种直径雨滴数量和体积占有率采用1.2节中Marshall-Palmer谱确定。对雨滴进行“面”释放，水平释放速度为0，竖向释放速度为-5 m/s，重力和阻力共同作用使得雨滴在足够大的高度范围内达到末速度。风雨耦合迭代完成后，可输出连续相流场结果及冷却塔表面捕捉到的雨滴信息，据此计算雨滴对冷却塔表面的撞击作用及附加荷载。

(a) v0=20 m/s

(b) v0=23.7 m/s

(c) v0=25.3 m/s

3.3 有效性验证

对3种风速下冷却塔喉部区域平均风压系数与国内外相关规范[19-20]及实测结果[21]进行对比，结果如图5所示。分析可知，3种风速下冷却塔喉部区域平均风压分布曲线的负压极值点和分离点对应角度与国内外规范及西热曲线完全一致；迎风和背风区域风压系数吻合较好，基本包络国内外规范及实测曲线，仅侧风区负压值略有差异；综上认为本文数值模拟结果具有一定的有效性。

4 结果对比分析

4.1 风场分析

图6和图7分别给出了加入雨滴前三种风速下冷却塔三维风速流线和涡量分布，由图可看出：

图5 3种风速下冷却塔喉部区域压力系数与国内外规范及实测结果对比曲线

Fig.5 Comparison among pressure coefficient in the throat area of the cooling tower under three wind speeds, domestic and foreign codes, actual measured results

①来流流经塔筒在迎风面产生分流，进而沿塔筒两侧外壁加速流动至背风面形成尺寸大小不同的涡旋，部分气流透过百叶窗进入塔筒内部，在塔体内表面附着流动、撞击并向上爬升，由于双曲线型冷却塔在喉部位置的颈缩，上升气流受到一定程度的阻碍，同时塔顶位置来流加速掠过改变了部分上升气流方向，气流产生回流并在喉部形成较大尺度的涡旋；②随着风速的增大，塔筒气流加速运转，涡旋脱落现象更为显著，风速流线更为密集，上升气流与来流在塔顶作用加剧，因此在背风面涡旋脱落更显著；③湍动能强度随着风速增大而增强，在塔筒背风区中下部尤为显著，峰值主要位于塔筒迎风区顶部、背风区喉部、百叶窗迎风区内部及背风区外部。

(a) v0=20 m/s

(b) v0=23.7 m/s

(c) v0=25.3 m/s

(a) v0=20 m/s

(b) v0=23.7 m/s

(c) v0=25.3 m/s

4.2 雨场分析

基于颗粒合速度对雨滴轨迹进行追踪，图8给出了9种工况风雨场中雨滴运动轨迹示意图，并对雨滴的密集程度进行了等比例粗化处理。由图可以看出：①受风力驱动作用，雨滴改变其竖直运动方向，并以一定的斜向速度降落，其运动轨迹倾斜率随风速的增大更为显著，而受降雨强度影响微弱；②冷却塔近壁面雨滴运动形式复杂，部分雨滴在风力作用下击打至冷却塔表面，其余雨滴沿气流在结构两侧发生分离，少量附着在冷却塔侧壁面，多数随风进入尾流区域；③冷却塔壁面收集的雨滴随雨强的增强而增多，随风速的增加而减少。

(a) 工况1

(b) 工况2

(c) 工况3

(d) 工况4

(e) 工况5

(f) 工况6

(g) 工况7

(h) 工况8

(i) 工况9

图9给出了9种工况下塔筒壁面收集到的各直径雨滴数量、撞击速度及速度占有率对比曲线。由图可知：

1)各工况塔筒表面收集到的雨滴直径均主要分布在3～6 mm之间，且以5 mm直径雨滴占比最大，这是由于在同等风力驱动下，小直径雨滴速度增幅较快，竖向位置未及冷却塔高度范围时，水平方向已随风略过塔体进入尾流区。

2)结构表面雨滴捕捉数量随风速的增大而减少，随雨强的增大而增加，但各工况塔筒表面收集到的雨滴水平速度占有率分布规律基本一致，速度占有率随着水平末速度的增大先增加后减小，其中5～15 m/s范围的的速度占有率最大，为60%左右。

3)各工况塔筒表面收集到的雨滴平均水平速度均远小于基准风速，塔筒表面雨滴水平速度主要分布在3～10 m/s范围内，且雨滴水平撞击速度大体上均随着雨滴直径的增加而增大。

(a) 雨滴数量

(b) 速度占有率

(c) 水平末速度

图10给出了塔筒不同高度范围雨荷载特征值(图中箭头指向为雨荷载沿此方向增大)。由图可知：①各工况下冷却塔表面雨荷载整体上随着高度的增加先减小后增加，均在0.23～0.31H高度范围内最小，在塔底或塔顶处最大；②各风速下塔筒雨荷载随雨强的提升而显著增大，固定雨强下风速的提高将导致表面中下部雨荷载减小，但塔筒中上部雨荷载基本呈先增大后减小的趋势，且差别较小，主要原因为风速的提高会显著减少塔筒表面雨滴数量，在塔筒中下部雨滴数量的影响大于风速的影响，但在塔筒中上部雨滴速度的增加会显著增大雨荷载。

图10 各工况塔筒表面不同高度范围雨荷载示意图

为更清晰展示各高度和角度雨滴位置、数量和对应的压力系数，图11给出了各工况冷却塔表面雨滴以及雨致压力系数三维分布示意图，为清晰显示雨滴撞击位置，将风雨场坐标系逆时针旋转90°并对雨滴进行粗化处理。由图可知：

1)各工况雨滴撞击位置多集中分布在冷却塔表面迎风区域，受气流漩涡驱动作用，侧风区壁面和背风区壁面有少量雨滴附着，且塔筒表面收集到的雨量随风速的增加迅速减少，随雨强的增大逐渐变多，其中以工况3最多；

2)各工况雨致压力系数均主要集中于迎风面两侧各60°范围内，其余范围数值基本为0，雨致压力系数最大值为0.184，发生在工况3的0.15～0.23H高度范围内，且除了工况2和3的雨致压力系数最大值发生在塔底处，其它工况的雨致压力系数最大值均发生在塔顶位置。

(a) 工况1

(b) 工况2

(c)工况3

(d) 工况4

(e) 工况5

(f) 工况6

(g) 工况7

(h) 工况8

(i) 工况9

4.3 等效压力系数分析

为定量比较不同工况组合下的塔筒表面压力分布，定义等效压力系数，见式(8)～(10)，其中：Cpei为风雨耦合作用下冷却塔第i个监控点等效压力系数，Cpwi为监控点风致压力系数，Cpri为雨致压力系数；Pri为雨压，Pwz0为参考高度处风压，本文参考高度为塔顶210 m；Fri为雨荷载，Si为计算面积。

Cpei=Cpwi+Cpri

(8)

(9)

(10)

图12给出了不同工况下塔筒表面6个典型断面的等效压力系数对比曲线，由图分析可得：

1)各工况相同高度截面处等效压力系数分布规律及数值基本一致，仅迎风区与背风区数值略有差异，其中压力系数从迎风面到背风面均呈现先减小后增大再减小直至平稳的分布规律；

(a) 0.27 H(56 m)

(b) 0.51 H(104 m)

(c) 0.73 H(152 m)

(d) 0.80 H(168 m)

(e) 0.87 H(184 m)

(f) 0.95 H(200 m)

2)每个工况6个典型断面等效压力系数数值均略有差异，但均关于风轴基本呈现较好的对称性，其中最大负压数值整体上随着高度的增加先增大后减小，在喉部区域附近约为-1.5，且背风区域负压数值也基本呈先增大后减小的趋势，且背风面负压平稳区域逐渐变宽。

图13给出了塔筒0°、75°、120°及180°四条典型子午线的等效压力系数对比曲线，分析可得：

1)各工况相同角度子午向压力系数分布规律及数值基本一致，仅迎风区和背风区数值略有差异，每个工况不同角度等效压力系数分布规律差异显著；

2)各工况塔筒0°子午线等效压力系数基本为1.0左右，随着高度的增加数值略有差异；75°和180°子午线等效压力系数随着高度的增加基本呈现先增大后减小的趋势；120°子午线等效压力系数随着高度的增加基本不变，但在塔顶区域处数值减小。

4.4 等效压力系数二维拟合曲面

为方便工程设计人员精确获得此类超大型冷却塔等效压力系数，本文基于非线性最小二乘法原理，以子午向高度和环向角度为目标函数，拟合给出超大型冷却塔等效压力系数的计算公式，其中塔筒沿环向均分为N1段，沿子午向均分为N2段，令N=N1×N2，公式具体定义为：

Mz,θ=(ai×I+bi×Z+ci×Z·2+di×θ·×Z+ei×θ·3+fi×Z·3+gi×θ·4+hi×Z·4+ii×θ·×Z·3+ji×θ·2·×Z·2+ki×θ·5+li×Z·5)·÷(I+mi×exp(ni×θ+oi×Z))

(11)

式中：I为元素全为1的N×1矩阵，θ为以N1个角度为循环单位且循环N2次的N×1矩阵，Z为以每N2个相同的高度为循环单位且循环N1次的N×1矩阵，·×为矩阵对应元素相乘，·÷为矩阵对应元素相除，·n为矩阵对应元素的n次方，exp()为返回括号内矩阵每个元素作为以e为底的指数的矩阵，MZ,θ表示以N1个环向角度对应的压力系数为单位且沿子午向高度变化N2次的N×1矩阵，ai、bi，…，oi(i=1，2,…，9)为不同工况下的拟合系数。表4给出了最不利工况下(工况3)冷却塔等效压力系数拟合公式系数。

(a) 0°

(b) 75°

(c) 120°

(d) 180°

系数数值系数数值系数数值a3-0.317f3-9.56×10-6k31.22×10-11b3-0.047g3-6.91×10-9l3-9.60×10-11c30.001h35.08×10-8m3672.758d35.85×10-5i3-1.83×10-10n3-0.137e39.74×10-7j3-4.51×10-10o30.016

图14给出了最不利工况下冷却塔等效压力系数二维实际及拟合曲面对比图，图中散点数值为冷却塔实际等效压力系数，曲面对应数值为根据二维拟合公式模拟得到的等效压力系数。由图可知冷却塔等效压力系数沿子午向及环向存在明显二维特征，从拟合曲面的整体分布来看，其沿子午向和环向的变化规律与实际等效压力系数分布基本一致，仅在负压最大值区域略有差异，但误差率在10%以内。对比结果表明本文提出的二维拟合公式可作为此类超大型冷却塔表面压力系数的取值依据。

图14 最不利工况下冷却塔等效压力系数二维实际及拟合曲面对比图

Fig.14 Bivariate distribution and fitting surface of equivalent pressure coefficient under the most unfavorable condition

5 结 论

(1) 雨滴水平方向作用力随风速提升而增大，推动雨滴加速运行，受冷却塔前缘和内部流出气体的驱动和包裹作用，大量雨滴掠过塔顶以及在塔筒两侧发生分离进入尾流区，仅有部分雨滴附着在塔筒表面；

(2) 结构表面雨滴捕捉数量随风速的增大而减少，随雨强的增大而增加，雨滴捕捉数量最多的为工况3(风速为20 m/s、雨强为200 mm/h)，同时表面收集雨滴直径均以5 mm占比最大，且雨滴水平速度主要分布在3～10 m/s范围内；

(3) 不同工况下雨滴撞击位置主要集中在塔筒表面迎风区域两侧各60°范围内，表面雨荷载和雨压系数整体上随着高度的增加先减小后增大，其中雨压系数最大值为0.184，发生在工况3的塔底附近；

(4) 不同工况下相同高度截面处等效压力系数在迎风区与背风区数值略有差异，最大负压值和背风区域负压值均随着高度的增加先增大后减小，其中喉部区域最大负压值约为-1.5，背风面负压平稳区域逐渐变宽；

(5) 本文基于非线性最小二乘法原理提出的冷却塔等效压力系数二维拟合公式可以很好地预测此类极端条件下的塔筒表面荷载取值，最大误差率均控制在10%以内。