钢结构冷却塔作为一种新颖的冷却塔类型，具有可回收、施工快且不受冬季温度限制等优点，与当前冷却塔高大化、多样化的发展趋势相适应。目前钢结构冷却塔在国外已有十余例成功应用于电厂的间接空冷建设^[1]，而在国内尚处于萌芽阶段。与传统的钢筋混凝土冷却塔^[2]结构形式相比，钢结构冷却塔阻尼更小、频率更低且振型复杂，其风致动力放大效应和受力性能问题更为突出。因此，表面风荷载数值模拟和分布特性研究成为钢结构冷却塔抗风设计亟待解决的首要问题。

自1965年英国渡桥电厂冷却塔风毁事故^[3]以来，国内外学者对大型冷却塔的风荷载特性进行了系统深入的研究，文献[4]早在1986年对大型冷却塔进行了塔群风洞试验，用以研究塔群的平面布局和来流湍流特性对塔群风荷载的影响；文献[5]于1992年对两座高120 m的双曲线型冷却塔群进行了风压分布的全尺寸测量；文献[6]于2007年针对单体和双塔情况下冷却塔的风荷载进行了CFD数值模拟，验证了CFD用于模拟冷却塔风荷载的可靠性；文献[7]于2012年提出了一种新的频域方法对超大型双曲冷却塔等效静风荷载进行了研究；文献[8]于2015年通过刚性模型测压风洞试验分析了某特大型冷却塔的表面风荷载三维效应。已有成果为大型双曲冷却塔的风荷载设计取值提供了很好的指导作用，但考虑到直筒-锥段型钢结构冷却塔气动外形与钢筋混凝土双曲冷却塔相差甚远，若直接采用现行冷却塔规范的平均风压系数分布八项式进行设计值得商榷，且国内外更缺乏此类结构形式新颖的直筒-锥段型钢结构冷却塔表面脉动风压分布随机特性研究。

鉴于此，以国内拟建的某新型超大直筒-锥段型钢结构冷却塔为例，采用LES方法对其进行了数值风洞模拟，基于计算得到的冷却塔外表面三维风压时程，对比国内外大型冷却塔实测及风洞试验结果验证了数值模拟方法的有效性，提炼出了此类直筒-锥段型钢结构冷却塔外表面平均、脉动和极值风荷载以及阻力系数等气动参数的分布特性，最后给出了平均和极值风压分布曲线和八项拟合式。研究结论可为我国此类超大型钢结构冷却塔的风荷载设计取值提供参考依据。

1 工程简介

该拟建超大直筒-锥段型冷却塔，整体塔高189 m，进风口高度32.5 m，直径144.5 m。塔体采用钢框架支承，由主筒、内部加强桁架和外部附属桁架组成，主筒共18层，加强桁架共5层，分别设置在32.5，67，107，148和189 m高度处，主筒和加强桁架均采用Q345钢材；附属桁架共30榀，钢材等级为Q235B。进风口高度以上结构外覆挡风钢板，表面光滑平整；塔筒底部设置百叶窗，用以调节进风量。冷却塔主要结构尺寸和示意图见表 1。

2 大涡模拟 2.1 流体控制方程

结构抗风研究中将流体视为黏性不可压缩，对瞬态的N-S方程进行空间平均，可得大涡模拟方法的控制方程为

$ \frac{{\partial {{\bar \mu }_i}}}{{\partial {x_i}}} = 0 $

(1)

$ \frac{{\partial {{\bar \mu }_i}}}{{\partial t}} + \frac{{\partial \left( {\bar \mu {\mu _i}{{\bar \mu }_j}} \right)}}{{\partial {x_j}}} = - \frac{1}{\rho }\frac{{\partial \bar p}}{{\partial {x_i}}} + v\frac{{{\partial _2}{{\bar \mu }_i}}}{{\partial {x_j}\partial {x_j}}} - \frac{{\partial {\tau _{ij}}}}{{\partial {x_j}}} $

(2)

式中：ρ为空气密度；t为时间；v为空气运动黏性系数；${\bar \mu _i}$，${\bar \mu _j}$表示滤波后空间3个方向的速度，i和j均取为1，2和3；x_i和x_j为坐标分量；p为滤波后压力；τ_ij为空间平均后N-S方程的不封闭项即亚格子应力，其表达式如下

$ {\tau _{ij}} = {{\bar \mu }_i}{{\bar \mu }_i} - {{\bar \mu }_i}{{\bar \mu }_j} $

(3)

根据Smagorinsky提出的基于涡旋黏度假设的亚格子模型，引入Boussinesq假设，亚格子应力可转化为

$ {\tau _{ij}} - \frac{1}{3}{\tau _{ij}}{\delta _{ij}} = - 2{\mu _t}{{\mathit{\boldsymbol{\bar S}}}_{ij}} = - {\mu _t}\left( {\frac{{\partial {{\bar \mu }_i}}}{{\partial {x_j}}} + \frac{{\partial {{\bar \mu }_j}}}{{\partial {x_i}}}} \right) $

(4)

式中：S_ij为可解尺度应变率张量；δ_ij为Kronecker delta函数；μ_t为亚格子湍流涡黏系数，通常采用Smagorinsky假设

$ {\mu _t} = {\left( {{C_{\rm{s}}}\Delta } \right)^2}\left| {\mathit{\boldsymbol{\bar S}}} \right| $

(5)

式中：C_s为Smagorinsky常数，位于0.1~0.23间，文中取为0.1；Δ为空间网格尺度；$\left| {\mathit{\boldsymbol{\bar S}}} \right|$应变张量。此为标准的Smagorinsky亚格子模型。一些学者提出动态确定C_s值的方法，可以更好地考虑钝体周围流场的冲撞、分离、自由剪切层和涡脱落等现象，即Dynamic Smagorinsky。

然而Smagorinsky与Dynamic Smagorinsky均属于代数模型，其假定亚格子湍动能的产生和耗散局部平衡，不能很好地模拟冷却塔此类高雷诺数结构。为解决该问题，本文基于文献[9, 10]采用一种新的亚格子模型，基本方程如下

$ \begin{array}{*{20}{c}} {\frac{{\partial {k_{{\rm{sgs}}}}}}{{\partial t}} + \frac{{\partial {\mu _j}{k_{{\rm{sgs}}}}}}{{\partial {x_j}}} = - {\tau _{ij}}{S_{ij}} - C\varepsilon \frac{{{k_{{\rm{sgs}}}}^{3/2}}}{\Delta } + }\\ {\frac{\partial }{{\partial {x_j}}}\left[ {\left( {{C_d}\Delta v\sqrt {{k_{{\rm{sgs}}}}} + v} \right)\frac{{\partial {k_{{\rm{sgs}}}}}}{{\partial {x_j}}}} \right] - \varepsilon w} \end{array} $

(6)

该模型无需采用试验滤波、计算量小且工程应用性高，适用于超高雷诺数冷却塔结构的风场模拟。

2.2 几何建模与网格划分

几何模型采用实际尺寸建立，按30%透风率^[11]考虑正常运行状态下的百叶窗开启效应。塔筒外表面沿环向间隔12°布置一个测点，沿子午向不同高度布置16层测点，共计480个测点。由于主筒由环向对称的三角斜杆单元组成，故主筒表面奇数层测点均相对偶数层逆时针旋转6°，几何模型详细测点布置如图 1所示(其中0~29 m范围内示意点为三维模型中用以模拟塔底百叶窗30%透风率而开设的等效孔洞)。

为保证冷却塔尾流得到充分发展，数值风洞尺寸定为24D×15D×4H(顺风向X×横风向Y×高度方向Z，D为底部直径，H为塔高)，冷却塔模型位置距离出口17D，阻塞率为1.21%。为了兼顾计算效率和精度，将整个计算域划分为外围区域和局部加密区域。外围区域形状规整，采用高质量的结构化网格进行划分，局部加密区域内含冷却塔模型，采用非结构化网格进行划分。核心区最小网格尺寸为0.5 m，总网格数量约为560万个，网格最小正交质量大于0.42(要求大于0.1以上，杜绝负体积出现)。计算域及模型网格划分如图 2所示。

2.3 边界条件与参数设置

计算域入口采用速度入口，按照B类地貌设置大气边界层指数风速剖面和湍流度剖面，参考风速取10 m高度处的基本风速33.8 m/s，通过UDF文件定义上述脉动风场，如式(7)和式(8)所示。

$ {v_z} = {v_0}{\left( {z/{z_0}} \right)^\alpha } $

(7)

$ {I_z}\left( z \right) = {I_{10}}{\left( {z/10} \right)^{ - \alpha }} $

(8)

式中：v₀为冷却塔所在地区基本风速；z为计算高度与地面间的距离；I₁₀为10 m高名义湍流度，其值取为0.14；地面粗糙度系数为α，取为0.15。

计算域出口采用压力出口边界，顶部和侧面采用等效于自由滑移壁面的对称边界条件，地面及冷却塔表面采用无滑移壁面。选用LES湍流模型进行计算时启用了(增强)壁面函数模型，壁面y⁺值为37.3，可保证底层网格对数分布律成立，符合此类高雷诺数结构的工程计算要求。数值计算采用3D单精度、分离式求解器，空气风场选用不可压缩流场，模型采用2.1节所述一种新的亚格子模型，压力速度耦合方程组求解采用SIMPLEC格式，计算过程设置了网格倾斜校正以提高混合网格计算效果。在进行非定常计算之前先进行RANS的定常计算，通过瞬态化处理使初始流场达到具有合理统计特征的状态。

3 数值模拟结果及有效性验证 3.1 典型测点风压系数时程

在空气动力学中，物体表面的压力通常用量纲一压力系数C_pi表示为

$ {C_{{\rm{p}}i}} = \frac{{{p_i} - {p_\infty }}}{{{p_0} - {p_\infty }}} $

(9)

式中:C_pi为测点i处的压力系数，p_i为作用在测点i处的压力，p₀和p_∞分别为参考高度处的总压和静压。

图 3给出了冷却塔中部迎风面、侧风面及背风面3个典型区域测点压力系数时程曲线。由图可见压力系数脉动特性十分明显，时程分布不对称且带有大幅值的脉冲特性。

3.2 有效性验证

对冷却塔锥段和直筒段的平均风压系数和规范^[12-13]曲线进行了对比，结果如图 4所示。图 5给出典型截面的脉动风压系数与国内外大型冷却塔实测^[14-16]及风洞试验结果^[17]对比。由图可知：(1)下部锥段平均风压分布曲线的负压极值点和分离点对应角度与火工规范中双曲冷却塔一致，但在背风区负压数值略大于规范值，上部直筒段风压系数分布曲线与荷载规范中圆截面构筑物分布曲线基本吻合；(2)脉动风压沿环向分布趋势与国内外双曲冷却塔的实测结果较为相似，沿着环向均呈现先增大再减小的变化规律，但其数值与实测结果稍有差异，这是由于大型冷却塔脉动风压分布与所处地形、相应高度、来流湍流和周边干扰等因素均密切相关，且直筒-锥段型钢结构冷却塔自身体形与双曲冷却塔也有差异，综上认为本文基于大涡模拟得到的脉动风压具有一定的有效性，可为此类直筒-锥段型钢结构冷却塔的脉动风荷载设计取值提供参考。

4 风荷载特性 4.1 平均风压

图 6和图 7分别给出了冷却塔平均压力系数分布云图和不同高度处的风速流线图。由图 6可知，塔顶和塔底的背风区存在着比中部更加明显、范围更大的风压平台，结合图 7中的对比可以发现冷却塔不同高度处风的流动形式决定了平均风压的分布规律：(1)冷却塔下部断面的风速流线在环向70°左右的位置达到速度最大值，从而导致风压最小值提前出现。冷却塔中部断面在环向90°左右达到风速的最大值，随着断面高度的增大，风速最大值出现位置又将提前，最终达到在环向70°左右出现风速最大值。由于风速与风压存在着平方的关系，风速最大值出现区域即对应平均体形系数最大绝对值所在的环向位置；(2)冷却塔下部断面较早出现了流动分离，中部断面流动分离点的环向角度明显推迟至环向160°之后，塔顶断面流动分离重新前移至环向100°，对比压力系数云图可见流动分离的区域对应了冷却塔体形系数分离角度对应的位置。

基于模拟结果，分析冷却塔外表面各测点的平均风压系数并与规范进行对比。根据钢结构环向节点位置将冷却塔外表面测点分为：下部测点(1~4层)，上部奇数层测点(第5，7，9，11，13，15层)，上部偶数层测点(第6，8，10，12，14，16层)3个部分。图 8给出了冷却塔各层的体形系数分布与规范对比曲线，由图可见：(1)CFD数值模拟得出的冷却塔下部锥段压力系数与双曲冷却塔的数据基本吻合，锥段底部受气流影响风压略偏离于规范数据；(2)上部直筒段的数值模拟结果与圆截面构筑物的压力系数十分接近，由于端部绕流特性使得塔顶压力系数偏离较大，但沿环向的分布趋势与规范相似。

通过上述对风压分布规律进行整合，去除对冷却塔风压分布影响明显的干扰点，平均后给出冷却塔整体体形系数的推荐取值。基于最小二乘法原理，以富氏级数展开式对冷却塔整体体形系数进行拟合，当m≥7时，能取得良好的拟合效果。现给出m=7时的整体体形系数拟合公式如式(10)所示。

$ {C_{\rm{p}}} = {a_0} + \sum\limits_{i = 1}^m {{a_i}\cos \left( {i \times \theta } \right)} $

(10)

式中：a₀=-0.536 2，a₁=0.362 0，a₂=0.877 6，a₃=0.345 3，a₄=-0.067 0，a₅=-0.006 2，a₆=0.015 9，a₇=-0.018 4。图 9给出了直筒圆锥型冷却塔体形系数拟合曲线和两种规范的对比，可以发现直筒圆锥型冷却塔体形系数分布介于双曲线冷却塔和圆截面构筑物之间，吻合相对较好。

4.2 脉动风压

图 10和图 11分别给出了冷却塔瞬时流线图和涡量图，由图可知：(1)气流在冷却塔顶部前缘发生流动分离，局部出现加速效应；(2)气流发生流动分离后在冷却塔顶部和侧面产生了漩涡脱落，冷却塔背风面形成尾流涡旋以及回流，涡旋作用于冷却塔背风面，形成了吸力。

图 12和图 13分别给出了冷却塔表面脉动风压随环向角度和子午线高度变化云图和脉动风压最值沿高度分布曲线。由图可知：(1)脉动风压关于来流方向不对称分布，各层最大值分布在冷却塔侧风面50°~75°范围内；(2)最大整体脉动风压0.34位于子午向77 m高度环向72°角度处，最小0.17位于附属桁架第1~3层的背风面处；(3)直筒段的奇数层和偶数层脉动风压具有明显的差别，脉动风压沿子午向和环向没有与平均风压类似的规律性分布。

4.3 极值风压

在验算冷却塔风荷载的强度和局部稳定性时需采用体形系数极值结果进行分析，极值风压系数是基于一定的保证率得到的，故可兼顾冷却塔的安全性和经济性。通常采用峰值因子法^[18-19]计算极值风压系数，该法假定脉动风压服从高斯分布，并基于零穿越理论给出峰值因子，将其乘以脉动风压，再加上平均风压即可得到风压峰值，如式(11)所示

$ {C_{{\rm{pext}}}} = {C_{{\rm{pmean}}}} \pm g{C_{{\rm{prms}}}} $

(11)

式中：C_pext，C_pmean和C_prms分别为极值、平均和脉动压力系数；g为峰值因子，取值参考作者在文献[18]中对冷却塔结构峰值因子的研究结果，在这里取为3.5。

图 14给出了冷却塔下部锥段、上部直筒段奇数层和偶数层的极值风压沿环向和子午向变化云图。由图可知：(1)极值风压分布关于风轴方向比较对称，最大正值和负值位于冷却塔的迎风面和侧风面，最大负向极值风压集中在子午向60~125 m高度处；(2)下部锥段侧风区风压负值大小及分布面积明显小于上部直筒段，而迎风区风压正值大小及分布面积相对上部直筒段较大，背风区域极值风压分布也有很大区别。

综合下部锥段和上部直筒段的数值模拟结果并基于最小二乘法原理，以富氏级数展开式对极值风压系数分布曲线进行拟合如式(12)所示，其中n=10，参数b_k取值见表 2。图 15给出了极值风压拟合曲线与数值模拟结果的对比曲线，由图可见拟合效果很好。

$ {C_{{\rm{pext}}}}\left( \theta \right) = \sum\limits_{k = 0}^n {{b_k}\cos k\theta } $

(12)

4.4 阻力系数

塔身在顺风向所受的风压力可用塔筒阻力系数^[20] C_D来表示，如式(13)所示。

$ {C_{\rm{D}}} = \frac{{\sum\limits_{i = 1}^n {{C_{{\rm{p}}i}}{A_i}\cos {\theta _i}} }}{{{A_{\rm{T}}}}} $

(13)

式中：A_i为第i测点压力覆盖面积；θ_i为第i测点压力方向与风轴方向的夹角；A_T为结构向风轴方向的投影面积。

图 16给出了钢结构冷却塔的层阻力系数与规范对比示意图。由图可知，钢结构冷却塔上部直筒和下部锥段的结构特点使得其层阻力系数大于规范中圆截面构筑物和双曲线形冷却塔的阻力系数值，下部锥段的阻力系数随高度的增大逐渐减小，上部直筒段阻力系数变化范围相对较小，但受顶部极小风压系数的影响其层阻力系数骤减。

5 结论

基于大涡模拟方法研究了国内拟建的新型超大直筒-锥段型钢结构冷却塔流场作用和气动力分布，得到如下主要结论：

(1) 通过与国内外双曲冷却塔实测和风洞试验结果对比验证了数值方法的有效性。其中下部锥段平均风压分布曲线的负压极值点和分离点对应角度与火工规范一致，但在背风区负压数值稍大于规范值，上部直筒段平均风压分布曲线与荷载规范中圆截面构筑物基本吻合。

(2) 钢结构冷却塔风速流线最大值点和流动分离点在下部锥形断面较为靠前，随着断面高度的增大位置逐渐后移，到喉部以上又重新提前。基于最小二乘法给出了整体平均风压系数的拟合公式和分布曲线。

(3) 脉动风压关于来流方向呈现不对称分布，各层最大值均位于钢结构冷却塔侧风面50°~75°范围内，最大整体脉动风压0.34位于子午向77 m高度环向72°角度处，最小值0.17位于附属桁架第1~3层的背风面处。

(4) 极值风压分布沿风轴方向对称，最大正值和负值分别位于冷却塔的迎风面和侧风面，最大负向极值风压集中在子午向60~125 m高度处，采用富式级数形式给出了冷却塔表面极值风压拟合公式

(5) 钢结构冷却塔下部锥段的阻力系数随高度的增大逐渐减小，上部直筒段阻力系数变化范围相对较小，层阻力系数大于规范数值，但在顶部由于风压系数极小而使其阻力系数骤减。