高层建筑是现代城市风灾的主要承灾载体，相对于主体承重结构，高层建筑围护结构(尤其是玻璃幕墙)往往更容易在风灾中发生破坏^[1-2]。早在20世纪70年代，研究者们即基于Lubbock风暴的调查结果，指出应重视飞掷物对高层建筑窗户的不利影响^[3]；在飓风Alicia、Hugo和Andrew的灾后综合研究报告，也明确指出多数幕墙破坏是由飞掷物冲击造成^[4-7]。中国东南沿海是台风通过高频区，也是经济发达的城市密集区^[8]，台风York(1999)、超强台风“莫兰蒂”(2016)、强台风“天鸽”(2017)等均导致登陆城市大量幕墙受飞掷物冲击破坏^[9-11]。

飞掷物的冲击破坏问题涉及较多影响因素，比如：流场特征、不同形式的冲击物在风场中的飞行轨迹，以及冲击破坏能力等；而飞掷物飞行轨迹的确定，是开展飞掷物风灾易损性研究的基础^{[1, 12]}。20世纪80年代，研究者们进行了大量富有成效的理论与试验研究，通过系统的碎片基础性试验研究，有针对性地测量了不同形状碎片的气动力系数，建立了考虑马格努斯效应的板状碎片定常气动力模型^[13-16]，但研究中没有考虑来流湍流及流场三维特性的影响，未分析碎片的三维轨迹等。21世纪初，国外研究者们开始考虑湍流的影响^[17]，在风洞中进行了不同形式飞掷物的轨迹试验，并设计了低矮房屋周边的飞掷物模拟试验^[18-20]。国内在这方面的研究虽起步较晚，但由于近几年风致碎片导致的破坏问题越来越突出，相关理论分析和数值模拟研究工作也陆续开展^[21-24]。总体上，目前国内外学者的相关研究主要围绕低矮房屋的飞掷物运动破坏问题，研究对象主要是片状飞掷物；实际上，城市飞掷物类型更具多样性，钢化玻璃碎粒、景观砂砾等粒状物是导致大量幕墙破坏的主要碎片类型^[5]；此外，在流场影响方面，城市高层建筑群相互干扰效应导致周围流场更为复杂，飞掷物在其中的飞行轨迹与二维流场结果存在显著差异。

针对上述问题，本文基于Fluent软件平台，结合某实际高层建筑群风灾案例，通过大涡模拟方法获取建筑群周围的非定常三维绕流场，采用确定性方法分析颗粒飞掷物的三维运动轨迹。

1 研究对象概况

某塔楼群由两栋方形截面塔楼和裙房组成，如图 1所示，两塔楼的宽度和高度分别为D₁×H₁=45.6 m×116.2 m和D₂×H₂=35.4 m×89.2 m，两塔楼的水平间距L=14.4 m，裙房高度H₃=23.5 m。该塔楼群距离台风“莫兰蒂”登陆路径的垂直距离约7 km，登陆期间最大地面观测风速的风向为NNW^[10]，与该双塔楼立面朝向基本一致，图 1给出了相应的来流方向。

根据灾害现场统计，本案例在最大地面观测风速NNW风向下(图 1所示来流方向)，T2塔楼迎风面玻璃破碎率超过70%，在所有立面中破坏最严重，经初步评估，飞掷物的冲击破坏是导致大面积玻璃损坏的主要原因^[10]。由于两楼群周边为空旷场所，其周边高处产生碎片的可能性较低，因此，本研究中主要考虑在台风登陆产生最大风速的NNW风向(图 1所示)，处于上游位置的T1塔楼局部钢化玻璃碎粒产生的飞掷物颗粒，对T2塔楼迎风面玻璃幕墙的影响。

2 数值模拟方法及参数设置 2.1 飞掷物风场模拟

为研究流场对飞掷物运动轨迹的影响，本文考虑了4类流场类型：1)理想均匀流，即不考虑来流湍流度和平均风剖面；2)绕流场，即以图 1所示塔楼群为对象，进行紊流边界层风场内的三维非定常绕流大涡模拟(large eddy simulation, LES)计算得到；3)三维湍流场，即不考虑建筑影响的三维自由来流湍流场，为确保流场对飞掷物飞行轨迹影响结果的可比性，该工况的风速场来自第2类流场的建筑群三维非定常绕流大涡模拟的计算域入口处记录的脉动风速值；4)三维平均风场，即仅考虑第3类流场的平均风速剖面，不考虑湍流。表 1给出了上述4类流场的主要参数特征。

非定常绕流大涡模拟中，空气可视为粘性不可压缩流体，对瞬态的Navier-Stokes(N-S)方程进行空间平均，可得大涡模拟方法的控制方程为^[25]

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

(1)

$ \frac{\partial \bar{u}_i}{\partial t}+\frac{\partial}{\partial x_j}\left(\bar{u}_i \bar{u}_j\right)=-\frac{1}{\rho} \frac{\partial \bar{p}}{\partial x_i}+\nu \frac{\partial^2 \bar{u}_i}{\partial x_j \partial x_j}-\frac{\partial \tau_{i j}}{\partial x_j} $

(2)

式中：ρ为空气密度，ν为空气运动粘性系数，u_i、u_j表示空间平均滤波后的速度，x_i、x_j为空间坐标分量，(i，j=1, 2, 3)分别代表水平x、y轴和竖向z轴分量，p为空间平均滤波后的压力，τ_ij为N-S方程空间平均后出现的不封闭项——亚格子应力，需要建立模型来使方程组封闭。

本文首先对宽高比为1∶ 1∶ 6方形截面的单体建筑进行大涡模拟，并将数值模拟结果与风洞试验进行对比，以验证大涡模拟方法及参数设置的准确性。单体建筑数值模拟中，模型缩尺比为1∶ 500，缩尺模型高度H=0.6 m；为保证阻塞率不大于3%，以及绕流后的流场充分发展，计算域大小取为：16H(流向y)×9H(展向x)×6H(竖向z)；采用区域分块非均匀结构化网格离散计算域，近壁面处网格加密，远离壁面处采用稀疏网格，边界条件设置及局部网格划分如图 2(a)、(b)所示，采用2种网格划分方案以考察网格分辨率的影响，详细参数见表 2。

塔楼群大涡模拟中计算域大小设置为16H₁(流向y)×9H₁(展向x)×6H₁(竖向z)，采用与单体建筑相同缩尺比以及网格划分方案，最小网格尺度与单体建筑Mesh-1网格相同(下文分析中解释)，对应近壁面第一层网格无量纲距离y⁺＜5.0；网格总数592万，边界条件设置及网格划分如图 2(a)、(c)所示。

采用速度入口边界条件，大涡模拟入流脉动的合成采用基于大气边界层自保持边界条件的改进的涡合成法^[26]。出流面采用压力出口边界条件；两侧面和顶面采用对称边界条件，等价于自由滑移的壁面；计算域底面和研究对象表面采用无滑移壁面。压力速度的耦合采用SIMPLEC算法，动量方程采用Bounded Central Differencing离散格式，时间离散为二阶全隐格式，时间步长0.000 5 s，模拟总时长10.0 s，其中后9.5 s用于记录风速和风压统计计算。亚格子模型采用Dynamic Smagorinsky动态亚格子模型^[27-28]，该模型可更好地考虑钝体周围的流场的冲撞、分离、自由剪切层、涡脱落等现象。

2.2 飞掷物运动轨迹计算

飞掷物在风场中的飞行轨迹，一般将受到升力、阻力、升力矩等气动力的影响，而对于颗粒飞掷物可仅考虑阻力的影响^[13-16]，即

$ F_{\mathrm{D}}=\frac{1}{2} \rho\left(\sqrt{\sum\limits_i\left(u_i-v_i\right)^2}\right)^2 C_{\mathrm{D}} A $

(3)

式中：u_i、v_i分别表示飞掷物所在非定常风场的瞬时风速、飞掷物飞行速度; 下标(i=1, 2, 3)表示沿水平向x、y轴和竖向z轴的分量，均为随时间和位置变化的变量，其中非定常风场的瞬时风速u_i分别考虑飞掷物风场模拟的4种流场情况; $\sqrt{\sum\limits_i\left(u_i-v_i\right)^2}$表示作用于颗粒飞掷物的相对速度合矢量; C_D为飞掷物的阻力系数，钢化玻璃碎粒、景观砂砾为形状不规则的多面体，根据文献[19]研究，对于常见的正方体和球体，值分别取为0.8和0.5;A是飞掷物横截面面积。

采用牛顿第二定律，对颗粒飞掷物在3个坐标方向可建立三维非定常风场中的运动方程，整理后如下所示:

$ \frac{\mathrm{d}^2 x_i}{\mathrm{~d} t^2}=C\left(u_i-v_i\right) \sqrt{\sum\limits_i\left(u_i-v_i\right)^2}-g_i $

(4)

式中：t为时间; x_i为飞掷物的空间坐标位置; g_i为加速度分量，这里仅考虑竖向加速度分量，即水平分量g₁=g₂=0，竖向分量g₃=g(其中g为重力加速度); C表示颗粒的气动力参数与质量参数之比，即:

$ C=\frac{\rho C_{\mathrm{D}} A}{2 m} $

(5)

式中m是飞掷物颗粒质量。飞掷物运动轨迹方程(4)的求解采用五阶Runge-Kutta法。需要说明的是，参数C主要考虑颗粒物固有特性的相对关系，与无量纲化系数Tachikawa数(K数)^[13-16]存在对应关系

$ K=\frac{\rho U^2 A}{2 m g}=\frac{U^2}{g C_{\mathrm{D}}} C $

(6)

3 结果与分析 3.1 大涡模拟方法验证

为说明本文大涡模拟非定常风场的可靠性，首先将宽高比1∶ 6方形截面单体建筑大涡模拟结果与文献风洞试验^[29]结果进行了对比。

图 3为大涡模拟所得B类风场与风洞试验^[29]结果的比较，图中“LES”、“Exp”和“规范”分别表示本文大涡模拟结果、文献风洞试验^[29]和《建筑结构荷载规范(GB 50009—2012)》^[30]建议风剖面；大涡模拟风场结果来自于空计算域(其网格由图 2中建筑模型位置采用网格填充所得)中拟放置建筑模型位置处，脉动风速谱选取该位置处距离地面H高度处监测点，U_H为单体建筑高度H处的来流平均风速。由图 3可见，大涡模拟所得平均风剖面与试验^[29]和规范值^[30]具有较好的一致性；湍流度剖面与试验^[29]剖面也能较好吻合，二者均在近地面处高于规范值；大涡模拟所得脉动风速谱也基本能够反映Karman谱的分布情况，模拟风速谱在无量纲频率超过1.0后的高频处产生了一定衰减，这是由于大涡模拟离散网格的滤波效应所致。

图 4为单体建筑2/3H高度表面风压系数与试验结果^[31]对比，图中风压系数采用模型高度H处来流平均风速U_H无量纲化，风压系数均值C_{p, mean}和根方差值C_{p, rms}分别表示测点平均与脉动风压系数。

由图 4可见，Mesh-1和Mesh-2两种网格划分下，大涡模拟结果所得结构表面平均和脉动风压分布与试验结果具有较好一致性；在数值上，两种工况均与试验吻合较好，进一步加密网格Mesh-2工况在模拟精度方面并未明显改善。可见，本文大涡模拟方法及参数设置能够较好地预测结构周围的三维非定常绕流风场。

3.2 塔楼群大涡模拟结果分析

综合考虑计算精度和计算效率，塔楼群的大涡模拟计算选取上述单体建筑Mesh-1网格划分方案；来流风向为最大地面观测风速的NNW风向(图 1)。

为说明塔楼间三维非定常流动的复杂性，图 5给出了塔楼群周围距地1/2H₁高度水平截面的平均和脉动风速云图，以及瞬时涡量分布，图中U_mean和U_rms分别指平均和根方差风速值。由图 5可见，T1和T2塔楼两侧面流动分离区存在风速加速现象，在两塔楼中间位置也因“狭管效应”而引起了风速加速；两塔楼和底部裙房的顶部、侧面均存在不同尺度的瞬时分离涡，在流动分离后在塔楼群周围及尾流区相互混合掺杂，使得塔楼群周围流动情况分布复杂。可见，塔楼群周围的绕流场具有非常明显的三维特性和非定常特性，预计将对塔楼群周围颗粒飞掷物的运动轨迹产生影响。

图 6为大涡模拟所得T1塔楼迎风面和侧面的风压系数，并结合瞬时风压时序分析其偏度与峰度分布，图中风压系数均采用塔楼距地10 m高度处来流平均风速无量纲化，偏度(S)与峰度(K)分别为测点风压系数时程的三阶矩和四阶矩^[32]。总体来看，塔楼T1迎风面整体均为正压力作用，侧面和顶面为风吸力，其中侧面流动分离现象较为剧烈(图 5)，风吸力显著。具体来看，T1塔楼左侧面边角局部区域(黑圈所示位置处)瞬时压力较不利，其风压系数偏度和峰度值分别可达-0.84和4.89(图 6(c)、6(d))，按照文献[32]建议将|S|＞0.5且|K|＞3.5作为脉动风压非高斯区域的划分标准来判断，则该区域脉动风压具有非高斯特性。

在上述风压非高斯性区域，选择局部区域监测点P0和P1(距离地面高度分别为110 m和100 m)，从统计角度分析对于工程设计有意义的极值风压峰值因子和极值风压系数，并以区域监测点结果作为代表值，按照规范方法得到玻璃表面最大应力设计值，结果见表 3。表中C_{p, ext}为极值风压系数^[33-34]，q为风荷载组合设计值(此处仅考虑风荷载单工况)，σ为截面最大应力设计值，具体计算公式分别如下：

$ q=\frac{1}{2} \rho U_{\mathrm{r}}^2 C_{p, \mathrm{ext}} \times 1.4 $

(7)

$ \sigma=\frac{6 m_0 q a^2}{t^2} \eta $

(8)

式中：U_r为10 m高度来流平均风速，取40 m/s; 可变荷载分项系数取为1.4;η和m₀分别为玻璃幕墙折减系数和弯矩系数，按照《玻璃幕墙工程技术规范(JGJ 102—2003)》^[35]，分别取为0.525 3和0.044 2;a和t分别为幕墙分隔边长和幕墙玻璃厚度，按照工程情况分别为1.5 m和5 mm。

由表 3可见，P0和P1两监测点位置的幕墙表面最大应力设计值σ在90 N/mm²以上。按照规范^[35]查阅5 mm厚钢化玻璃幕墙的容许应力设计值f_g=84 N/mm²，显然σ>f_g，这说明该区域幕墙抗风安全储备不足，是风毁破坏风险点。

3.3 与飞掷物轨迹试验结果对比

钢化玻璃破碎后将产生离散碎片，《建筑用安全玻璃第2部分：钢化玻璃》(GB 15763.2—2005)对碎片有严格限定：在任何50 mm×50 mm区域内的最少碎片数在40个以上，允许有少许长条形碎片，但其长度不超过75 mm^[36]。从这一规定来看，作为常用的安全破璃，钢化玻璃破碎后不允许产生棱角特别突出的碎片，且要求均匀度较好，基本形状是介于规则正方体与球状之间的颗粒状。

文献[19]在风洞中测量了典型颗粒状飞掷物在均匀流场中的轨迹。为验证本文分析方法的可靠性，现参考其试验基本参数(来流平均风速U=40 m/s，球状和正方体颗粒飞掷物阻力系数分别取C_D=0.5和0.8，K数8.0)，按照风掷物运动轨迹计算方法模拟飞行轨迹，并与文献[19]试验结果进行对比，比较结果见图 7。

由图 7可见，随着时间的增加，两种颗粒的飞行轨迹均呈现非线性增加现象。本文计算结果与试验结果(图 7中“Lin试验”)能够较好地吻合，只是在飞行距离上略低于试验值，其中球状颗粒和正方体颗粒在两种结果下最大误差分别为13.78%、8.9%。本文计算结果与试验结果的一致性在一定程度上说明了本文运动轨迹计算方法的可靠性。

3.4 塔楼群幕墙碎粒飞掷物轨迹研究

从塔楼群大涡模拟结果知，气流遇到建筑物后发生冲撞、分离等复杂流动现象(三维非定常绕流场)，流场的时空分布与实验室模拟的理想均匀流场或规范研究中采用的特定大气边界层风场(三维自由来流湍流场)等，存在显著的区别。

为评估流场差异对颗粒物飞行轨迹的影响，本节固定颗粒飞掷物的K数为8.0，入流10 m高度处平均风速值U₁₀=40 m/s，在T2塔楼迎风面前端水平距离60 m、离地高度100 m范围内释放颗粒物，以接触至塔群墙面或裙房顶面为截止模拟时刻，对比不同流场条件颗粒全过程飞行轨迹及最终飞行速度。图 8和图 9为相应的对比结果，主要考虑了前述4种流场条件：理想均匀流场、基于大涡模拟所得三维自由湍流场(图中“三维湍流场”)、三维自由平均风速场(图中“三维平均风”)，以及塔楼群三维非定常绕流场(图中“绕流场”)。具体对比了3个典型释放点(P1、P2和P3)结果，其中P1为前文T1塔楼玻璃幕墙的风致破碎风险点，P2和P3为对比工况，目的在于说明在绕流场中颗粒物初始位置对飞掷物运动轨迹的影响。从图 8和图 9可以看出：

1) 总体上，在相同的流场条件和释放位置情况下，当接触到T2迎风面时(y=60 m)，球状颗粒(图 9)的下降高度略大于正方体颗粒(图 8)。在绕流场中，P1释放的正方体颗粒最终接触到T2迎风面的73 m高度处，比球状体高20.7 m。

2) 对比不同的流场结果来看，理想均匀流场与三维自由湍流场的模拟结果比较接近；对于三维自由来流湍流场，若忽略三维湍流成分时(图中“三维平均风”)，颗粒的飞行速度降低，飞行总距离更长。在塔楼群三维非定常绕流场(即前述第2种流场)中，释放位置影响显著，当释放位置靠近上游干扰建筑侧面时(P1位置)，颗粒物初始下落，随后向远离壁面方向偏转(图 8(c)、图 9(c))，到达下游建筑立面之前，飞行速度要大于其他2个释放位置结果。

3) P1释放点在三维绕流场中的飞行轨迹模拟结果对于飞掷物所致风灾模拟具有参考意义。从图中可看出，P1释放颗粒物运动表现为2个阶段。a)初始阶段，颗粒物刚脱离建筑近壁面，此时其飞行轨迹主要受重力控制，该阶段轨迹以竖直分量为主，阻力影响较小，飞行速度要大于其他工况(图 8(b)、图 9(b))。b) 第二阶段，颗粒物下降至一定高度后，受回流区水平风速分量的影响进入了T1塔楼的强分离流动区域，开始加速往下游飞行，此时其飞行轨迹主要受流场控制，水平运动分量为主，随后阻力与重力达到相对平衡状态，此时飞行速度变化缓慢，直至接触到T2立面。从这一规律可看出，在塔楼群玻璃幕墙的风致飞掷物破坏评估中，应考虑塔楼群周围三维非定常绕流的影响，理想均匀流或是三维湍流场中的模拟结果可能会错误估计颗粒物在建筑群中的撞击速度，并会误判建筑受飞掷物影响位置。

4 结论

对比研究了不同风场中颗粒状飞掷物的运动轨迹，并结合某一实际案例进行了高层建筑群三维非定常绕流场下飞掷物轨迹模拟，主要结论如下：

1) 本文飞掷物运动轨迹计算方法结合大涡模拟获得的流场，能够有效地对飞掷物飞行轨迹进行预测。在保持颗粒释放位置相同的情况下，与三维自由来流场相比，理想均匀流场下颗粒的飞行距离及飞行速度较大；三维自由来流场忽略来流湍流成份时会使颗粒的飞行距离增大，飞行速度降低。

2) 理想均匀流和自由来流场得到的飞行轨迹与三维非定常绕流场结果相差较大；此外，在三维非定常绕流场中，不同初始位置下的颗粒物轨迹相差较大，并可能对于最终的评估结果产生显著影响。因此总体来看，风致飞掷物的飞行轨迹预测，应充分考虑建筑群三维非定常绕流场及可能的初始破坏位置的影响。