考虑地表粗糙影响的下击暴流数值仿真

2018-11-09罗啸宇谈正铠吴晓波刘震卿涂元刚

土木工程与管理学报 2018年5期

罗啸宇，谈正铠，吴晓波，刘震卿，涂元刚

(1. 广东电网有限责任公司电力科学研究院，广东广州 510062；2. 华中科技大学土木工程与力学学院，湖北武汉 430074)

随着观测技术水平的提高，人们逐渐识别出一大部分气象灾难的罪魁祸首是下击暴流。下击暴流是气流迅速下降冲击地面后，形成的从冲击中心向四周沿地表方向辐散的一种极具破坏性的强风。下击暴流在航空航天和建筑结构两个领域对人类生产生活有重大威胁，严重情况可能导致灾难。

目前，学者主要采用三种研究方法对其展开研究，分别是基于历史记录的气象数据的风场实测研究、基于下击暴流产生和发展原理的气象理论分析、以及依靠专用的计算机进行实验室仿真模拟。在实验方面，Mason等[1]开发并采用了下击暴流模拟器，重点研究了下击暴流风时程特性。仿真模拟方面，Aboshosha等[2]采用大涡模型研究了下击暴流风场特性，着重研究了下击暴流的积分尺度、风谱特性等。翟伟廉等[3,4]都采用了二维计算域进行数值模拟，在计算域的边界条件限制了流场发展。陈波等[5]采用三维计算域进行数值模拟，但是只重点研究了风场对于平屋面的风荷载特性。Aboshosha等[2]进行实际尺寸的数值模拟但是计算资源消耗过于巨大。因此，本文进行了实验室缩尺的下击暴流风场仿真模拟，以探明下击暴流风场特性以及地表粗糙程度对风场的影响。

1 下击暴流风场模型

目前，下击暴流仿真模拟的风场模型主要归结为两种，一种是环形涡流模型(图1a)，另一种是冲击射流模型(图1b)。冲击射流模型是将射流垂直喷射向地表，气流下沉过程中同时也产生环形涡流，协同效果下模拟下击暴流的风场特性效果良好[6,7]，弥补了环形涡流模型对风场内核以及外围区域模拟效果不理想的缺点。因此，本文选取冲击射流风场模型进行研究。

图1 下击暴流风场模型

2 数值模型

2.1 控制方程

本文选用大涡模拟的方法进行仿真模拟。如上文所述，仿真模拟中直接模拟计算大尺度涡旋，对于低于网格解析度的部分小尺度涡旋采用亚格子模型考虑其对大尺度涡旋的影响，采用Boussinesq假说以及标准Smagorinsky-Lilly模型计算亚格子应力(Sub-Grid Scale，SGS)。过滤亚格子涡旋后的时间依赖Navier-Stokes方程(笛卡尔坐标系)如下所示：

(1)

(2)

(3)

(4)

(5)

(6)

为了与光滑地面情况的风场作对比，本文选取了浸没边界法(Immersed Boundary Method，IBM)来实现对存在建筑物的地面区域的模拟。这种方法直接在Navier-Stokes方程中添加地面粗糙引起的阻力项，进而可以有效避免建立粗糙元的几何外形的步骤，公式如下：

(7)

(8)

空间离散采用有限体积法，二阶中心差分格式用于对流项与粘性项，二阶隐式格式用于非稳态项的时间推进，SIMPLE算法用于压强速度解耦，见文献[9]，求解器采用Fluent 6.3.26[10]。

2.2 几何模型

选取一个完整的三维圆柱体作为计算域，如图2所示。这使流场演变的过程不会受到限制，模拟结果更为精确。选取几何缩尺为1∶1000，取冲击射流的直径Dj为1 m。所选用的圆柱体计算域的尺寸如图2a所示，计算域直径为8Dj，竖直方向高度为4Dj。其中冲击射流入口位于竖直高度2Dj处。冲击射流通道也为圆柱体，延伸至圆柱体计算域上部圆面。根据下击暴流的风速特征以及选取的计算域尺寸，本文采用Mason等[11]推荐的速度缩尺1∶3，冲击射流的速度取为20 m/s。对于粗糙地表条件，动量平衡方程中的阻力项将加入到距离地表10 cm的范围内。

图2 大涡模拟计算域及其网格分布

2.3 网格分布与边界条件

使用GAMBIT划分三维结构化网格，网格数量为594000。在网格划分过程中，首先在水平截面中心按照正方形划分网格，再在圆周方向上划分出面网格，最后完成体网格的划分，以确保网格结构质量良好。贴地网格高度为0.1 mm，整个计算域最大网格高度为10 cm，最小水平网格尺寸为1 mm，最大水平网格尺寸为10 cm，竖直与水平网格增长率均小于1.05。冲击射流通道与计算域相交的圆柱面为无剪切应力滑移墙面，排除冲击射流在进入计算域形成下击暴流前对其的影响。圆柱体计算域中，底面为壁面边界条件，圆柱面及上表面为压力出口。计算时间步长为0.0001 s，总计算时间为50 s。

3 光滑地表流场特性

3.1 数据对比

将本文数值模拟结果与以往实测风场以及仿真模拟结果以最大风速做归一化处理，绘制风速剖面并进行对比，见图3(图中，z为高度；ur为径向速度；urpeak为最大径向速度)。可以清楚地看出本数值模拟与以往的研究结果[1,2,4,12,13]吻合良好，由此验证了本文湍流模型、网格分布、边界条件以及时间步长等设置的可靠性。

图3 风速对比

3.2 风场随时间的演变

由于风场的对称性，本文选取一个竖向截面进行研究。图4(图中Tn为归一化时间)为所选取的特征时刻下击暴流风场涡量云图，以此可以清楚地展现出流场演变过程。图4a中，冲击射流产生的主环形涡流开始冲击地面，地面边界处已经产生次生涡流；图4b中环形涡流正冲击地面，地面边界处产生次生涡流；图4c中，主环形涡流脱落，在径向上产生位移；图4d中，主环形涡流径向移动；图4e中，主环形涡流径向移动较大，并竖直向上移动；图4f为风场末期，涡流逐渐消散。

图4 风场涡量云图

3.3 风场分析

为了降低坐标系对数据的影响，将速度分解为轴向速度和径向速度两部分，并且仅着重分析变化情况复杂的径向速度特征，着眼四个特征位置处的风速剖面特征进行分析。在径向风速变化最为复杂的半径1Dj～2.5Dj区域中选取了四个特征位置绘制风速剖面图，如图5，半径R分别为1Dj，1.5Dj，2Dj，2.5Dj，而绘制风剖面图所选取的时间点是四个位置处分别达到最大径向速度的时间点，Tn分别为7.8，8.8，11，12.8对应上文中相同时间点的涡量图。可以清楚地看出当R=1Dj达到最大径向风速时，主环形涡流中心已经过该位置，但涡流并未完全离开该位置。此时剖面处上部接近涡流中心，因此风速较低；剖面下部受主环形涡流冲击地面的后续涡流作用，速度较大。在R=1.5Dj处达到最大径向速度时，剖面正好经过涡流中心，这也正是剖面上端风速出现负值的原因。剖面下端存在次生涡流，其附近风速增强到1.6Vj(Vj为Dj处风速)。下击暴流冲击地面后，会产生沿地面边界流动的气流。该气流受主环形涡流的负压力影响，剧烈变动后，形成了与主涡流方向相反的次生涡流。主环形涡流与次生涡流协同作用下，近地面风速得到较大增强。而R=2Dj处达到最大径向速度时，主环形涡流发生了较明显的上移，同时次生涡流强度逐渐减弱，因此该处最大风速低于R=1.5Dj处。同时，主环形涡流的上移也导致未被次生涡流增强的位置处风速略小，也就表现为图中风速曲线凹陷部分。而R=1Dj位置处，风速又得到增强是副环形涡流冲击地面的结果。由于R=2Dj处距风场中心较远，能量耗散，最大风速也较低，次生涡流明显消散，因此近地面风速没被明显增强。

图5 最大径向风速剖面

4 地表粗糙影响

图6为统计了模拟的全部时刻0～0.3Dj高度内的风速后绘制的最大风速包络图。可以清楚地看出地面粗糙对最大径向风速的影响。不加设地面粗糙时，最大径向风速达到1.8Vj并且发生在0.25Dj高度处；加设地面粗糙后，近地面风速减小至0.2Vj，并在近地面范围内出现陡坡，然后风速随高度增加迅速增大至1.4Vj，在0.1Dj～0.18Dj之间分布均匀。结合模型建立的部分，可以分析出速度陡坡是由于地表加设的粗糙层引起的，其高度与粗糙层高度一致。径向风速达到最大值后，无论地表粗糙与否，径向风速沿高度的衰减率基本一致。综上，在加大地表粗糙后，最大径向风速减小，并且出现的位置上移。这与Mason等[1]采用方法物理模拟的实验结果一致。

图6 最大径向速度包络图

图7是统计竖向0≤Z≤0.2Dj、径向0.5Dj≤R≤3Dj范围内各点最大径向速度后绘制的等高线图。图7a中出现了两个峰值，分别位于径向0.9Dj～1.0Dj和1.2Dj～1.6Dj之间，并且前者峰值分布范围小于后者。结合前文涡量图可以对此做出解释，两峰之间即为主环形涡流冲击地面的中心，因此较大的径向风速出现在该位置两侧；后者分布范围更大的原因正是次生涡流的增强作用以及涡流的竖向移动。加设粗糙情况下最大径向风速整体上低于光滑情况，这与前文结论一致，峰值分布在1.1Dj～1.4Dj和0.07Dj～0.18Dj之间，呈狭长形状。如此不同的分布结果是两种地面情况中主环形涡流的移动路径不同引起的。