下击暴流作用下超大型冷却塔风场驱动机理与风荷载极值模型

韩光全，柯世堂，杨 杰，李文杰，任贺贺

(南京航空航天大学 民航学院 土木与机场工程系，南京 211106)

下击暴流[1]又称雷暴冲击风，其主要特征是强下沉气流和近地面爆发状强风。雷暴天气中，下落雨滴伴随着湿冷空气产生局部强下行气流，撞击地面后沿四周迅速扩散，在近地表面形成放射状的强风，该流动模式通常被认为反向的龙卷风。不同于大气边界层良态风场，下击暴流最大风速产生在距离地面约50～100 m 高度处[2]，且气象记录表明，其最大风速可达75 m/s[3]。同时，下击暴流在核心射流区域和突出前缘均存在明显的垂直分量，在与强切变的水平风速共同作用下，造成了建构筑物的严重损伤甚至破坏[4-5]。超大型冷却塔作为火/核电厂重要的基础设施，是典型的高耸旋转薄壳建筑，风荷载是其结构设计的控制荷载之一。然而，GB/T 50102—2014《工业循环水冷却设计规范》[6]仅给出良态风作用下的荷载，对下击暴流等特异风作用下冷却塔风荷载分布模式未给予考虑。因此，对下击暴流作用下冷却塔风场作用机理和风荷载特性进行研究具有现实意义，可为超大型冷却塔设计和建设提供重要的参考依据。

对于下击暴流三维风场，限于其短时突发性[7]与移动性[8]，尚且难以全面即时地进行现场实测，目前国内外学者主要采用物理试验和数值模拟等手段对其进行研究。物理试验一般采用主动控制风洞[9-10]或冲击射流装置[11-13]，对结构表面风压和空气动力学特性进行研究，但是由于难以满足缩尺比等试验条件，物理试验难以对复杂结构形式的建构筑物进行下击暴流作用下的准确模拟。随着计算机技术的发展，数值模拟逐渐被应用到下击暴流三维风场研究中，利用数值模拟技术，可对风场及其对结构的作用机理[14-16]进行较为全面的探究。文献[17-18]采用冲击射流模拟装置，对下击暴流作用下高层和双坡屋面建筑风场流动特性进行了分析，结果发现在风场一定位置处建筑表面风压显著高于ASCE 7-05标准[19]。文献[20-21]分别借助可移动的冲击射流装置和滑移网格技术对移动型下击暴流进行试验和数值模拟，研究表明，移动的下击暴流作用于建筑物上会使建筑表面出现整体变化，甚至超出规范良态风风压限值。文献[22-23]研究发现下击暴流会引起输电塔间纵向拉力急剧增大或者关键结构构件的局部失效，进而引起失稳倒塌。对下击暴流作用下大跨度桥梁[24-25]、风力机[26-27]、屋盖厂房[28-30]等结构的研究也显示相比于良态风，下击暴流会不同程度增加结构物的潜在破坏性。而目前鲜有涉及下击暴流三维脉动风场对超大型冷却塔的风场作用机理与风荷载特性，亟需开展特异风作用下超大型冷却塔风荷载取值研究。

鉴于此，本文基于冲击射流模型与大涡模拟(large eddy simulation,LES)技术，首先对下击暴流三维风场进行了全生命周期非定常数值模拟，将模拟结果与经验模型和实测数据进行对比验证，并阐释其风场特性；在此基础上，详细分析了超大型冷却塔在下击暴流不同径向位置处的绕流特性，升/阻力系数分布特征，以及瞬态风压和极值风荷载分布规律，结论可为此类超大型冷却塔抗雷暴设计提供一定的参考。

1 工程概况

该在建世界最高超大型冷却塔位于我国内蒙地区，设计塔高突破规范190 m的限定，冷却塔主体结构由塔顶刚性环、塔筒、支柱、环基及加劲肋构成。其中，塔筒采用双曲线型钢筋混凝土壳体结构，其厚度呈指数变化，最小与最大壁厚分别为0.42 m和2.45 m，塔筒外表面沿环向均匀布置120条梯形加劲肋，60榀X型支柱通过60个支墩与底部钢筋混凝土环板型基础承台连接，单个支墩下布24根长30 m直径0.8 m的钢筋混凝土灌注桩，冷却塔主要结构尺寸参数如表1所示。

表1 冷却塔主要结构尺寸Tab.1 The main structural sizes of cooling tower

2 数值模拟计算概况

2.1 数值模型与工况设置

为全面监测冷却塔各处风压，在冷却塔模型内、外表面沿子午向布置14层测点，每层沿环向逆时针间隔12°布置一个风压系数测点，内、外表面测点一一对应，则环向共布置60个测点，塔筒内、外表面共布置840个测点。

CFD数值模拟采用三维足尺模型，冲击射流速度入口直径Djet=600 m，射流高度为2Djet，为保证流动的充分发展，计算域尺寸为10Djet×10Djet×3Djet。边界条件设置如下：射流入口采用速度入口，射流初始风速Vjet=29 m/s，湍流强度为1%；计算域四周及顶面采用压力出口边界条件，回流湍流强度同为1%；地面采用无滑移壁面，射流口上部圆柱面设为滑移壁面。计算域及边界条件如图1(a)所示。

考虑到下击暴流风场与冷却塔结构特性，为探究冷却塔在风场不同位置处的风荷载分布特性，共设置4个典型计算工况。其中：工况1为冷却塔置于射流入口正下方；工况2～工况4分别为塔筒底部迎风面前端距射流中心径向距离r=1.0Djet，r=1.5Djet，r=2.0Djet。工况设置如图1(b)所示。

图1 计算域边界条件及工况设置示意图Fig.1 Calculation domain boundary conditions and working condition settings

2.2 网格划分与参数设置

假定冲击射流沿高度方向均为尺寸统一的圆形，而O型拓扑适合此类两端均为圆形截面的物体，故在射流中心区域采用O-Block二次拓扑生成双层O型网格。同时为满足冷却塔结构周围的模拟精度，采用混合网格离散形式，将整个计算域分为内、外两个部分：核心区域采用四面体网格，并对冷却塔加劲肋周围网格进行局部加密，外围区域采用规则拓扑的高质量六面体结构网格，细部及整体网格划分如图2所示。

图2 细部及整体网格划分示意图Fig.2 Diagram of detailed and overall grid subdivision

不同网格数量下r=1.0Djet径向位置处的风剖面，如图3所示。由图3可知，当网格数量达到1 200万时，在上部区域风剖面随网格数量的增加无明显变化，即此时数值模型已基本满足与网格数量无关的要求。考虑到计算效率和准确性，后续研究中使用的网格数为1 600万。

图3 不同网格数量下风剖面对比示意图Fig.3 Comparison of wind field results with different gird numbers

此网格方案中近壁面网格高度为0.001 m，壁面法向网格的增长率为1.08，宽度为15 m，故近壁面最大网格高宽比为6.67×10-5。模拟采用的时间步长为Δt=0.03 s，冷却塔表面最小网格尺寸为Δx=1.2 m，由

Courant=Δt·Δx/Vjet

(1)

计算可得Courant number为0.725，其值小于1，可判定时间步长与网格是匹配的。

基于Fluent程序进行计算，在LES模拟中，亚格子模型采用algebraic wall-modeled LES(WMLES)模型以弥补近壁面网格数量不足的问题，从而实现雷暴对流系统的高精度模拟需求。为生成随时间变化的入口边界条件，采用Vortex method，Number of Vortices设置为190；雷诺应力定义方法采用k or turbulent intensity。速度-压力耦合方法采用压力耦合方程组的半隐式方法(semi-implicit method for pressure linked equations，SIMPLE)，梯度插值方案采用least squares cell based，压力采用standard格式，对动量的计算采用second order upwind格式以实现高阶离散化。质量残差设置为10-6，x，y，z方向速度收敛标准设置为10-3。在计算时，采用Data Sampling for Time Statistics，sampling interval设置为1。

3 下击暴流风场模拟

3.1 有效性验证

首先对没有冷却塔结构的下击暴流风场进行模拟，采用无量纲时间T=Δt·Vjet/Djet，为确保风场模拟的精确和完整性，无量纲时间步长为0.001 45，共模拟20 000步，总模拟时长T=29。

r=1.0Djet径向位置处径向风速竖向风剖面和z=0.2Djet高度处径向风速径向风剖面，分别如图4和图5所示。其中：图4横坐标为径向风速u与最大径向风速umax的比值，径向风速u为下击暴流发展成熟阶段(T=8.7～14.5)瞬时风速的平均值，纵坐标为离地高度z与最大径向风速对应高度zmax的比值；图5横坐标为径向距离r与最大径向风速对应径向距离rmax的比值。从图5中可知，数值模拟得出的归一化风剖面与各经验模型[31-34]及实测数据吻合度较高，说明本文采用冲击射流模型和LES技术可以实现对下击暴流风场的有效模拟。

图4 径向风速归一化竖向风剖面示意图Fig.4 The vertical profile of the radial velocity

图5 径向风速归一化径向风剖面示意图Fig.5 The radial profile of the radial velocity

3.2 风场特性

下击暴流部分时刻沿流域纵剖面的风速云图，如图6所示。T=5.8时刻为下击暴流发生初期阶段，此时冲击风自射流入口向下方运动，气流前端尚未到达地面，极值风速为33.51 m/s；随着时间的推移，在T=8.7时刻，冲击风前缘已经抵达地面，与地面接触并开始径向对流，触发边界层气流的非定常分离，形成反向旋转、径向移动的第一个涡环，极值风速增大到41.00 m/s；在射流入口逐渐有稳定气流进入的情况下，于第一个涡环之后产生一系列尺度较小的涡环，极值风速也增大到约47.30 m/s左右，如图6(c)、图6(d)所示。

图6 下击暴流发展阶段不同时刻的纵剖面风速云图Fig.6 Wind speed nephogram of longitudinal section at different times of downburst development stage

不同径向位置处近地面风速时程曲线，如图7所示。在r=0径向位置处，竖向风速为主导分量，随着冲击风的持续输入，近地面各测点风速从0逐渐增大，在约T=7.25时刻达到最大值，10～60 m高度处各测点极值风速分别约为1 m/s，2 m/s，3 m/s，4 m/s，5 m/s，6 m/s，而后略有减小并维持稳定波动状态。

在r=1.0Djet，r=1.5Djet，r=2.0Djet径向位置处，水平风速为主导分量，由于气流撞击地面并向四周逐渐扩散，各径向位置处近地面测点风速首先达到首个峰值，且达到首个峰值的时间随径向距离的增加而延迟，首个峰值的最大值分别约为38 m/s，50 m/s，22 m/s。随着径向移动、反向旋转涡环的不断生成，风速相应呈现波动趋势，形成一系列“波峰”与“波谷”。其中r=1.0Djet，r=1.5Djet两处一系列风速“波峰”与“波谷”的数值较为稳定，稳定“波峰”的峰值均维持在40 m/s左右，r=1.5Djet处“波谷”风速小于r=1.0Djet处，可见r=1.5Djet处涡环运动相较于r=1.0Djet处更为剧烈。r=2.0Djet径向位置由于远离入流中心，湍流发展紊乱且能量削弱，一系列风速“波峰”与“波谷”呈现逐渐衰减的趋势。

图7 近地面风速时程曲线Fig.7 Time history curves of near-surface wind speed

4 结构流场与风荷载特性

4.1 绕流特性

由图7可知，各工况在T=15.0之后风速基本维持稳定变化，即下击暴流在此时刻后为发展成熟阶段，故图8给出了在T=17.4时刻沿y=0纵剖面冷却塔周围风压云图与空气流线图。工况1中，下冲气流在塔筒顶部产生分离，进入塔筒内部的气流从X型支柱处流出，与外部下行气流产生对流，从而在塔筒外部形成大尺度涡旋，涡旋处风压相对周围区域较小，而塔筒内部风压较高且分布均匀；工况2～工况4中，流动的复杂性相比工况1愈加明显，近地面气流一部分沿塔筒外表面爬升，与上部下行气流在迎风面塔顶冲击，形成涡旋；另一部分近地面气流从X型支柱处进入塔筒内部，撞击塔筒内壁后形成涡旋，此3个工况中塔筒迎风区外表面和背风区内表面为高风压区。

图8 沿y=0剖面冷却塔周围风压云图与流线图Fig.8 Pressure cloud and streamline diagram around the cooling tower along the y=0 section

为研究在水平方向上冷却塔结构内、外空气流线变化，图9给出了在T=17.4时刻沿z=40 m横剖面冷却塔周围风压云图与空气流线图。由图9可知：工况1中由于冲击风竖直流下，故塔筒内部几乎没有水平方向气流，塔筒外部气流呈现辐射状发展，且存在大小不一的漩涡；工况2～工况4中，高压区同样主要分布在塔筒迎风面外部和背风面内部，随着冷却塔迎风面与入流口的距离越来越远，高压区的范围和大小随之逐渐变小，并且塔筒前端气流由集束状态逐渐转向水平状态，塔筒内部和背风面尾流区域均存在漩涡。

图9 沿z=40 m剖面冷却塔周围风压云图与流线图Fig.9 Pressure cloud and streamline diagram around the cooling tower along the z=40 m section

4.2 风压系数瞬态分布

定义风压系数为

(2)

式中：P为测点的风压；Pref为参考静压力；ρ为空气密度；Vref为参考风速，在常规大气边界层风场中，通常取为模型顶部平均风速，而在下击暴流风场模拟中，由于径向风速最大值一般出现在近地面附近，且大小随径向位置的变化而变化，因而参考风速一般取为射流入口的初始风速Vjet。

工况1中冷却塔位于冲击风正下方，塔筒内、外表面风压系数沿环向分布规律基本一致。因此，图10给出了顶部、喉部以及底部3个典型断面测点风压系数时程曲线。由图10可知：塔筒内、外表面3个不同断面时程风压系数变化趋势基本相同，在初始时间范围内，风压系数先减小再增大，在约T=7.25时刻达到最大值，而后趋于近似正弦曲线的波动状态，外表面底部、喉部、顶部测点风压系数分别在0.8，0.7，0.9附近波动，波动趋势稍显杂乱，而内表面相应测点风压系数分别在0.9，0.7，0.8附近波动，波动趋势更为平稳；由于气流到达塔筒各位置的时间不同，故不同断面测点风压系数在正弦曲线变化阶段存在相位差，此现象在外表面更为明显。不同于良态风作用下冷却塔风压分布特征，此工况下喉部断面风压系数整体最小。

工况2～工况4中塔筒底部、喉部、顶部典型断面在迎风区、侧风区、背风区测点的风压系数时程曲线，分别如图11～图13所示。由图11～图13可知，各工况下相应测点时程风压系数变化规律相似。以工况2为例分析如下：

底部外表面迎风区测点的风压系数由0.8减小至0.3左右，随后增加达到首个峰值，而侧风区和背风区测点则是由0.8不断减小至负压峰值，各测点达到首个峰值的时刻与r=1.0Djet处风速达到首个峰值的时刻相吻合，此后风压系数分别在1.0，-0.2，-0.1左右维持稳定小幅震荡。

图10 工况1典型断面测点风压系数时程曲线Fig.10 Time history curves of wind pressure coefficient at typical cross-section with working condition 1

图11 工况2典型断面测点风压系数时程曲线Fig.11 Time history curves of wind pressure coefficient at typical cross-section with working condition 2

图12 工况3典型断面测点风压系数时程曲线Fig.12 Time history curves of wind pressure coefficient at typical cross-section with working condition 3

图13 工况4典型断面测点风压系数时程曲线Fig.13 Time history curves of wind pressure coefficient at typical cross-section with working condition 4

喉部测点风压系数由0.7左右逐渐减小，而后在0左右维持小幅震荡。但由于涡环主要发生在近地面，故峰值现象明显弱于底部，同理可为顶部测点时程风压系数变化作出解释。由于塔内受气流湍流扰动程度较小，故在首个峰值时刻后，各断面内表面测点风压系数震荡较外表面更为平缓。

另一方面，对比不同工况测点时程风压系数差异可知：达到首个峰值的时刻随着径向距离的增加而延迟，分别对应该径向位置处风速达到首个峰值的时刻；由于r=1.5Djet径向位置处湍流发展紊乱，故工况3中风压系数时程曲线震荡幅度较大，且底部外表面测点峰值绝对值较大。

4.3 升/阻力系数

冷却塔升/阻力系数由风压系数积分计算得到[35]

(3)

(4)

式中：Ai为第i测点压力等效覆盖面积；θi为第i测点压力与来流方向夹角；A为整体结构沿顺风方向投影面积。由于工况1的特殊性，在此不对其进行升/阻力系数计算。

工况2塔筒底部、喉部、顶部3个典型测点层的升/阻力系数时程曲线，如图14所示。由图14可知，各测点层的升/阻力系数在下击暴流初始时间范围内均维持在0左右，随后产生波动，波动的幅度有所差别。整体上，外表面各层的升/阻力系数波动幅度大于内表面，底部测点层的阻力系数波动幅度最大。其中：外表面均为正值，最大可达到1.5；内表面均为负值，最大可达到-0.5。工况3、工况4时程升/阻力系数变化规律与工况2相似，在此不再附图。

对内、外表面各测点层的时程升/阻力系数进行平均，如图15所示。由图15可知，除工况2中冷却塔喉部以下外表面测点层外，内外表面各层的平均升力系数基本为0，外表面阻力系数自塔顶沿塔高方向逐渐增大，在塔底达到正值最大值，内表面阻力系数变化趋势相似，在塔底达到负值最大值。随着冷却塔距离射流中心越来越远，塔底阻力系数数值由大变小，即工况2中阻力系数数值最大，外表面阻力系数为0.5，内表面阻力系数为-0.2。

图14 典型测点层升/阻力系数时程曲线Fig.14 Time history curves of lift/drag coefficient of typical measuring point layer

4.4 极值风荷载

下击暴流对超大型冷却塔具有强烈的冲击性，将冷却塔各典型断面测点内、外表面时程风压系数进行叠加，并取其峰值，如图16所示。同时给出了考虑内、外压的规范良态风作用下风压分布曲线，并以此计算得到考虑极值风效应的风压系数包络值：-0.92～2.89。由图16分析可知：

(1) 工况1中各断面测点总风压系数峰值基本在-0.5上下小幅波动，处在包络范围内，究其原因为冷却塔壁厚较薄，此工况中气流沿筒壁近似竖直流下，故内、外表面瞬时风压系数相差较小。此工况下按规范进行结构设计偏于安全。

(2) 工况2～工况4中各断面测点总风压系数峰值呈现明显的对称分布特性，沿环向分布规律与规范良态风曲线相似，均未超出包络正值2.89，但均存在超出包络负值-0.92的区域。其中：工况2、工况3底部背风区总风压系数峰值分别为-1.46，-1.52，超出包络范围分别为58.70%，65.22%；另外工况3、工况4中塔筒中部和喉部断面在侧风区同样存在小幅度超出包络范围的区域。

下击暴流短时突发形成的涡环冲击作用在冷却塔内、外表面，极有可能引起瞬时极值风荷载超出规范限值，进而易引起结构的破坏，需引起设计人员的注意。

图15 各测点层平均升/阻力系数分布曲线Fig.15 Distribution curve of average lift/drag coefficient of measuring point layers

图16 塔筒典型断面测点总风压系数峰值曲线图Fig.16 Curves of peak values of total wind pressure coefficient of typical tower sections

5 结 论

(1) 采用冲击射流模型和大涡模拟技术模拟得出的径向风速归一化风剖面与经验模型、实测数据吻合度较高；下击暴流冲击气流与地面接触碰撞后，产生一系列径向移动，反向旋转的涡环，在r=1.0Djet，r=1.5Djet，r=2.0Djet径向位置处风速随之呈现波动变化趋势。在射流初始风速为29 m/s的情况下，r=1.5Djet处极值风速最高可达50 m/s。

(2) 在下击暴流发展成熟阶段，工况1中冷却塔塔筒外表面附近区域存在大尺度涡旋，生成涡旋部位风压小于周围区域；而在工况2～工况4中涡旋主要存在于塔筒内部、迎风面塔顶区域和背风面近地尾流区，迎风区外表面和背风区内表面风压高于其他区域。

(3) 塔筒内、外表面风压系数时程分布三维效应明显，工况1中典型区域测点风压系数随时间先减小后增大，随后在下击暴流发展成熟阶段维持稳定的正弦分布状态；而工况2～工况4风压系数脉动趋势更为显著，由大幅震荡发展成为维持在固定数值附近的小幅震荡，其中底部区域测点受涡环影响震荡尤为明显。

(4) 层平均升力系数较小，工况2中塔底层平均阻力系数最大，在内、外表面分别可达到-0.2，0.5；工况1处在规范极值风压包络范围内，工况2～工况4均有区域超出包络范围，其中工况2、工况3的底部区域超出幅度分别达到58.70%，65.22%。

综上所述，下击暴流对超大型冷却塔的风场驱动主要表现为近地面涡环的冲击作用，其瞬态极值风压发生位置从规范良态风作用下的喉部侧风区转移到底部背风区，同时最大数值较规范良态风作用增大65.22%。主要结论可为此类特异风作用下超大型冷却塔结构抗风设计提供科学依据。