非稳态雷暴冲击风场的瞬态数值模拟

汪之松， 王 超， 刘亚南， 李正良

(1.重庆大学 土木工程学院，重庆 400045; 2.重庆大学 山地城镇建设与新技术教育部重点实验室，重庆 400045;3.湖南省建筑设计院，长沙 410001)

非稳态雷暴冲击风场的瞬态数值模拟

汪之松1,2， 王 超3， 刘亚南1， 李正良1,2

(1.重庆大学 土木工程学院，重庆 400045; 2.重庆大学 山地城镇建设与新技术教育部重点实验室，重庆 400045;3.湖南省建筑设计院，长沙 410001)

冲击射流是雷暴冲击风场研究中最常用的模拟方法。大部分物理试验和数值模拟都是假定射流速度入口风速不随时间变化，而实际雷暴冲击风的下沉气流速度都是随时间连续变化的。在雷暴冲击风的全生命周期内一般会经历逐渐增大到最大值后再逐渐减小的过程。对于重要的工程结构抗风设计而言，得到雷暴冲击风全生命周期的风速时程信息十分有必要。基于冲击射流三维足尺模型模拟雷暴冲击风风场，在入口处引入一个更符合雷暴冲击风实际演化过程的衰减函数，采用大涡模拟数值分析获得了雷暴冲击风的非稳态风场，并得到其随时间衰减的瞬态演化过程。结果表明，该模拟方法可以较好地再现雷暴冲击风场的非稳态过程，为进一步讨论非稳态雷暴冲击风场中的结构风荷载特性奠定了基础。

雷暴冲击风；非稳态；冲击射流；大涡模拟；数值模拟

雷暴冲击风下沉气流冲击地面后，其水平风速在近地面区域产生剧烈的低空风切变，对建筑物和输电塔等结构有强烈的破坏作用。近年在国内发生了大量的输电线塔风致倒塌事故，有记录的75起输电线塔灾害事故中绝大多数都是由雷暴天气产生的雷暴冲击风造成[1-2]。国内外众多学者已采用实际观测、理论研究、物理实验和数值模拟等方法对雷暴冲击风进行了研究，并有一些学者提出了冲击风径向和竖向风剖面的解析、经验模型[3-6]，研究重点从最早的集中于实测记录和数据收集，到后来的稳态和瞬态的冲击风理论研究和实验、数值模拟。SHEHATA[7]研究了一个跨度为480 m的输电线路冲击风荷载，得出冲击风的整个生命周期和强度衰减是非常重要的参数。ABD-ELAAL[8]基于数值模拟研对非稳态雷暴冲击风的影响因素进行参数分析，表明获取雷暴冲击风的非稳态特性及其整个生命周期和强度衰减的信息十分必要。

关于冲击风的非稳态研究，FUJITA[9]记载了发生于1983年7月1日美国华盛顿圣安德鲁斯空军基地(AAFB)的一次强烈的冲击风风速时程表现出强烈的非平稳特性。HJELMFELT[10]总结了(Joint Airport Weather Studies,JAWS)研究的数次冲击风记录与出流强度随时间变化情况，表明冲击风从开始到最大强度几乎是呈线性增加的，然后再迅速衰减或者在一段时间内强度保持恒定。HOLMES等[11]通过一个衰减函数乘以径向速度，第一个提出了描述雷暴冲击风的时程曲线；CHAY等[12-13]更进一步发展了之前的概念，并且引入了几个强度衰减函数来描述风速时程曲线。CHAY等[14]利用HJELMFELT等[15]的观测结果给出了一个分析模型的强度衰减函数，他们认为强度在前5～9 min内线性增加，然后在接下来的5～9 min内呈指数衰减(如图1)。李锦华等[16]对雷暴冲击风非平稳脉动风速进行了数值合成，引入非均匀调制函数实现非平稳脉动风速的时变功率谱。

本文基于冲击射流三维足尺模型模拟雷暴冲击风场，在入口处引入一个更符合雷暴冲击风实际演化过程的衰减函数，采用大涡(LES)数值模拟获得了雷暴冲击风的非稳态风场，并得到其随时间衰减的瞬态演化过程。

图1 雷暴冲击风出流强度随时间变化的实测记录[15]Fig.1 Normalized radial velocity versus time[15]

1 CFD数值模拟

用计算流体力学软件Fluent14.0来模拟不可压缩的雷暴冲击风场，采用冲击射流模型作为雷暴冲击风的基本风场模型。为了尽可能反映真实物理尺度的雷暴冲击风现象，建立三维足尺计算域来模拟非稳态冲击射流场，模型参数包括：初始出流直径Djet=600 m，计算入口到地面的距离Hjet=2Djet，边界到喷口中心的径向距离R=6Djet，计算域示意图如图2(a)，(b)所示。

图2 计算域示意图Fig.2 The view of the computational domain

为了得到精确的数值模拟结果，采用结构化网格划分计算域，在冲击射流中心区域采用4层O型网格，在近壁面区域满足无量纲距离y+≈1，y+定义为[17]：

(1)

式中:Δy为首层网格到壁面的最小距离，单位m；ν指空气的运动黏性系数，单位m2/s；τw为壁面切应力，单位Pa；ρ为空气的密度，单位kg/m3。

计算域网格数量约为1.2×107；速度入口取湍流强度I=1%；顶部及侧边均为压力出口，压力为1.013×105Pa；近地面壁面采用增强壁面处理模型；滑移壁面剪应力取为0。

选用大涡模型(LES)对冲击射流计算域进行求解。该方法将湍流的大涡和小涡分开处理，大漩涡直接求解Navier-Stokes方程，而小涡采用湍流模型来解决。冲击射流模拟计算的时间步长为0.1 s，总迭代步数6 000步，总计算时长T=10 min。空气密度为1.25 kg/m3，黏性系数为1.789 4×10-5Ns/m2。

2 速度入口函数

实际的雷暴冲击风其速度入口出流风速随时间是连续变化的，一般有一个先逐渐增大到最大值然后再逐渐减小的过程。Abd-Elaal[14]定义了一个速度入口分段函数为：第一段为恒定值，第二段直接减小到接近于0，来描述冲击风的衰减过程。本文为更加真实地重现雷暴冲击风的衰减过程，定义依然采用两段的分段函数来表示速度入口，具体为：第一个时间段Δt1内仍然为恒定值30 m/s；第二个时间段Δt2内为呈指数逐渐减小到0，函数表达式为：

(2)

为了系统地考虑冲击风衰减周期对整个风场的影响，这里第一个时间段Δt1选取2 min，3 min和4 min三种情况，对应的第二个时间段Δt2=T-Δt1分别为8 min，7 min和6 min，速度入口速度变化函数曲线如图3所示。

图3 入口速度变化曲线Fig.3 Inlet velocityversus time

3 非稳态计算结果及分析

3.1 稳定阶段平均风速剖面

上述三种情况CFD数值模拟结果的稳定阶段(即Δt1段的后30 s风速)水平向平均风剖面与各个分析模型，经验模型和实测数据竖直风剖面和径向风剖面进行对比如下图4、5所示。CFD模拟结果选取r=1.5Djet处的水平风竖直剖面和z=0.01Djet高度处的水平风径向剖面，与Oseguera和Bowles模型[2]、Vicroy模型[3]、Wood和Kwok模型[4]、Hjelmfelt实测数据[10]和Letchford和Illidge[18]的物理实验模型等结果进行对比，其中z为竖向高度，zm为最大水平风速对应的竖向高度，u为沿径向的水平风速，um为沿径向水平风速最大值，r为相对喷口中心的径向距离，rm为最大水平风速对应的径向距离。

图4 CFD数值模拟与各模型竖直风剖面比较Fig.4 Comparison of vertical wind profile of CFD numerical simulation and analytical models

图5 CFD数值模拟与各模型径向风剖面比较Fig.5 Comparison of radial wind profile of CFD numerical simulation and analytical models

从图中可以看出各个分析模型，经验模型和实测数据与本文CFD模拟结果都给出了雷暴冲击风水平风竖直剖面和径向剖面的主要特征。在水平风竖直剖面方面，由于本文是基于冲击射流模型的模拟结果，所以和基于冲击射流模型的物理实验得出的Wood和Kwok模型基本一致。在水平风径向剖面方面，CFD数值模拟结果对达到最大水平风速后风速衰减段的风速值过高估计，主要原因是实际的冲击风经历形成、持续并最终扩散的过程，并没有停留在稳定的状态，CFD数值模拟是提取的风场稳定阶段的数据。总体来说，基于LES瞬态风场数值模拟平均水平风剖面同经验风剖面及实测数据吻合较好，在缺乏实测数据的情况下，除了风洞试验可以提供较近似的结果外，基于LES的雷暴冲击风场数值模拟有利于了解其发展全过程尤其是达到最大水平风速前流场的特性。

3.2 流场风速云图和矢量场

选取Δt1=3 min时的数值分析结果，冲击射流下沉，发展，稳定和衰减阶段的速度云图和矢量图如图6(a)～(f)所示。从图中可以看到，在t=1 min时，速度入口下沉气流与周围气体之间拖拽卷吸，在下沉气流周围形成环形漩涡(a)；在t=2 min时，入口速度保持30 m/s不变，环形漩涡冲击地面后沿径向水平发展，风场迅速向周围辐射扩散，在近地面达到最大出流强度(b)；在3 min过后，入口速度开始减小，到3.5 min时，入口速度下降到最初的50%，整个风场从入口位置开始逐渐衰减，近地面出流强度逐渐减小(c),(d)；从4.5 min到5.5 min，入口速度迅速减小到接近为零，此时整个风场都迅速衰减，随着近地面出流强度的进一步减弱，雷暴冲击风逐渐消退(e),(f)。

图6 冲击射流各时刻风度云图和矢量场Fig.6 The contour and vector plot of impinging jet at different time

3.3 水平风速的瞬态分布

选取Δt1=3 min时的数值分析结果来考察雷暴冲击风整个过程水平风速的瞬态分布情况。

3.3.1 竖直分布

雷暴冲击风水平风速的竖直分布表征了水平风速沿高度的变化规律，对结构风工程有重要参考价值。为了细致的考察冲击风整个过程，图7给出了冲击射流发展过程中不同径向位置上水平风速的竖直分布情况。从图7中可以看到在雷暴冲击风冲击地面前后和衰减过程中不同径向位置水平风速的竖直分布呈现较大差异。随着下沉气流冲击地面并沿地面像四周扩展，水平环形涡流一次经过距离冲击中心由近至远的各个径向位置，导致超过一定高度后出现负向的速度。总体而言，各个径向位置的最大水平风速在2～3 min左右迅速增加到最大值，然后随着时间增加逐渐减小，到5 min左右以后在整个竖向位置都下降到最大风速的40%以内。

3.3.2 径向分布

雷暴冲击风水平风速的径向分布表征了水平风速沿径向的变化规律，图8给出了冲击射流发展过程中不同高度位置上水平风速径向分布情况。从图8中可以看到，雷暴冲击风冲击地面后在结构风工程关心的高度范围内水平风速的径向分布具有相似的规律，最大水平风速出现在2 min左右，位于r=1.5Djet位置附近。然后随着时间的增加，水平风速逐渐减小，到5 min左右以后在整个径向位置都下降到最大风速的40%以下，这点跟竖直分布所观察到的一样，说明雷暴冲击风在最大风速衰减2 min左右整个风场水平风速都迅速减小到最大值的40%以下。

3.4 竖向风速的瞬态分布

选取Δt1=3 min时的数值分析结果来考察雷暴冲击风整个过程竖向风速的瞬态分布情况。

3.4.1 竖直分布

雷暴冲击风竖向风速的竖直分布表征了竖向风速沿高度的变化规律，为了细致的考察冲击风整个过程，图9给出了冲击射流发展过程中不同径向位置上竖向风速的竖直分布情况。从图9中可以看到在雷暴冲击风冲击地面前后和衰减过程中不同径向位置竖向风速的竖直分布呈现较大差异。随着下沉气流冲击地面并沿地面像四周扩展，水平环形涡流一次经过距离冲击中心由近至远的各个径向位置，在水平环形涡流经过的径向位置，无论是风速值还是风向，竖向风速的竖直分布呈现出较大的变化。总体而言，靠近喷口中心的r=0.6Djet在1 min内就很快达到了最大值；r=1.0Djet和r=1.5Djet位置都是在2 min内达到最大值，然后随着时间增加而逐渐减小；而离喷口中心较远的r=2.0Djet位置各个时刻的竖向风速都比较小，基本都在0.2Vjet范围以内，由此可知，竖向风速的影响范围主要到达1.5Djet左右。

图7 冲击射流全过程各径向位置水平风速的竖直分布Fig.7 Vertical profiles of radial velocity at different time and different radial distances

图8 冲击射流过程中各高度位置水平风速的径向分布Fig.8 Radial profiles of radial velocity at different time and different heights

3.4.2 径向分布

雷暴冲击风竖向风速的径向分布表征了竖向风速沿径向的变化规律，图10给出了冲击射流发展过程中不同高度位置上竖向风速径向分布情况。从图10中可以看到，雷暴冲击风冲击地面后在结构风工程关心的高度范围内竖向风速的径向分布形状上具有相似的规律，最大竖向风速出现在2 min左右，位于r=1.5Djet位置附近，但最大值的峰值大小有较大差异。然后随着时间的增加，竖向风速逐渐减小，到5 min左右在整个径向位置的风速值都下降到0.1Vjet以下。

3.5 水平速度和竖向速度时变特性

数值模拟结果水平速度在Δt1=2,3,4 min三种情况下的时程曲线如图11所示(径向位置r=1.0Djet，高度位置z=0.03Djet处)。可以看到三种情况下都有环形涡流经过，导致在2 min左右达到了出流的最大强度，随着Δt1的增加出流强度没有更进一步的增加。当入口速度开始下降的时候，监测点速度并没有突然的下降，而是保持了一段时间的稳定值然后再开始下降进入衰减阶段。

图9 冲击射流全过程各径向位置竖向风速的竖直分布Fig.9 Vertical profiles of vertical velocity at different time and different radial distances

图10 冲击射流过程中各高度位置竖向风速的径向分布Fig.10 Radial profiles of vertical velocity at different time and different heights

选取Δt1=3 min时的数值分析结果来考察雷暴冲击风整个过程风速的时变特性。在径向r=1.0Djet处结构工程比较关心的高度范围内(如z=0.002、0.005、0.025、0.1Djet)的水平风速和竖向风速时程曲线如图12、13所示；在高度z=0.01Djet处各个径向位置(如r=1.0、1.5、2.0、2.5、3.0Djet)的水平风速和竖向风速时程曲线如图14和15所示。

图12展示了速度强度时程曲线随高度的改变，速度强度刚开始随着高度的增加而增大，直到最大速度所在高度位置然后减小，其变化趋势和稳态模型水平风速的竖直风剖面的变化趋势一致；图13展示了竖向风速时程曲线随高度的改变，速度强度随高度连续的增加，其变化趋势和稳态模型竖向风速的竖直风剖面的变化趋势也一致；图14、15展示了速度强度时程曲线随径向距离的改变，速度强度随径向距离的增加而增大，到r=1.5Djet处达到最大值然后减小，其变化趋势与稳态模型水平风速的径向风剖面也基本一致。

图11 数值模拟时程曲线Fig．11Temporalprofilesoftheradialspeedforthreedifferentages图12 r=1．0Djet各高度位置水平风速时程曲线Fig．12Temporalradialspeedatr=1．0Djetfordifferentheights图13 r=1．0Djet各高度位置竖向风速时程曲线Fig．13Temporalverticalspeedatr=1．0Djetfordifferentheights

图14 z=0.01Djet处各径向位置水平风速时程曲线Fig.14 Temporal radial speed at z=0.01Djetfor different radical locations

图15 z=0.01Djet处各径向位置竖向风速时程曲线Fig.15 Temporal vertical speed at z=0.01Djet for different radical locations

4 结 论

利用大涡数值模拟方法，基于冲击射流模型完成了三维足尺的雷暴冲击风非稳态模拟。在CFD计算结果与相关经验、分析模型和实测数据进行比较以验证CFD数值模型及参数准确性的基础上，进行了风场特性的全面研究。在速度入口处引入一个更符合雷暴冲击风演化过程的速度衰减函数，得到了冲击风场整个生命周期及衰减过程的详细信息。结构风工程关心的冲击风近地面风场的水平风速和竖向风速总体而言都在短时间内迅速达到最大值，随后逐渐衰减或保持一个恒定值波动后再迅速衰减至零。基于冲击射流模型开展的CFD瞬态数值模拟可以较好地体现雷暴冲击风整个衰减过程的基本特性，为进一步讨论非稳态雷暴冲击风场中的结构风荷载奠定了基础。

[ 1 ] 张勇. 输电线路风灾防御的现状与对策[J]. 华东电力, 2006,34(3):28-31. ZHANG Yong. Status quo of wind hazard prevention for transmission lines and countermeasures [J]. East China Electric Power, 2006,34(3):28-31.

[ 2 ] 瞿伟廉，吉柏锋. 下击暴流的形成与扩散及其对输电线塔的灾害作用[M]. 北京：科学出版社，2013.

[ 3 ] OSEGUERA R M, BOWLES R L. A simple analytic 3-dimensional downburst model based on boundary layer stagnation inflow [R].NASA Technical Memorandum 100632, 1988.

[ 4 ] VICROY D D.A simple, analytical, asymmetric microburst model for downdraft estimation [R]. NASA Technical Memorandum 104053, 1991.

[ 5 ] WOOD G S, KWOK K C S, An empirically derived estimate for the mean velocity profile of a thunderstorm downbursts [C]//Proceedings of 7th AWES Workshop. Auckland, 1988.

[ 6 ] HOLMES J D, OLIVER S E. An empirical model of a downburst [J].Engineering Structures, 2000, 22: 1167-1172.

[ 7 ] SHEHATA A Y, EL DAMATTY A A, Savory E. Finite element modeling of transmission line under downburst wind loading [J]. Finite Elements in Analysis and Design, 2005, 42(1):71-89.

[ 8 ] ABD-ELAAL E, MILLS J E, MA X. A coupled parametric-CFD study for determining ages of downbursts through investigation of different field parameters [J]. Journal of Wind Engineering and Industrial Aerodynamics, 2013, 123: 30-42.

[ 9 ] FUJITA T T. The downburst report of Projects NIMROD and JAWS [R]. Chicago: University of Chicago, 1985.

[10] HJELMFELT M R. Structure and life circle of microburst outinflows observed in Colorado[J]. J. ApplMeteorol, 1988, 27(8): 900-927.

[11] HOLMES J D, OLIVER S E. An Emprical Model of a Downburst[J].Engineering Structures, 2000, 22:1167-1172.

[12] CHAY M T, ALBERMANI F, WILSON R. Numerical and analytical simulation of downburst wind loads [J]. Engineering Structures, 2006, 28(2): 240-254.

[13] ABD-ELAAL E, MA X, MILLS J E. A newly developed analytical model of transient downburst wind loads [C]//Proceedingsof the 22th Australian Conference on the Mechanics of Structures and Material, Sydney, Australia, 2012: 425-430.

[14] ABD-ELAAL E, MILLS J E, MA X. Empirical models for predicting unsteady-state downburst wind speeds [J]. J. Wind Eng. Ind. Aerodyn, 2014, 129: 49-63.

[15] WILSON J W, ROBERTS R D, KESSINGER C, et al. Microburst wind structure and evaluation of Doppler radar for airport wind shear detection [J]. Journal of Climate and Applied Meteorology, 1984, 23(6):898-915.

[16] 李锦华，吴春鹏，陈水生.下击暴流非平稳脉动风速数值模拟[J]. 振动与冲击，2014, 33(14):54-60. LI Jinhua，WU Chunpeng，CHEN Shuisheng. Simulation of non-stationary fluctuating wind velocity in downburst[J]. Journal of Vibration and Shock,2014, 33(14):54-60.

[17] Fluent Inc. FLUENT User’s Guide[Z]. Fluent Inc., 2003.

[18] LETCHFORD C W, ILLIDGE G. Turbulence and topographic effects in simulated thunderstorm downdrafts by wind tunnel jet[C]//Proceedings of the Tenth International Conference on Wind Engineering, Denmark, June 1999:1907-1912.

Transient numerical simulation for unsteady-state downburst winds

WANG Zhisong1,2,WANG Chao3, LIU Yanan1, LI Zhengliang1,2

(1. School of Civil Engineering, Chongqing University, Chongqing 400045, China;2. Key Laboratory of New Technology for Construction of Cities in Mountain Area, Ministry of Education, Chongqing University, Chongqing 400045, China;3. Hunan Provincial Architectural Design Institute, Changsha 410001, China)

Impinging jet is one of the most commonly used simulation methods in studying the wind field of thunderstorms. Most of the previous physical tests and numerical simulations assume that the inlet velocity of the jet does not change with time, while the sinking velocity of downburst winds is continuously changing with time actually during its whole lifecycle, it gradually increases to the maximum value and then decreases in general. For wind-resistant design of significant engineering structures, it’s necessary to acquire the temporal profiles of downburst wind speed in its whole lifecycle. Here, a three-dimensional full-scale impinging jet was taken as a wind field model of downburst. The decay functions for the inlet jet velocity were given for the actual downburst，and numerical simulation was performed by adopting the large eddy simulation(LES) to acquire an unsteady-state wind field where the intensity decay information of downburst was presented. The simulation results showed that the unsteady wind field of downburst can be well simulated with this proposed method, and the results are correlated well to Andrews AFB recorded ones. These results laid a foundation for further studying wind load characteristics of unsteady downburst.

downburst; unsteady-state; impinging jet; LES; numerical simulation

国家自然科学基金(51208537)

2015-10-09 修改稿收到日期：2016-01-27

汪之松 男，副教授，1980年生

TH212；TH213.3

A

10.13465/j.cnki.jvs.2017.03.009