基于大涡模拟超大直筒-锥段型钢结构冷却塔极值风压研究

王 浩, 柯世堂

(南京航空航天大学 航空宇航学院， 南京 210016)

作为一种新颖的冷却塔结构形式，大型钢结构冷却塔属于典型的高耸空间网架结构，较混凝土冷却塔自重更轻、柔性更大、阻尼更小，风荷载已成为此类冷却塔结构设计中的控制荷载。同时，直筒-锥段型钢结构冷却塔在用钢量、加工制造和施工安装等方面较双曲线型钢结构冷却塔具有优势，在国际上已有十余例成功应用于电厂间接空冷建设，而在国内此类钢结构冷却塔尚处于预研阶段[1]。

对于直筒-锥段型钢结构冷却塔这类典型风敏感结构，研究人员和工程师们最关注的往往是风荷载的极值，尤其是风压极值作用下冷却塔的局部和整体稳定性能验算[2-4]。而现有冷却塔设计规范[5]和建筑结构荷载规范[6]中仅给出了此类圆柱形断面沿环向的平均风压分布曲线；柯世堂等[7]基于改进的峰值因子法探讨了大型双曲冷却塔表面极值风压的取值问题；李鹏飞等[8]在考虑相关性的基础上给出了双曲线冷却塔外表面极值风压建议取值；罗尧治等[9]基于计算流体动力学(Computational Fluid Dynamics, CFD)方法进行了屋盖结构的风压极值和相关性分析。由于直筒-锥段型钢结构冷却塔特殊的结构形式和气动外形，导致难以寻找相应的设计参考标准，现有规范和文献均没有明确给出此类冷却塔的风压极值。

鉴此，以国内拟建的首座超大直筒-锥段型钢结构冷却塔为例，运用大涡模拟(Large Eddy Simulation, LES)法模拟出该塔外表面三维气动力时程；并通过与国内外大型冷却塔实测及风洞试验结果对比验证了数值模拟结果的有效性。在此基础上，基于三种常用极值估算方法(峰值因子法[10]、改进峰值因子法[11]和Sadek-Simiu法[12])对比分析了此类钢结构冷却塔表面峰值因子和风压极值的分布规律，并基于最小二乘法分别给出了直筒和锥段部分的风压极值一维/二维拟合公式及取值建议。相关研究结论可为此类直筒-锥段型钢结构冷却塔的抗风设计提供参考依据。

1 工程概况与大涡模拟

国内某拟建的直筒-锥段型钢结构冷却塔塔高189 m，由主筒、内部加强桁架、外部附属桁架和外表面的围护层组成。进风口高32.5 m，直径123.2 m，进风口高度以上结构外覆挡风钢板，表面光滑平整。直筒-锥段型钢结构冷却塔主要结构尺寸与构件信息，如表1所示。

表1 直筒-锥段型钢结构冷却塔主要结构尺寸与构件信息表

该冷却塔几何模型按实际尺寸建立，以30%透风率[13]考虑正常运行状态下的百叶窗开启效应。塔筒外表面沿环向逆时针每间隔12°布置一个测点，沿子午向不同高度布置16层测点，共计16×30=480个测点。为保证尾流充分发展，数值风洞尺寸定为24D×15D×4H(顺风向X×横风向Y×高度方向Z，D为底部直径，H为塔高)，冷却塔模型位置距离出口17D，阻塞率为1.21%。将整个计算域划分为外围区域和局部加密区域以兼顾计算效率和精度，外围区域形状规整，采用高质量的结构化网格进行划分，局部加密区域内含冷却塔模型，采用非结构化网格进行划分。加密区最小网格尺寸为0.5 m，总网格数量约640万个，网格质量>0.4。冷却塔测点布置及网格划分，如图1所示。

定义进口边界条件为速度入口，按照B类地貌设置大气边界层指数风速剖面和湍流度剖面，地面粗糙指数为0.15，10 m高度处的基本风速为26.83 m/s，通过用户自定义函数(User Defined Function)实现上述入流边界条件与FLUENT的连接，如图2所示。出口采用压力出口边界条件；顶部和侧面采用对称边界条件；地面及冷却塔表面采用无滑移壁面，如图3所示。基于LES法对冷却塔周围复杂的流场进行模拟。亚格子尺度选用Smagorinsky-Lilly模型，计算残差设置为1×10-6，时间步长定为0.05 s。

2 大涡模拟结果与有效性验证

2.1 数值模拟方法有效性验证

由于现阶段尚无直筒-锥段型冷却塔的风荷载研究见于文献，为验证本文数值模拟方法的有效性，首先进行了相同塔高的双曲线型冷却塔风荷载CFD模拟。图4给出了该塔塔筒喉部外表面体型系数的平均值和根方差沿环向分布曲线，并分别在图中与水工规范[10]和大型双曲冷却塔实测及风洞试验结果[14-15]进行对比。由图4可知，数值模拟得到体型系数平均值与规范曲线吻合良好，且根方差分布曲线与Sageau进行的实测得到的曲线较为接近，考虑到风压信号的脉动程度与实测塔所处的地形、来流湍流和周边干扰密切相关，因此可认为LES法具有一定的有效性。

(a) 测点布置

(b) 整体网格

(c) 局部网格

图2 风速剖面及湍流强度分布

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

2.2 直筒-锥段型冷却塔

前文双曲线冷却塔风荷载大涡模拟结果准确可信的基础上，采用相同的参数设置进行了直筒-锥段型钢结构冷却塔风荷载数值模拟。直筒-锥段型冷却塔典型截面表面绕流及涡量分布图，如图5所示。由图5可知，气流在冷却塔迎风面发生流动分离且局部出现加速效应，气流流动分离后在冷却塔顶部和侧面产生漩涡脱落，持续发展后在冷却塔背风面形成尾流涡旋及回流。

(a) 平均值验证

(b) 根方差验证

分别对直筒和锥段外表面各层体型系数平均值和根方差进行平均，并与水工规范、建筑结构荷载规范和国内外大型双曲冷却塔实测及风洞试验结果进行对比，如图6所示。对比分析可知：① 锥段部分体型系数分布曲线的负压极值点和分离点对应角度与水工规范中双曲线冷却塔风压分布曲线较为一致，直筒段体型系数分布曲线与荷载规范中圆截面构筑物分布曲线基本吻合；② 脉动风压沿环向分布规律与国外双曲冷却塔的实测结果较为接近，但在侧风面和背风面的数值要稍大于实测结果，考虑到直筒-锥段型钢结构冷却塔气动外形与双曲冷却塔差异明显，本文基于大涡模拟得到的脉动风压具有一定的可信度，可用于后续的极值风压研究。

(a) 直筒顶部流线图

(b) 侧立面流线图

(c) 直筒顶部涡量图

(d) 侧立面涡量图

需要说明的是，本文后续研究给出的直筒和锥段峰值因子和极值风压等数据均是对直筒部分和锥段部分所有测点层进行平均的结果。

(a) 平均值验证

(b) 根方差验证

3 极值估计方法

工程设计中脉动风荷载和风振响应极值都应考虑一定的保证系数，计算公式如下

*g*σr

(1)

3.1 峰值因子法

早期Davenport等以高斯分布假定为基础，利用零值穿越理论给出一峰值因子，以该峰值因子乘以风压的根方差后再加上平均风压得到极值风压，此即为峰值因子法。峰值因子按下式计算

(2)

式中:r=0.577 2，为欧拉常数；v表示单位时间内数据穿越平均值的次数；T为计算的时间长度。

这类方法为不少国家的荷载规范所引用，但对于时间T的取值方式略有不同，我国规范的风压时距均为10 min。按照此方法得到的峰值因子偏小，约为3.0～3.5。

3.2 改进峰值因子法

Winterstein等[16]在Davenport的工作基础上，将高斯随机变量表示为考虑高阶统计量(偏度和峰度)的非高斯随机变量的Hermite多项式，从而将Davenport仅仅适用于高斯过程的基于零值穿越理论的峰值因子法扩展到适用于估算非高斯过程的风压极值，此即为峰值因子法。给出非高斯随机变量x(t)和标准高斯过程u(t)的关系式

x=α{u+h3(u2-1)+h4(u3-3u)}

(3)

通过Hermite多项式，可得非高斯过程的峰值因子

(4)

3.3 Sadek-Simiu法

随着统计学的发展以及在工程上的应用，Sadek和Simiu在零值穿越理论的基础上提出了非高斯过程的极值计算方法(Sadek-Simiu法)。时距T内，时程y(t)的极值ypk,T的概率分布函数为

(5)

(6)

式中：v0,y为高斯过程y(t)的零值穿越率，与峰值因子法的含义相同。

因为在映射过程中并不能得到实际的高斯时程y(t)，因此上式在计算时仍然由x(t)的谱Sx(n)代入计算。峰值累积分布函数Fyoi,r(ypi,T)可由式(5)得到，将高斯过程映射到非高斯过程上即可获得非高斯分布时间历程x(t)的峰值累积分布。

4 极值风压研究

4.1 非高斯特性

通过绘制频率分布直方图获得测点风压概率密度分布，并与标准高斯分布曲线进行比较发现环向断面上风压既有高斯分布又包含非高斯分布。限于篇幅，本文仅给出环向典型非高斯测点的风压时程曲线和概率密度分布曲线，如图7和图8所示。由图可知，冷却塔侧风面位于来流风分离区，风压时程曲线具有明显的不对称性，伴随有大量间歇性脉冲信号并伴随有很大的风吸力。从更多测点体型系数的概率分布看，结构表面有大部分区域的风压信号概率分布属于非高斯分布。

(a) 负压极值区

(b) 分离区

(a) 负压极值区

(b) 分离区

4.2 峰值因子取值探讨

由三种极值风压算法计算得出的直筒段和锥段沿环向峰值因子平均值，如图9所示。对比可以发现三种不同方法计算得出的结果差异较大，基于高斯分布假定的峰值因子法计算得出的数值普遍偏低，且沿环向角度变化微弱；锥段部分背风区的峰值因子较直筒段背风区更大，说明此类钢结构冷却塔在锥段区域的背风面可能会产生更为明显的风压极值。

(a) 锥段部分

(b) 直筒部分

不同算法计算得到的峰值因子二维分布图，如图10所示。对比发现峰值因子分布呈明显的二维分布特征，改进峰值因子法和Sadek-Simiu法计算得到的峰值因子沿高度变化明显，尤其是在锥段和直筒顶部等三维绕流明显的区域峰值因子变化显著。三种算法计算得到的峰值因子最大值均出现在背风区，这是由于背风区风压时程的均值较小，且对应脉动值相对较大。

(a) 峰值因子法

(b) 改进峰值因子法

(c) Sadek-Simiu法

4.3 风压极值对比研究

三种方法计算得出的直筒段和锥段体型系数极值曲线，并与规范曲线进行对比，如图11所示。对比发现锥段迎风面和侧风面的极值风压接近于双曲冷却塔推荐曲线，而直筒段迎风面和侧风面的极值风压与大型圆截面构筑物风压极值曲线更接近。三种方法计算得到的锥段背风区的极值风压分别比规范数值大35.3%、58.9%和39.5%，直筒段计算结果分别比规范数值大37.4%、54.9%和45.6%。

(a) 锥段部分

(b) 直筒部分

Fig.11 The comparison diagram of extreme shape coefficient under different methods with standard value

直筒-锥段型冷却塔体型系数极值二维分布图，如图12所示。由图12可知，① 极值风压分布关于风轴方向比较对称，最大正值和负值分别位于冷却塔的迎风面和侧风面，最大负值集中分布于60～125 m高度区间内；② 锥段侧风区风压负值大小及分布面积明显小于上部直筒段，背风区域极值风压分布也与直筒段有很大区别。对比表明改进峰值因子法计算得出的风压极值基本包络住其它两种方法的结果和规范值，后续极值拟合研究均基于改进峰值因子法的计算结果。

4.4 极值风压拟合

上述研究表明直筒-锥段型钢结构冷却塔极值风压沿环向和子午向分布呈现出典型的二维特征，而相关设计规范只给出了双曲线冷却塔和圆截面构筑物的的一维风压极值曲线，该方法对直筒-锥段型钢结构冷却塔而言不尽合理。为方便工程研究与设计人员精确获得此类冷却塔极值风压，本文基于非线性最小二乘法原理，以子午向高度比n和环向角度θ为目标函数，拟合给出此类冷却塔二维风压极值的计算公式，同时给出了以环向角度θ为目标函数的直筒段和锥段分段一维平均拟合公式，供设计参考。

锥段和和直筒段原始数据与一维拟合曲线对比示意图，如图13所示。由图13可知，公式拟合结果与实际数据一致，拟合公式定义如式(7)所示。具体参数见表2和表3。

(7)

(a) 峰值因子法

(b) 改进峰值因子法

(c) Sadek-Simiu法

(a) 锥段

(b) 直筒部分

参数a1a2a3b1b2b3c1c2c3数值2．40300．0093-3．24105．96900．04200．28875．61900．04393．0860

表3 直筒部分目标拟合公式参数表

直筒-锥段型钢结构冷却塔二维极值风压拟合公式定义为

Cpm(n,θ)=a00+a10θ+a01n+a20θ2+a11θn+a02n2+

a21θ2n+a12θn2+a03n3+a22θ2n2+a13θ1n3+

a04n4+a23θ2n3+a14θn4+a05n5

(8)

式中：n=z/H，其中z为测点对应高度；H为塔顶总高度；aij(i=1～2,j=1～5)为拟合系数；具体数值，见表4。

直筒-锥段型钢结构冷却塔极值风压实际分布与二维拟合曲面对比图，如图14所示。图14中散点值即为冷却塔真实极值风压，曲面对应数值为根据二维拟合公式计算得到的极值风压。从拟合曲面的整体分布来看，其沿子午向和环向的变化规律与实际极值风压基本一致，数值上可以包络住真实数值，误差率<10%。对比结果表明本文提出的二维极值风压拟合公式可作为此类大型直筒-锥段型钢结构冷却塔极值风压取值的参考依据。

表4直筒-锥段型冷却塔二维极值风压目标拟合公式系数表

Tab.4Twodimensionalextremewindpressurefittingformulacoefficienttableofcylindrical-conecoolingtower

拟合系数a001．245a0532．4400a140．003036a0113．61a10-0．1440a200．001526a02-64．49a110．0332a21-0．000109a03123．30a12-0．0690a220．000219a04-104．30a130．0408a23-0．000129

图14 极值风压原始数据与二维拟合曲面对比图

Fig.14 Comparison between the raw data and the two dimensional fitted surface of extreme wind pressure

5 结 论

本文系统详尽地研究了直筒-锥段型钢结构冷却塔表面风压极值分布特性，主要涉及大涡模拟、峰值因子计算方法、风压极值计算方法和非线性最小二乘法拟合公式。得到如下主要结论：

(1) 通过与国内外实测及风洞试验结果对比验证了LES方法和结果的有效性，风荷载时程可用于后续风压极值研究。

(2) 直筒-锥段型钢结构冷却塔负压极值区和分离区附近风荷载呈现明显的非高斯特性，基于高斯分布假定的峰值因子法计算得出的极值风压数值偏小，改进峰值因子法计算结果基本包络住其他两种方法，Sadek-Simiu法计算得到的风压极值介于二者之间。

(3) 锥段和直筒段背风面极值风压均大于规范数值，锥段迎风面和侧风面极值风压与水工规范极值接近，而直筒段迎风面和侧风面极值风压与荷载规范极值接近。

(4) 本文提出的二维风压极值拟合公式可包络直筒和锥段的实际风压极值，最大误差控制在<10%，可对此类钢结构冷却塔风荷载设计取值预测提供参考。

[1] 邹云峰, 李寿英, 牛华伟,等. 双曲冷却塔等效静力风荷载规范适应性研究[J]. 振动与冲击, 2013, 32(11):100-105.

ZOU Yunfeng， LI Shouying, NIU Huawei, et al. The hyperbolic cooling tower equivalent static wind load code for the adaptive study[J]. Journal of Vibration and Shock, 2013, 32 (11): 100-105.

[2] 柯世堂, 葛耀君, 赵林. 大型双曲冷却塔表面脉动风压随机特性—非高斯特性研究[J]. 实验流体力学, 2010, 24(3): 12-18.

KE Shitang, GE Yaojun, ZHAO Lin. Research on festures of fluctuating wind pressure on large hyperbolic cooling tower: features of non-Gaussian[J]. Journal of Experiments in Fluid Mechanics, 2010, 24(3): 12-18.

[3] 全涌, 顾明, 陈斌,等. 非高斯风压的极值计算方法[J]. 力学学报, 2010, 42(3):560-566.

QUAN Yong, GU Ming, CHEN Bin, et al. Study on the extreme value estimating method of non-Gaussian wind pressure[J]. Chinese Journal of Theoretical and Applied Mechanics, 2010, 42(3): 560-566.

[4] GIOFFRE M, GUSELLA V. Damage accum ulation in glass plates[J]. Journal of Engineer Mechanics, 2002(7): 801-805.

[5] 工业循环水冷却设计规范:GB/T 50102—2014[S]. 北京: 中国计划出版社, 2014.

[6] 建筑结构荷载规范:GB 50009—2012[S]. 北京: 中国建筑工业出版社, 2012.

[7] 柯世堂, 葛耀君, 赵林. 大型双曲冷却塔表面脉动风压随机特性——风压极值特性探讨[J]. 实验流体力学, 2010, 24(4): 7-12.

KE Shitang, GE Yaojun, ZHAO Lin. Research on festures of fluctuating wind pressure on large hyperbolic cooling tower: discussions on extreme wind pressure[J]. Journal of Experiments in Fluid Mechanics, 2010, 24(4): 7-12.

[8] 李鹏飞, 赵林, 葛耀君, 等. 超大型冷却塔风荷载特性风洞试验研究[J]. 工程力学, 2008, 25(6):60-67.

LI Pengfei, ZHAO Lin, GE Yaojun, et al. Wind tunnel investigation on wind load characteristics for super large cooling towers[J]. Engineering Mechanics, 2008, 25(6):60-67.

[9] 罗尧治, 丁慧. 波纹表面月牙形悬挑屋盖风荷载特性[J]. 浙江大学学报(工学版), 2013, 47(6):1063-1071.

LUO Yaozhi, DING Hui. Wind load characteristics of crescent-shaped cantilevered roof with wavy surface[J]. Journal of Zhejiang University(Engineering Science), 2013, 47(6):1063-1071.

[10] DAVENPORT A G. Gust loading factors[J]. Journal of Structural Division, 1967, 93(ST3): 11-34.

[11] KAREEM A, ZHAO J. Analysis of non-Gaussian surge response of tension leg platforms under wind loads[J]. Journal of Offshore Mechanics and Arctic Engineering, 1994, 146: 137-144.

[12] SADEK F, SIMIU E. Peak non-Gaussian wind effects for database-assisted low-rise building design[J]. Journal of Engineering Mechanics, 2002,128: 530-539.

[13] KE S T, JUN L, ZHAO L, et al. Influence of ventilation rate on the aerodynamic interference between two extra-large indirect dry cooling towers by CFD[J]. Wind & Structures an International Journal, 2015, 20(3):449-468.

[14] NISHIMURA H, TANIIKE Y. Aerodynamic characteristics of fluctuating forces on a circular cylinder[J]. Journal of Wind Engineering & Industrial Aerodynamics, 2001, 89(1): 713-723.

[15] CHENG X X, ZHAO L, GE Y J, et al. Wind pressures on a large cooling tower[J]. Advances in Structural Engineering, 2015, 18:201-220.

[16] WINTERSTEIN S R. Nonlinear vibration models for extremes and fatigue[J]. Journal of Engineering Mechanics, 1988, 114(10): 1772-1790.