平屋面摩擦速度的数值模拟研究

辛林桂，周晅毅，顾 明

(同济大学土木工程防灾国家重点实验室，上海 200092)

平屋面是一种典型的屋盖形式，其跨度对建筑周围的流场有着重要的影响[1-5]。风经过平屋面时会出现复杂的绕流、再附现象，在气流作用下雪颗粒将发生复杂的飘移运动，从而造成屋面和建筑周围积雪的不均匀分布[6-9]。大跨度的平屋面结构对于不均匀分布的雪荷载很敏感，雪荷载的不均匀分布形式往往是结构的最不利荷载分布之一。因此，准确预测不同跨度屋面上的积雪分布和传输特征有利于保障结构的安全。

目前对平屋面积雪分布和传输的研究方法主要有现场观测[10-11]、风洞或水槽试验[12-13]和数值模拟[14-15]三种。现场观测记录了平屋面在自然环境中的积雪分布和传输规律，反映了最为真实积雪飘移的运动。THIIS 和GJESSING[10]在挪威对平屋面进行了屋面积雪的实测研究，结果表明屋面尺寸对屋面积雪的飘移和分布有显著的影响；ZHANG 等[11]采用现场观测研究了不同比例的平屋面和其他屋面模型的缩尺比对屋面积雪分布的影响，观测结果表明几何缩尺比对屋面积雪分布影响甚微。与现场测量相比，风洞或水槽试验可以提供一个稳定可控的试验平台来研究屋面的积雪飘移。O’ROURKE 等[12]在水槽中利用核桃壳对平屋面积雪传输率进行了研究，结果表明，水槽模型与全比例模型的输运率吻合较好；QIANG 等[13]在低温风洞中使用人工雪颗粒对积雪在平屋顶上的飘移进行了模拟，研究了降雪期间与无降雪期间积雪飘移的异同。不同于现场实测和风洞或水槽试验，对屋面积雪的数值模拟更偏向于机理研究。LIU 等[14]利用拉格朗日方法建立了一个数值模型研究了平屋面的积雪飘移，探讨了击溅夹带对平屋顶积雪输运的作用及其贡献。QIANG 等[15]介绍了一种积雪传输模型，用于计算平屋面的风致积雪传输量，在该模型中考虑了降雪量和阈值摩擦速度等变量的差异，模拟结果与前人的实测结果较为一致。上述的研究采用不同方法分析了平屋面积雪的分布特征，但忽略了屋面跨度可能对积雪传输造成的影响；同时，O’ROURKE 等[16]通过研究发现饱和积雪传输率应按传输距离进行折减。由于目前对平屋面的积雪研究中，鲜见有学者考虑屋面跨度对屋面雪荷载的影响，因此结合屋面积雪传输率与传输距离分析跨度对平屋面雪荷载的影响是非常必要的。

积雪野外实测的研究结果表明，传输率与摩擦速度之间存在指数关系[17-18]。并且一些研究结果表明，建筑屋面的摩擦速度是决定其风致迁移雪荷载的重要因素[19]，如：康路阳等[20]采用分阶段的欧拉-欧拉方法模拟了有休止角的积雪表面的摩擦速度，并结合摩擦速度预测了平屋面的积雪分布；CAO 等[21]通过风洞试验和数值模拟研究了低坡度屋面的积雪质量传输率与摩擦速度之间的关系。因此，准确获得平屋面的摩擦是研究屋面积雪传输的前提条件。

本文首先采用CFD 方法模拟了平屋面的摩擦速度，并且为了研究不同湍流模型的模拟精度，利用Irwin 探头风洞试验与数值模拟结果进行了对比。接着，根据CFD 模拟的结果分析了不同跨度平屋面流场和摩擦速度的分布特征。最后，基于数值模拟得到的摩擦速度，计算了平屋面积雪传输率，并讨论了不同来流风速下屋面跨度对积雪传输率和侵蚀通量的影响。

1 数值模型

1.1 研究对象

为了获得不同跨度平屋面的摩擦速度，本文采用RANS 方法对二维平屋面的周围的流场进行模拟。研究对象的跨度L分别为30 m、60 m、90 m和120 m，高度H为15 m，屋面跨度与高度之比L/H分别为2、4、6 和8，如图1 所示。

图1 研究对象Fig.1 Research object

1.2 计算域及网格方案

CFD 数值模拟采用ANSYS FLUENT 软件进行，数值模型的几何缩尺比为1∶100，图2 为计算域和网格划分示意图(未按比例画出)，计算域的长度为16L，高度为20H。建筑迎风面到入口(velocity inlet)边界的距离为5L，出口(outflow)边界距离建筑背风面为10L， 如图2(a)所示。以模型L/H=4 为例，高分辨率的结构化网格如图2(b)所示，最小网格尺寸约为0.03H，网格的横向和纵向的增长率均为1.05，总网格数量约为5.5 万。

图2 计算域和网格划分示意图Fig.2 Schematics of computational domain and computational grid

1.3 边界条件和参数设置

计算域入口边界采用速度入口，参考《建筑结构荷载规范》对B 类地貌的规定，并按照几何缩尺比(1∶100)缩放后在数值模拟中进行应用。入口平均风剖面U(y)如下所示[22]：

式中：风剖面指数α 取0.15；为了探究风速对屋面摩擦速度和积雪传输率的影响，本文取U0.1(距地面0.1 m 高度处风速)的风速大小分别为6 m/s、8 m/s 和10 m/s。

湍流强度剖面I(y)为[22]：

根据风速剖面U(y)和湍流强度剖面I(y)，湍动能剖面k(y)和湍流耗散率ε(y)剖面为[23]：

式中，Cµ为0.09。

下边界和建筑表面的壁面边界条件采用无滑移标准壁面函数。上边界设置为对称边界，出口边界设置为0 梯度的出流边界。采用稳态RANS方法，对Standardk-ε、 RNGk-ε 和Realizablek-ε三种湍流模型进行了比较分析。所有的传输方程均采用二阶迎风离散格式，压力速度耦合方法为SIMPLE 算法，收敛准则为10-6。

1.4 网格无关性分析

采用三种尺寸的网格进行网格敏感性分析，分别为细网格、中网格和粗网格。细网格的最小网格尺寸为0.01H，屋面y+的范围为25～50；中网格的最小网格尺寸为0.03H，屋面y+的范围为31～70；粗网格的最小网格尺寸为0.05H，屋面y+的范围为38～95。图3 为模型尺寸L/H=4 时，采用Realizablek-ε 模型和三种网格划分方案模拟的屋面摩擦速度对比。从图3 可以看出，粗网格与细网格的差异较大，而中网格与细网格的差异较小，因此本文在下文中选择中网格的模拟结果进行分析。

图3 网格无关性分析Fig.3 Grid-sensitivity analysis

2 湍流模型对模拟结果的影响

为了验证Standardk-ε、 RNGk-ε 和Realizablek-ε 三种湍流模型的模拟精度，并对比不同湍流模型对模拟结果的影响：首先，对Irwin 探头进行了标定；接着，在不同来流风速下采用标定后Irwin探头测量了平屋面的摩擦速度；最后，将不同湍流模型的数值模拟和风洞试验进行了对比。风洞试验的风压测量采用美国Scanvival 公司生产的ZOC33 电子压力扫描模块，压力扫描阀的采样频率为312.5 Hz，采样时间为30 s。

2.1 Irwin 探头的标定

本文采用的Irwin 探头为无方向探头[24]，尺寸如图4 所示，测压孔K1到K2的垂直距离h为5 mm。

图4 Irwin 探头尺寸 /mmFig.4 Dimensions of the Irwin sensor

FERREIRA 等[25]和FARIA 等[26]认为壁面切应力与Irwin 探头测得压差存在指数关系，即：

式中：a、b为Irwin 探头测量壁面切应力时的标定系数；ΔP=P0-Ph，为Irwin 探头测得的压差。由摩擦速度u*=， ρa为空气密度，可得u*与ΔP的关系为：

式中，m、n为Irwin 探头测量摩擦速度时的标定系数。标定试验在同济大学TJ-2 大气边界层风洞的空风场中，根据WANG 等[27]的方法对Irwin 探头测量摩擦速度进行了标定，如图5(a)所示。假设标定试验的近壁面风速符合对数分布：

图5 Irwin 探头测量摩擦速度的标定示意图Fig.5 Calibration diagram of friction velocity measured by the Irwin sensor

式中：U(z)为z/m 高度处的风速，标定试验中不同风速下的风剖面U(z)分布如图5(b)所示；Von Karman 常数κ为0.42；z0为标定平台的气动粗糙度高度。通过式(7)的关系对风剖面进行拟合后，可得到z0为1.173×10-5m，五种不同风速下的u*分别为0.17 m/s、0.23 m/s、0.32 m/s、0.41 m/s 和0.50 m/s。于是再根据式(6)对高度h=5 mm 的Irwin 探头测得的 ΔP和u*进行标定，通过重复性和方向性检验后，可得到h=5 mm 的Irwin 探头在测量摩擦速度时的标定系数m=0.183、n=0.338。

2.2 Irwin 探头风洞试验

采用Irwin 探头在同济大学TJ-2 大气边界层风洞对平屋面的摩擦速度进行测量，风洞的试验段长15 m、宽3.0 m、高2.5 m，如图6 所示。屋面模型的几何缩尺比为1∶100，跨度为0.6 m，高度为0.15 m，阻塞比小于5%。在进行风洞试验时仅模拟风场的二维效应，模型的两侧增加了侧挡板，并且为了减小挡板对风场的影响，挡板端部进行了削尖处理。Irwin 探头在风洞试验模型上的布置如图7 所示。考虑到Irwin 探头相互遮挡可能会影响测量精度，故在屋面上对Irwin 探头进行了交错布置，两排Irwin 探头的间距为50 mm，Irwin探头距离迎风侧和背风侧屋檐的距离均为10 mm，其他的Irwin 探头在屋面上均匀排布。

图6 风洞试验试验段示意图 /mFig.6 Schematic diagram of wind tunnel test section

图7 风洞试验模型测点布置 /mmFig.7 Layout of measuring points of wind tunnel test mode

在测量平屋面的摩擦速度时，定义屋面高度H为参考高度，平屋面Irwin 探头风洞试验采用了三种风速：参考风速U(H)分别为6 m/s、8 m/s 和10 m/s，参考高度处的湍流强度I(H)约为16%，风速和湍流强度剖面如图8 所示。参考了康路阳等[20]对风洞试验风速的拟合方法，风速剖面可由式(7)进行拟合，气动粗糙度高度z0为5.362×10-5m，湍动能由式(2)进行计算，湍流耗散率的定义为：

图8 风洞试验风速和湍流强度剖面Fig.8 Wind speed and turbulence intensity profile of the wind tunnel test

式中，Cμ为经验常数，取0.09。

2.3 不同湍流模型的对比

为了对比Standardk-ε、 RNGk-ε、Realizablek-ε 三种湍流模型对数值模拟结果的影响，将CFD 模拟获得的屋面摩擦速度u*与风洞试验结果进行了对比。以参考风速U(H)为8 m/s 为例，数值模拟与风洞试验的对比如图9 所示。

图9 数值模拟与风洞试验的摩擦速度对比(U(H)=8 m/s)Fig.9 Comparison of friction velocity between numerical simulation and wind tunnel test (U(H)=8 m/s)

在三种风速下，采用FAC2(the factor of 2)、线性相关系数R和平均百分比误差APD 对数值模拟结果进行评价。FAC2、R和APD 的定义如下所示[28 - 30]：

式中：Pi和Mi(i=1, 2,…,N;N=48)分别为CFD 模拟的模拟值和风洞试验测量值； σP和 σM分别为CFD 模拟和风洞试验结果的标准差；上划线代表平均值，R理想值为1.0，APD 理想值为0%。

根据上述三种评价指标，不同湍流模型的数值模拟结果评价如表1 所示。在FAC2 评价中，Realizablek-ε 模型的摩擦速度的评价结果与理想值最为接近，这表明Realizablek-ε 模型的整体表现最好，很少发生明显过高或过低的预测。对比线性相关系数R可得，Realizablek-ε 模型模拟的摩擦速度与风洞试验测量结果相比相关性最好。在平均百分比误差APD 评价中，Realizablek-ε 模型的摩擦速度评价结果与理想值0%最接近。综上所述，Realizablek-ε 模型对摩擦速度的预测表现最好。因此，在后文的研究中选择Realizablek-ε模型的数值模拟结果进行分析。

表1 三种风速下不同湍流模型模拟结果对比Table 1 Comparisons of simulation results of different turbulence models under three wind speeds

3 结果讨论

3.1 不同跨度平屋面的流场和摩擦速度分布

以U10=8 m/s 为例，不同跨度平屋面的风速分布如图10 所示。由图10 可得，来流在屋盖迎风侧的前缘发生分离，迎风前缘最大风速约为来流风速U10(8 m/s)的1.5 倍。来流在迎风前缘发生强烈的分离，形成旋涡后发生再附，并在屋面后方形成一个回流区。屋面跨度L=2H时(图10(a))，分离产生的旋涡几乎覆盖了整个屋面；随着屋面跨度的增大(图10(b)～图10(d))，旋涡的长度也在不断增大，但旋涡的长度与屋面跨度的比值在不断的缩小。同时，旋涡后方近屋面的风速也随着跨度的增大不断增大。

图10 不同跨度平屋面周围的风速分布Fig.10 Wind speed distribution around flat roofs with different spans

摩擦速度可以确定雪颗粒运动状态，在数值模拟中屋面摩擦速度u*可由屋面剪切应力τw计算获得，即u*=。受近屋面风速的影响，不同跨度屋面的摩擦速度分布如图11 所示。由于来流在屋盖前缘分离的影响，摩擦速度在屋盖前缘处较大，并且不同跨度屋面的摩擦速度在屋盖前缘处大小分布基本相等；而后经过分离产生的旋涡覆盖区域，摩擦速度先增大后减小；最后旋涡覆盖区域后方风速平稳发展，摩擦速度也逐渐增大。参考孙晓颖等[31]对雪颗粒阈值摩擦速度的取值，本文取雪颗粒的阈值摩擦速度u*t为0.15 m/s。当屋面跨度与高度之比L/H=2 时，屋面摩擦在背风侧屋檐小于阈值，屋面其他位置均大于阈值，屋面整体呈现侵蚀状态。随着屋面跨度的增大，屋面摩擦速度的最小值所在的位置逐渐向迎风侧屋檐靠近，即屋面上产生积雪沉积的位置也逐渐靠近迎风侧屋檐。

图11 不同跨度屋面的摩擦速度分布Fig.11 Friction velocity distribution of roofs with different spans

图12 不同跨度的平屋面表面的平均摩擦速度Fig.12 Mean friction velocity of roofs with different spans

3.2 屋面跨度对积雪传输和侵蚀的影响

本节采用CFD 模拟获得的平屋面摩擦速度的分布讨论屋面跨度对屋面积雪传输的影响。屋面积雪传输率Qu*与数值模拟获得的平均摩擦速度的关系可由OWEN[17]和KIND[18]的研究结果获得：

根据屋面积雪传输率Qr，屋面积雪沿跨度的平均侵蚀通量qr可以表示为Qr与L比值：

图13 为不同风速下跨度对平屋面的积雪传输率和侵蚀通量的影响。由于来流风速为屋面积雪的传输和侵蚀提供了动力，对比图13(a)、图13(b)可得，屋面跨度相同时，风速越大屋面积雪的侵蚀就越剧烈，屋面积雪传输率和侵蚀通量也就越大。图13(a)显示屋面的积雪传输率随着屋面跨度的增大而增大。同时，从图13(b)可见，屋面积雪的沿跨度平均侵蚀通量随着屋面跨度的增大而减小。

图13 不同跨度平屋面的积雪传输率和侵蚀通量Fig.13 Snow transport rate and erosion flux of flat roofs with different spans

3.3 屋面跨度对雪荷载暴露系数的影响

为了进一步分析屋面跨度对屋面雪荷载的影响，本文参考ASCE 荷载规范[34]对平屋面雪荷载Pf的定义：

式中：Ce为暴露系数；Ct为热力系数；Is为重要系数；Pg为地面雪荷载，kg/m2。将Ct与Is取1.0 后，Ce可以表示为：

式中：Gg/(kg/(m2·s))为地面积雪的沉积通量，即单位时间内在单位面积上的降雪质量，新降积雪密度为100 kg/m3时，这里将Gg取为1.5×10-4kg/(m2·s)，对应的气象条件为暴雪天气[35]；T/s 为屋面风吹雪的持续时间；qr/(kg/(m2·s))为屋面积雪的侵蚀通量。图14 为不同风速下屋面跨度对雪荷载暴露系数的影响。由图可得，屋面跨度相同时，雪荷载暴露系数随着来流风速的增大而减小；并且在同样的气候条件(风速)下，跨度越大的屋盖结构，来流风速对其屋面雪荷载的折减就越小。

4 结论

本文采用CFD 方法模拟了平屋面的摩擦速度，并结合其分布特征分析了屋面积雪传输率、侵蚀通量和暴露系数随跨度的变化趋势，主要结论如下：

(1) 本文根据风洞试验对Irwin 探头测量摩擦速度进行了标定，结合Irwin 探头风洞试验对比不同湍流模型的CFD 模拟结果发现，Realizablek-ε模型的模拟结果与风洞试验最为接近。

(2) 对比不同跨度平屋面的CFD 模拟结果可得，随着屋面跨度的增大，屋面摩擦速度的最小值所在的位置和积雪沉积的区域均逐渐向迎风侧屋檐靠近。不同跨度平屋面的平均摩擦速度随着风速的增大而增大，并且屋面跨度与高度之比L/H=4 时，屋面的平均摩擦速度在本文选取的4 种跨度的屋面中最小。

(3) 屋面积雪的平均侵蚀通量随着跨度的增大而减小。并且在同样的气候条件(风速)下，跨度越大的屋盖结构，来流风速对其屋面雪荷载的折减就越小。