叶形屋面风致积雪重分布及干扰特性研究

Study on Redistribution and Interference Characteristics of Snow Caused by Wind on Leaf-Shaped Roofs

摘要

叶形屋面造型新颖多变，近年在大跨空间建筑中得到应用推广，但复杂的屋面造型给其雪荷载的不均匀分布取值带来一定挑战。本文基于Euler-Euler方法和风雪单相耦合，借助数值风洞模拟了不同风向下多叶形建筑屋面风致雪漂移的积雪不均匀分布。研究结果表明，叶形屋面风致积雪分布不仅要考虑屋面的倾斜，且应合理评估凹凸对积雪分布的影响。叶形屋面积雪通常表现为侵蚀为主，且凸起部分和叶片两端的侵蚀更为明显，同时下凹部分可能形成局部积雪的严重堆积，因此建议工程设计中屋面凹陷率不宜超过10%。叶形屋面之间的相互干扰对积雪有影响，低屋盖干扰影响主要集中于受遮挡作用屋面的迎风边缘区域，高屋盖干扰主要影响整体积雪分布系数，通常干扰变化可达到20%左右。

Abstract

Leaf-shaped roofs， known for their novel and varied forms， have been widely used in long-span spatial structures in recent years. However， their complex geometry poses certain challenges to the non-uniform distribution of snow loads. Using the Euler-Euler method and a single-phase wind-snow coupling approach， the wind-induced snow drift on multi-leaf building roofs under different wind directions was simulated with a numerical wind tunnel. The findings indicate that the wind-driven snow distribution on leaf-shaped roofs should consider not only the roof inclination but also the effect of roof surface irregularity. Snow on leaf-shaped roofs is predominantly subject to erosion， which is more pronounced on protruding portions and leaf tips. Meanwhile， concave areas may experience severe local snow accumulation； therefore， excessive concavity should be avoided in design. The snow distribution on multi-leaf roofs is significantly influenced by mutual interference among the leaf units. Such interference impacts snow accumulation patterns： interference from lower roofs mainly concentrates on the windward edge of the sheltered roof， whereas interference from higher roofs primarily affects the overall snow distribution coefficient， typically inducing variations of about 20%.

Graphical abstract

关键词

叶形屋面 / 风致雪漂移 / 积雪重分布 / 数值模拟

Key words

leaf shaped roof / wind induced snow drift / snow redistribution / numerical simulation

0 引 言

随着全球气候变暖，极端天气事件的频率和严重性显著增加，积雪灾害也在逐年加剧。例如2023年河南信阳体育馆、湖北长江大学食堂、黑龙江中学体育场等多起建筑倒塌事件，均与积雪灾害密切相关。降雪往往伴随着大风，积雪在风作用下导致的不均匀分布往往是灾害发生的主要原因。

积雪分布常用的研究方法有实地观测、风洞实验和数值风洞模拟。DeGaetano等^［1］结合雪堆指标的风速、温度以及雪水当量、降雪量和雪深等每日观测数据，得到建筑屋顶极端雪堆荷载的气候测量值。Zhou等^［2-3］分别使用低密度锯木灰、中密度泡沫塑料和高密度硅砂通过风洞试验模拟了阶梯式平屋顶上雪荷载的再分配。Zhang等^［4］采用细硅砂代替雪颗粒结合风洞实验研究了高层建筑对平屋顶积雪分布的影响。近年来随着计算机性能的显著提升和数值模拟方法理论的日趋完善，计算流体动力学（CFD）逐渐成为风致积雪模拟的重要工具之一。Beyers等^［5］将雪视为连续相，提出了单方程雪飘模型，结合有限体积法和侵蚀沉积公式模拟了格子结构的不均匀积雪分布，并与南极洲实测结果进行了对比验证。在对屋顶积雪分布的研究中，Choi等^［6］结合人工降雪实验，通过离散元方法对温室拱形屋面进行了模拟，结果表明背风面积雪较迎风面大，且随风速增加而减少。Sun等^［7］提出了一种考虑颗粒间特性的修正方法，并研究了颗粒间特性修正的雪廓线对屋顶周围流场和曲面屋顶积雪特性的影响。Zhang等^［8］研究了阶梯形平屋顶上的雪积累过程和雪荷载形状系数，并提出了一种计算雪荷载形状系数和下屋顶雪保护最佳坡度的方程。康路阳等^［9］通过引入积雪休止角改进了雪漂移模拟方法，并应用于平屋面屋面积雪重分布模拟研究。孙晓颖、洪财滨等^［10-11］通过对四种不同湍流模型在防雪栏背风区域积雪分布中的比较，发现k-kl-ω模型在模拟防雪栏背风区域积雪分布方面较优。李恺轩等^［12］充分比较了各国规范对漂移雪荷载的取值，并结合结构的减灾设计提出了相应的建议。

已有研究关注对象通常表现为结构在空间水平尺寸上相差不大，例如体育场、体育馆等结构形式。近几年随着建筑美学的要求，建筑外形更趋新颖个性化，叶形屋面逐渐应用在一些工程中。相较传统屋面形式，叶形屋面主要表现为两个方向的尺度通常相差较大，屋面可能同时存在上凸、下凹及翘起等多个复杂的外形要素，给该类屋面积雪分布的取值带来一定困难。现有叶形屋面的积雪分布研究较为缺失，特别是考虑屋面间干扰性方面的研究。

本文以某文化中心为工程背景，基于Euler-Euler两相流理论进行数值模拟。该理论是一种多相流模型，通常用于描述在同一控制空间中多种不可互溶流体之间的相互作用。本文对叶形屋面的大尺度风致雪飘移进行了数值模拟研究。叶形屋面具有显著的大尺度结构和复杂的几何特征，因此采用在CFD应用中计算效率高、收敛性好且适于大尺度湍流模拟的Realizable k-ε湍流模型。在计算流体软件Fluent基础上，结合UDF进行二次开发对叶形屋面的风致雪飘移进行数值模拟研究，对屋面风致积雪重分布及其相互干扰性进行分析探讨。

1 计算方法

本文根据Euler-Euler两相流理论，认为空气相和雪项为连续相，忽略雪对空气相影响，认为两相为单项耦合，即雪相在空气相作用下发生飘移，但雪相并不影响空气相运动^［9］，雪相采用Beyers的单方程雪飘模型^［5］。空气相控制方程通常考虑连续方程和质量方程^［10］，受篇幅有限此处省略。

1.1　雪相控制方程

雪相控制方程为

\frac{\partial {\rho }_{\mathrm{s}}f}{\partial t}+\frac{\partial {\rho }_{\mathrm{s}}f{u}_j}{\partial x_j}=\frac{\partial }{\partial x_j}\left({\mu }_{\mathrm{t}}\frac{\partial {\rho }_{\mathrm{s}}f}{\partial x_j}\right)+\frac{\partial }{\partial x_j}(-{\rho }_{\mathrm{s}}f{u}_{\mathrm{R}})（1）

式中：{\rho }_{\mathrm{s}}为积雪密度，本文取\mathrm{150}\mathrm{k}\mathrm{g}/{\mathrm{m}}^{\mathrm{3}}；f为单位体积内雪相所占比例；u_j为风的速度矢量；{\mu }_{\mathrm{t}}为空气相的湍流黏性系数；u_{\mathrm{R}}为雪相相对于空气相的相对运动，本文取\mathrm{0.3}\mathrm{m}/\mathrm{s}^［10］。

在求解空气相控制方程的同时计算雪相控制方程，计算出单位体积内雪相f所占比例。

1.2　积雪侵蚀与沉积

积雪的沉积与侵蚀根据Beyers的经验公式来看主要和摩擦速度u_{\mathrm{*}}、积雪阙值摩擦速度u_{\mathrm{*}t}有关^［5］。摩擦速度由空气相流动特性决定，计算公式如下：

u_{\mathrm{*}}=\sqrt[]{{\tau }_{\mathrm{0}}/{\rho }_{\mathrm{a}}}（2）

其中{\rho }_{\mathrm{a}}为空气密度，{\tau }_{\mathrm{0}}为壁面剪切应力。

积雪阙值摩擦速度u_{\mathrm{*}t}受温度、湿度、雪颗粒半径等多种因素影响，本文取0.2 m/s^［13］。当摩擦速度大于阙值摩擦速度时，表面积雪被刮起进入计算域，从而表现为雪侵蚀；当摩擦速度小于阙值摩擦速度时，计算域中的雪颗粒沉积在雪表面而表现为雪沉积。雪的侵蚀量q^{\mathrm{e}\mathrm{r}\mathrm{o}}和沉积量q^{\mathrm{d}\mathrm{e}\mathrm{p}}公式如下：

q^{\mathrm{e}\mathrm{r}\mathrm{o}}=A_{\mathrm{e}\mathrm{r}\mathrm{o}}(u_{\mathrm{*}t}^{\mathrm{2}}-u_{\mathrm{*}}^{\mathrm{2}})u_{\mathrm{*}t}\le u_{\mathrm{*}}（3）

q^{\mathrm{d}\mathrm{e}\mathrm{p}}=C\frac{(u_{\mathrm{*}t}^{\mathrm{2}}-u_{\mathrm{*}}^{\mathrm{2}})}{u_{\mathrm{*}t}^{\mathrm{2}}}u_{\mathrm{*}t}>u_{\mathrm{*}}（4）

其中A_{ero}为常数，取\mathrm{7.0}\times {\mathrm{10}}^{-\mathrm{4}}；C为雪的质量浓度，等于雪密度和体积分数的乘积。

积雪深度变化为

\mathrm{\Delta }h=\frac{q^{\mathrm{e}\mathrm{r}\mathrm{o}}+q^{\mathrm{d}\mathrm{e}\mathrm{p}}}{{\rho }_{\mathrm{s}}}\mathrm{\Delta }t（5）

其中\mathrm{\Delta }h为积雪变化高度，\mathrm{\Delta }t为积雪作用时间。

2 数值模拟

2.1　叶片布置及几何尺寸

叶形屋面共有8片，其空间布置及编号示意如图1所示。风向角定义以正北方向为0^o，顺时针方向旋转。由于叶形屋面形状较为复杂，可能包含上凸、下凹及翘起等外形特征，首先将屋面按照叶形的长短轴分别定义为X、Y方向，然后沿X、Y方向定义叶形的两个主轴方向断面，如图2所示。其中B对应叶形截面宽度，D对应截面两端的高差，决定了叶形的翘起程度，H_max、H_min分别对应截面表面最大上凸和下凹的程度。屋面几何尺寸见表1。

将叶形屋面建筑群的数字模型置入数值风洞中，主要保证数值风洞的四周边界离建筑足够远从而不产生明显的边界干扰效应，同时风洞阻塞率控制在3%以下。流场域的网格划分首先保证在屋面表面区域进行足够的细分加密，边界层网格距离控制在0.001 m左右，并逐步稀疏均匀过渡到其他边界，网格划分数量约1 200万，经检验数值计算的网格独立性满足要求。定义叶型屋面建筑群宽度184 m为B，计算域及网格划分如图3所示。

2.2　边界条件

入口边界条件：

空气相入口风速如式（6）所示：

v(\mathrm{z})=v_{\mathrm{0}}(z/z_{\mathrm{0}}{)}^{\alpha }（6）

其中\alpha取0.15；参考高度z_{\mathrm{0}}取10 m；参考高度风速v_{\mathrm{0}}取6 m/s^［14］。

雪相入口浓度在跃移层和悬移层由公式（7）定义^［11］：

（7）

h_{\mathrm{s}}=\frac{\mathrm{1.6}{u}_{\mathrm{*}}^{\mathrm{2}}}{\mathrm{2}g}（8）

式中：{\rho }_{\mathrm{a}}为空气密度；h_{\mathrm{s}}为跃移层高度；g为重力加速度；u_{\mathrm{p}}为雪颗粒的平均跃移速度，等于\mathrm{2.8}\times u_{\mathrm{*}t}；z为空间高度。

出口：远场出口边界，对应压力、速度梯度为0；

计算域顶部和侧壁：对称边界（symmetry）；

建筑物表面及地面：无滑移壁面（wall）；

积雪表面：忽略积雪表面对空气相的影响，对空气相和雪相为无滑移壁面（wall），但额外通过式（3）和式（4）计算其源项的设置方式来表示积雪表面对雪相的影响；

求解设置：压力-速度耦合算法选用SIMPLE算法，空间离散格式为Second-order upwind，流场各变量的收敛准则为残差达到10^-5。

2.3　结果处理

数值模拟结果通过定义不均匀积雪分布系数{\mu }_{\mathrm{r}}来讨论，其定义如式9所示。

{\mu }_{\mathrm{r}}=S_{\mathrm{h}}/S_{\mathrm{0}}（9）

式中：S_{\mathrm{h}}为风致雪漂移后的屋面积雪荷载；S_{\mathrm{0}}为不受风致影响下的初始屋面积雪荷载。

当不均匀积雪分布系数{\mu }_{\mathrm{r}}大于1时表示屋面此处积雪表现为沉积，小于1时对应侵蚀。

3 结果分析

针对8个屋面整体的风致积雪重分布进行了共12个风向角下（0°~360°，间隔30°）的数值风洞模拟，分别研究了不同风速、风向影响下的屋面积雪分布特性。

3.1　风速影响

风速对屋面积雪分布有直接影响，为研究其影响效应，分别考虑来流风速为1、3、6、9、12、15 m/s风速以及270°风向角作用下的屋面积雪分布进行研究分析，整体屋面的平均积雪系数如图4所示。

由图4可看出，屋面整体积雪分布系数随风速的增大而降低，同等情况下，无遮挡的房屋积雪分布系数随风速变化更明显。结合我国城市降雪期风速的概率分布在6 m/s以下的比例接近85%^［14］，因此后续论述重点针对6 m/s的来流风速进行进一步分析。

3.2　屋面整体分析

受篇幅限制，本文仅选取具有代表性的4个风向角下（60°、150°、210°和270°）和各风向角下最大最小积雪分布系数结果进行讨论，其分布云图如图5所示。

由图5可看出，风向角对屋面积雪分布情况有显著影响，对比来看迎风向的屋面边缘侵蚀比较严重，背风向屋面普遍沉积，且建筑迎风方向几何外形也对分布有影响，拱起区域侵蚀大于凹陷区域。从最大最小积雪分布系数云图可看出，侵蚀较多区域主要集中在屋面边缘及叶形屋面两端尖部区域，沉积较多区域则主要集中在叶面凹陷区域。

图6给出了每个屋面在不同风向角下积雪分布系数的均值曲线，可以看出：①屋面整体积雪分布系数通常小于1，说明屋面积雪以侵蚀为主，各屋面在风向角60°时积雪沉积量较低。②屋面1在0°~90°度的情况下受屋面2尾流和屋面3翘起头部的影响，导致临近区域侵蚀较大，其他情况下因为所处位置高不易受周围建筑影响，积雪分布基本稳定。③屋面2、3、4受迎风方向屋面1、5、6、7遮挡而导致在180°~270°积雪分布系数较大。通过计算比较不受遮挡的情况，受遮挡情况下积雪分布系数屋面2增加18.4%、屋面3增加16.3%、屋面4增加19.7%，基本控制在增加20%附近。④屋面5、6、7、8在90°~180°风向下积雪分布系数较大；此外屋面5因为屋面6、7、8侵蚀的积雪在屋面5后部分发生沉积，但前部分建筑间的空隙对风起到加速作用导致侵蚀严重，因此在210°和240°风向角下积雪分布较不均匀。

总体分析可看出，屋面积雪分布所受干扰主要分为三种：①屋面前方存在高屋面干扰时，屋面积雪侵蚀量降低、沉积量提高，主要是高屋面的遮挡导致屋面摩擦速度降低从而沉积以及屋面上部积雪浓度增加两方面原因；②屋面前方存在临近低屋面干扰时，屋面前部受低压区影响而导致沉积以及中部侵蚀；③当屋面受到前方建筑群空隙干扰时，狭管效应使得屋面摩擦速度增大从而提高侵蚀。

3.3　屋面局部分析

屋面局部几何外形对其积雪分布有直接的影响，取0°风向角下不受周围建筑影响的屋面3—6靠中部的三个位置截面分别进行讨论。屋面几何参数定义见图2叶形屋面外形几何尺寸定义示意图，横截面外形主要考虑D/B和H/B两个参数，分别代表截面的倾斜程度和凹凸程度，每个截面几何信息如表2所示。

根据截面形状，大致可将截面分为倾斜形（截面5-1、5-2、6-1等）和倾斜+凹凸形（截面3-1、3-2、4-2等）两类，其积雪分布系数如图7所示，可以看出：①屋面3具有较为明显的中部凸起特征，因此屋面积雪分布表现为典型的两侧较大而中部较小的分布特点，说明中部凸起区域风速得到加速效应，积雪侵蚀更为强烈；②屋面5、6的大部分区域整体倾斜较为一致，其积雪分布系数较为一致，可匹配建筑结构荷载规范建议取值；③屋面4一侧有一定的凸起，导致其积雪分布系数存在一定的降低；同时屋面4、6另一侧边缘存在局部凹陷，迎风状态下易形成较为强烈的侵蚀效果。

以3-1截面为例绘制截面积雪分布系数与曲面凹凸程度的关系曲线，如图8所示，可以看出对于凸起较为明显的曲面，积雪分布系数基本与凸起程度成反比关系。进一步绘制凹凸高度变化率与积雪分布系数变化率之间的关系，如图9所示，可以看出：二者呈现线性负相关，凹凸变化率越小说明越接近背风陡峭段，积雪分布系数提高得越快；凹凸变化率越大说明越接近迎风面陡峭区，积雪分布系数下降越快。因此，凹凸变化值应尽量保持在0.3以下，以避免积雪可能的过量堆积。

4 结 论

本文基于Euler-Euler方法与风雪单向耦合假定，借助数值风洞模拟技术，结合工程案例分析了叶形屋面风致雪漂移下的积雪分布，获得的主要结论如下：

（1） 叶形屋面几何造型较为复杂，其屋面风致积雪分布不仅要考虑屋面的倾斜，且应合理评估屋面的凹凸对积雪分布的影响。

（2） 叶形屋面积雪分布通常表现为侵蚀为主，且凸起部分和叶片两端的侵蚀更为明显，同时下凹部分可能形成局部积雪的严重堆积，应重视雪荷载验算过程。

（3） 对于一般倾斜叶形屋面，剖面长度越大，积雪分布越不均匀。当长度较小时，可根据倾斜程度取均值；而当长度较大时，由于前后部分分布差异过大，建议分段计算。

（4） 叶形多屋面的积雪分布相互影响较为明显，受低屋盖干扰影响区域主要集中于受遮挡作用屋面的迎风边缘区域，受高屋盖干扰下其积雪分布系数可增大20%左右。

（5） 对于即倾斜又凹凸叶形屋面，积雪分布系数随屋面拱起下降而提高，随着拱起的上升而减小，且下降变化与单位长度的拱起变化呈线性关系。实际设计中，为减少积雪堆积的风险，建议合理控制屋面的凹凸高度，工程设计中建议凹陷率不宜超过10%。

参考文献

[1] DEGAETANO A T，O'ROURKE M J.A Climatological measure of extreme snowdrift loading on building roofs［J］.Journal of Applied Meteorology，2004，43（1）：134-144.

[2] ZHOU X Y，HU J H，GU M.Wind tunnel test of snow loads on a stepped flat roof using different granular materials［J］.Natural Hazards，2014，74（3）：1629-1648.

[3] ZHOU X Y，LIU Z B，MA W Y，et al.Experimental investigation of snow drifting on flat roofs during snowfall：Impact of roof span and snowfall intensity［J］.Cold Regions Science and Technology，2021，190：103356.

[4] ZHANG Q，ZHANG Y，YIN Z，et al.Experimental study of interference effects of a high-rise building on the snow load on a low-rise building with a flat roof［J］.Applied Sciences，2021，11：11163.

[5] BEYERS J H M，SUNDSBO P A，HARMS T M.Numerical simulation of three-dimensional transient snow drifting around a cube［J］.Journal of Wind Engineering and Industrial Aerodynamics，2004，92：725-747.

[6] CHOI Young-Bae，KIM Rack-Woo，LEE In-bok.Numerical analysis of snow distribution on greenhouse roofs using CFD-DEM coupling method［J］.Biosystems Engineering，2024，237：196-213.

[7] SUN Xiaoying，ZOU Jialin，WU Yue，et al.Investigation of effects of repose angle and critical sliding angle of snow on curved roofs［J］.Journal of Wind Engineering and Industrial Aerodynamics，2022，227：105071.

[8] ZHANG Zhibo，MA Wenyong，LI Qiang，et al.Snow load shape coefficients and snow prevention method for stepped flat roofs［J］.Applied Sciences，2023，13：12109.

[9] 康路阳，周晅毅，顾明.考虑积雪休止角的屋面积雪漂移数值模拟方法［J］.同济大学学报（自然科学版），2016，44（1）：11-15.

[10] KANG Luyang，ZHOU Xuanyi，GU Ming.Numerical simulation mehod about snowdrift on roofs with consideration of repose angle of snow［J］.Journal of Tongji University （Natural Science），2016，44（1）：11-15.（in Chinese）

[11] 孙晓颖，洪财滨，武岳，建筑物周边风致雪漂移的数值模拟研究［J］.工程力学，2014，31（4）：141-146.

[12] SUN Xiaoying，HONG Caibin，WU Yue，et al.Numerical simulation of snow drifting around build model［J］.Engineering Mechanics，2014，31（4）：141-146.（in Chinese）

[13] 洪财滨.典型形式大跨度屋面风致雪漂移的数值模拟研究［D］.哈尔滨：哈尔滨工业大学，2012.

[14] HONG Caibin.Numerical simulation of snow drifting on typical form of long-span roofs［D］.Harbin：Harbin Institute of Technology，2012.（in Chinese）

[15] 李恺轩，马啸晨，李兴.漂移雪荷载的各国规范比较及减灾设计初探［J］.结构工程师，2015，31（5）：18-23.

[16] LI Kaixuan，MA Xiaochen，LI Xing.Drift snow loads in load codes of different countries and disaster reduction design［J］.Structural Engineers，2015，31（5）：18-23.（in Chinese）

[17] KINDＲJ.Mechanics of aeolian transport of snow and sand［J］.Journal of Wind Engineering and Industrial Aerodynamics，1990，36（Part 2）：855-866．

[18] ZHOU Xuanyi，ZHANG Yu，KANG Luyang，et al.CFD simulation of snow redistribution on gable roofs：Impact of roof slope［J］.Journal of Wind Engineering and Industrial Aerodynamics，2019，185：16-32.

基金资助

国家自然科学基金资助项目(51778365)