任超洋，尤本祥，宋志强，孙远辰，温阿强

(中国十七冶集团有限公司, 安徽 马鞍山 243000)

在工程实际中大跨度空间结构建筑得到了十分广泛的应用，其发展速度相对以前加快，如我国在2008年奥运会筹备期间已经建好了包括“鸟巢”“水立方”等举世闻名的大跨度空间结构建筑。大跨度结构是一种新颖的建筑结构形式[1]，通常跨度在30 m以上的建筑统称为大跨度建筑，其特点主要是质量轻、阻尼比小、内部可用空间较大、自振频率低且柔度较大，主要应用在展览馆、火车站、飞机场、大会堂等大型公共建筑上。目前对于大跨度屋盖结构的风荷载研究已取得大量数据，但还不能像高层结构和高耸结构那样在《规范》中给出一个具有代表性的风荷载分布形式。大跨度空间结构建筑的屋面平均风荷载的获取较为简单，研究的关键部分在于如何获取其脉动部分的风压。20世纪以来，多以现场实测和风洞试验获得数据，随着计算机技术的日渐成熟，利用数值方法获得数据的方法已经日趋流行。文献[2]通过现场实测方法对Rugby通讯站测量了4个不同高度的风荷载，并且在Cardington进行了观测，得出了不同高度的风速与湍流积分尺度变化的关系。文献[3]通过现场实测方法在1986年对悉尼的Belmore体育场进行了实测，证明了实测与风洞试验的风压系数平均值和极值较为接近，但是脉动风压系数较小。文献[4]通过风洞实验于1893年在澳大利亚利用喷气式飞机产生的气流测得了缩小后的房屋模型的风荷载分布情况。文献[5]在1955—1957年通过加长风洞来流试验段的长度成功模拟出了大气边界层。以上研究都是利用现场实测和风洞试验来解决建筑物的湍流来流问题，本文的中心思想是利用数值模拟方法对大跨度平屋盖、大跨度鞍形屋盖、大跨度球面屋盖以及大跨度悬挑屋盖进行数值模拟，主要目的是得出其平均风压系数和脉动风压系数。数值模拟技术的基本思想是：把原来在时间域以及空间域上连续的物理量的场用一系列有限个离散点上的变量值的集合来代替，并通过求解代数方程组获得场变量的近似值。国内诸多专家学者都利用此方法完成了一些计算工作，比如陈善群等[6]通过运用标准k-ε模型、重整化数群k-ε模型(RNGk-ε)和雷诺应力模型(RSM)对二维单山包和多山包绕流进行了数值模拟，确定了数值模拟对于山包数值模拟的可行性。卢春玲等[7]利用数值模拟中的大涡模拟方法对长沙机场扩建航站楼屋盖结构进行了计算，得出了屋盖表面的风压分布以及最不利风向角，并利用均方根误差公式得出脉动风荷载的研究方法。本文将在这一基础上对4种大跨度屋盖进行计算，本文所运用的方法是大涡模拟，这是一种能够在保证计算精度的情况下占用计算资源相对较少的方法。

1 研究方法与边界条件设定

本文采用大涡模拟方法对平屋盖、球面屋盖、鞍形屋盖和悬挑屋盖进行计算，并求得它们的平均风压系数和脉动风压系数。大涡模拟方法包括以下3个方面：

1.1 滤波器

1) 谱空间低通滤波器

(1)

在滤波空间的各个方向上用相同的滤波器，又称为各向同性低通滤波，其数学表达式为

Gl(k)=θ(kc-|k|)

(2)

式中θ表示台阶函数，当|x|<0时，θ(x)=0；当|x|>0时，θ(x)=1。用kc=π/l表示截断波数，其中l是相当的物理空间滤波单位长度。

2) 物理空间的盒式滤波器

(3)

物理空间的滤波器同样满足正则条件，此条件可以保证过滤器内各种物理量的守恒性，对于常数而言在此条件内是不变的。

(4)

式中：Ω是过滤的空间体积。1维盒式的滤波器可以表示为

(5)

式中：l是滤波器的单位长度，物理空间盒式滤波器的特点就是以空间长度为标准，当脉动长度小于滤波器的单位长度时将会被过滤掉。

1.2 亚格子尺度模型

此模型下的速度张量涡黏系数vt可以如下表示：

(6)

(7)

(8)

(9)

1.3 边界条件

为了更好地得到入口处的风速条件，并准确模拟入口处风的特征，在设置速度时采用已编辑好的UDF，使用UDF可以实现呈现指数形式的平均风剖面：

(10)

图1 风速剖面

入口处湍流强度边界条件除了速度外还有两个条件需要用UDF进行定义与描述， 这两个条件分别是湍动能k和湍流耗散率ε。湍动能k和湍流耗散率ε的数学表达式如下：

k(z)=1.2[I(z)*Vz]2

(11)

(12)

(13)

上表面、下表面和前后壁面条件设置为固壁无滑移边界条件(Wall)，左边界的边界条件设置为速度入口(Velocity-inlet)，右边界的边界条件设置为自由出流(Outflow)，方柱表面同样采用固壁无滑移边界条件(Wall)。

2 数值验证

本算例利用ICEM CFD的几何建模功能生成大跨平屋盖结构的三维模型，如图2所示。数值模型为一个平屋盖结构形式，平屋盖结构被置于一个矩形的流场中，流场四周封闭。对于所需计算的平屋盖结构来说，计算区域是一个足够大的空间，整个流场区域尺寸为80H×10H×15H。其中模型的长为80H、宽为10H、高为15H，模型的原点O定在方柱的左下侧顶点处。计算模型平屋盖结构的尺寸的长和宽为600 mm×600 mm，高定义为基本尺寸H=200 mm，因此在数值模拟中采用与风洞试验一致的参数。

图2 平屋盖计算域

2.1 平均风压系数验证

针对大跨度平屋盖选取一系列的观测点，大跨度平屋盖的上表面是受风力影响最重要的部分，将平屋盖面上沿着半个屋面的轮廓路径a→b→c→d→a按照每段之间取10个点来研究其平均风压系数，如图3所示。由于本算例需证明大涡模型(LES)的计算精确度，所以不仅提取了该模型下的数据点，同样也提取了利用k-ε模型计算的数据点用以对比。利用前面两种计算模型提取平均风压系数测点上的数据之后，与试验数据[8]进行对比，结果如图4所示。结果显示，大涡模拟(LES)较实验数据相接近，而k-ε模型的数据所展现的趋势较实验数据相差较大。相比较之下，大涡模拟(LES)的数据结果较好，平均风压系数的风洞试验结果和数值模拟结果如图5、6所示。

图3 观测点的选取

图4 数据对比

图5 风洞试验结果

图6 数值模拟结果

2.2 脉动风压系数验证

利用LES计算得出的脉动风压系数与试验结果的对比如图7所示。与平均风压系数的取法一样，沿着半个屋面的轮廓a→b→c→d→a提取脉动系数的值，作图与试验数据进行比较，如图8所示。

图7 LES法与风洞试验对比

图8 观测点数据对比

通过对结果的分析，可以得出如下说明：

1) 脉动风压系数表达式为

(14)

均方根误差表达式为

(15)

通过上式可以看出：脉动风压系数的表达式与均方根误差的表达式是一样的，均方根误差是观测值与真值偏差的平方和观测次数n比值的平方根，均方根误差可以很好地反映测量的精度。同理，脉动风压系数也是观测点风压和平均风压的差求均方根，所以可以得出一个结论：脉动风压系数与试验对比的准确度反映了所求风场的精度。

2) 利用LES计算得出的结果与风洞试验结果相比存在误差，但是最大误差相差在10%以内。根据平均风压系数Cp和脉动风压系数RMSE数据与风洞试验数据对比，发现两者的基本发展趋势是一样的，数值存在误差，但误差不大，因此利用LES方法对建筑物在大气边界层中湍流进行数值模拟是可行的。

3 大跨度屋盖结构的LES数值模拟

3.1 球面屋盖

杭州国际会议中心以举办大型国际会议和白金五星级酒店为标准进行功能设计，其采用钢结构建设，建筑整体高度为85 m，主体由13 m高椭圆形裙房和直径85 m球体组成。如图9所示，杭州国际会议中心是一个大跨度球面屋盖结构，球体部分可以简化为如图9右图所示的模型，其中OA=OC=50 m、DA=DB=40 m、OD=30 m。整个模型的计算域为1 220 m×920 m×340 m，将球体迎风面与左侧面交点投影在底面设为0点，如图10所示，从速度入口到建筑本体迎风面切平面的距离为400 m。

图9 杭州国际会议中心

图10 计算域示意图

3.2 球面屋盖的平均风压系数和脉动风压系数

得到计算结果后，在杭州国际会议中心模型表面取X=0 m和Y=0 m这两条线，如图11所示，并且在这两条线上各取22个点。由于国家规范并没有规定球面屋盖体型系数参考点的选定，本算例取球体顶点为参考点，球体顶点坐标为(110,110,85)，分别对这两条线上的点求量纲为1值，保留小数点后3位，取得点上的平均风压系数，如图12所示。根据脉动风压系数公式(14)(15)取得测点上的脉动风压系数，如图13所示。

图11 杭州国际会议中心风压系数选取线

图12 杭州国际会议中心平均风压系数

图13 杭州国际会议中心脉动风压系数

3.3 鞍形屋盖

大跨度鞍形屋盖不仅具有类似方柱绕流的流动特性，还具有特殊的负高斯曲率弧形屋面的形状特性。模型参考文献[9]，使用前处理软件ICEM对0°、45°和90°这3个风向角下的大跨度鞍形屋盖进行建模和网格划分，整个模型的计算域为9.35 m×4.25 m×1.2 m，从速度入口到鞍形屋盖的距离为3.4 m，其整体流域如图14所示。

图14 鞍形屋盖整体流域

3.4 3个风向角下大跨度鞍形屋盖的平均风压系数和脉动风压系数

对计算结果进行处理，在鞍形屋盖表面各取2条交于中点的线段，对于0°的鞍形屋盖取迎风处最高的顶点为参考点O点，参考点坐标为(3.4,2.125,0.38)；对于45°的鞍形屋盖取迎风面下侧最高点为参考点O1点，参考点坐标为(3.4,1.7,0.38)；对于90°的鞍形屋盖取迎风处的最低点为参考点O2点，参考点坐标为(3.4,2.125,0.1)，如图15所示。在每条线上取10个点并且求它们的无量纲值，保留小数点后3位，取得每个点上的平均风压系数，如图16所示。再根据脉动风压系数公式(14)(15)求得每个点上的脉动风压系数，如图17所示。

图15 鞍形屋盖参考点和数据提取线的选取

图16 3个风向角下的平均风压系数

图17 3个风向角下的脉动风压系数

3.5 悬挑屋盖

在目前诸多的大跨度空间结构中，悬挑结构造型独特且具有遮阳避雨的功能，属于典型的风敏感结构，而风荷载是悬挑结构的控制荷载之一。本文的大跨度悬挑屋盖形式与文献[10-12]中试验的几何形式保持一致，将看台简化成一个梯形，如图18(a)所示，整个模型采用缩尺比例，具体尺寸为：从地面到屋盖顶部H=0.25 m，从地面到看台的距离H1=0.2 m，从看台顶部到屋盖顶部的距离H2=0.05 m，看台顶部d=0.02 m，看台底部的宽度D=0.22 m。也就是说，简化的三角形看台的倾角度数α=45°，屋盖宽度D=0.25 m，屋盖厚度h=0.01 m，整个屋盖和看台的长度取相同值L=L1=1 m，屋盖倾角取0°。对于闭口式的悬挑屋盖，在看台和屋盖之间还有厚度为0.02 m的1个隔墙，2种模型剖面图见图18(b)和18(c)。为了使计算结果更准确以及避免四周边界对流场的影响，2个模型的计算域均取10 m×4.5 m×1.6 m。为了使来流充分发展，气流入口处距离模型迎风面1.25 m，来流风向沿x轴方向并且与结构侧面保持平行，见图19。目前对于悬挑屋盖网格的划分一般是以非结构性网格为主，也有混合型网格。通过前面的介绍已经知道了结构性网格的计算精度要高于其他形式的网格，所以本节悬挑模型同样采用结构性网格划分。本算例采用了“Y”型剖分使结构性网格更贴合悬挑的表面，网格划分完成之后悬挑结构表面的网格如图20所示。

图18 计算模型的剖面图

图19 2种悬挑模型及相应的计算域

图20 2种悬挑结构表面结构性网格划分

3.6 2种形式大跨度悬挑结构的平均风压系数和脉动风压系数

得到计算结果之后，先在悬挑结构表面选取1个参考点，然后再选取2条参考线，并在每条线上提取10个数据点，通过计算最后可得平均风压系数和脉动风压系数的分布。其中，参考点选取屋盖迎风面最左下角的点，坐标为(1.2,1.75,0.25)。取屋盖表面对称于X轴和Y轴的两条线Line1和Line2,如图21所示，并在这两条线上各取10个数据点，取其无量纲值，保留小数点后3位，根据公式求得每条线上的平均风压系数和脉动风压系数，如图22～25所示。

图21 悬挑结构参考点和数据提取线的选取

图22 2种悬挑结构Line1处的平均风压系数

图23 2种悬挑结构Line2处的平均风压系数

图24 2种悬挑结构Line1处的脉动风压系数

图25 2种悬挑结构Line2处的脉动风压系数

4 结论

通过对4种大跨度屋盖结构的大涡模拟计算，可以得出如下结论：

1) 首先通过对大跨度平屋盖的计算，并对比与设置条件相同的实验结果，发现通过大涡模拟可以得出较准确的结果。

2) 球面屋盖计算模型采用等比例简化过的杭州国际会议中心模型。通过观察可得，与X=110 m处的平均风压系数和脉动风压系数相比，Y=110 m处较为稳定，特别是Y=110 m处的脉动风压系数振荡幅度较大；鞍形屋盖的计算采用矢跨比为1/6的模型，对0°、45°和90°这3个不同风向角的计算，通过平均风压系数和脉动风压系数可知45°时的振荡最小，对鞍形屋盖的现场布置时需要考虑风向角的影响；对悬挑屋盖的计算采用最不利风向角位置(0°风向角)，设置悬挑形式为闭口和开口两种。相比之下悬挑顶部垂直于来流方向区域的平均风压系数和脉动风压系数振荡较大，出于安全考虑在建设时需要加固。