1. 引言

随着科学技术的发展，建筑工程技术、材料性能等方面有了长足的进步，各种复杂多样的结构形式层见叠出，建筑结构对于抗风性[1][2][3][4]的要求也越来越高。我国《建筑结构荷载规范》(GB 50009-2012)[5](简称荷载规范)中参考了以往的设计经验对一些简单的结构进行了风荷载计算方式的规定。但当复杂结构在面对风荷载作用时，其结构响应无法通过经验进行定量分析，需要对结构周围风场进行风荷载研究。

风洞试验[6][7][8]是目前结构风工程领域最主要的研究方法，该方法通过制作实体建筑物缩尺模型，利用动力装置使气流在管道内流动，以此来实现模拟真实建筑结构在风场中的风致振动现象。但由于风洞试验存在缩尺模型的物理量相似问题，无法模拟高湍流强度模型，同时风洞试验存在造价高、功耗大、模型利用率低、无法考虑风荷载与结构的耦合作用、易受洞壁干扰等问题。以计算流体动力学(Computational Fluid Dynamics, CFD)为基础的数值风洞逐渐成为了结构风工程领域新发展的计算方法。与传统的风洞试验相比，计算流体力学的模拟成本较低，高效灵活，并且可以进行结构全尺寸、高雷诺数计算。近年来，运用CFD技术对结构风工程领域相关问题进行的研究越来越多，2018年，吴颖等[9]基于CFD技术，以珠海华发艺术馆为工程背景，基于RANS湍流模型对复杂树状结构的各树冠部分表面风压进行研究，成功模拟出了结构表面的平均风压分布并给出了类似工程的最优湍流模型；2021年，郑得乾等[10]基于CFD技术，对萧山国际机场T4航站楼主楼屋盖表面风荷载进行了分析研究，得到了屋盖周围的流场，找到了屋盖较容易发生局部破坏的位置并给出了相关的抗风结构设计建议；2022年，刘行等[11]基于CFD技术，以厦门市复杂体型高层建筑为研究对象，得出了单体建筑在不同周边环境的条件下的风压分布规律，并用传统风洞试验验证了结论的正确性。然而，湍流控制方程是一个复杂的非线性偏微分方程，模拟计算时严重依赖计算机资源，以目前的算力，难以实现湍流流动的直接数值模拟。因此，在模拟过程中需对模型进行简化，使其适应特定的工程问题，常见的湍流计算模型有Spalart-Allmaras模型[12]、Standardk-ε模型[13]、RNGk-ε模型[14][15]、Realizablek-ε模型[16][17]、五方程雷诺应力(RSM)模型[18][19]等。此类湍流模型均为半经验模型，针对不同类型的工程问题求解精度差异很大。

本文以典型建筑国家大剧院[20]为实例(图1)，基于FLUENT软件平台，对国家大剧院周围风场的分布特性进行了数值模拟研究。应用雷诺平均法中定常的RNGk-ε湍流模型，探讨不同风向角时国家大剧院表面及周围风场，并研究周围建筑物，即人民大会堂对国家大剧院风场的影响。

Figure 1.National Grand Theatre

图1.国家大剧院

2. 数值模拟

2.1. 计算模型与网格划分

国家大剧院外部为钢结构壳体，外形呈半椭球形，东西方向长轴长度为216.25 m，南北方向短轴长度为146.25 m，建筑高度为46.125 m，其外形特征方程为：

${\left(\frac{x}{108.125}\right)}^{2.2}+{\left(\frac{y}{73.125}\right)}^{2.2}+{\left(\frac{z}{46.125}\right)}^{2.2}=1$(1)

与标准椭球面特征方程不同，其指数为2.2，模型相对复杂。综合运用GAMBIT和MATLAB软件，建立国家大剧院计算模型。将式(1)中x、y、z的值依次设为0，分别得到y-z平面、x-z平面、x-y平面椭圆方程，将生成的椭圆方程代入MATLAB中，将椭圆线分离成有限个坐标点，计算出各点的坐标值，将所有点的坐标值汇总后导入GAMBIT中，生成多个坐标点，将坐标点连接起来即生成对应的椭圆线。将z值分别设为5、10、15、20、25、30、35、40，采用相同方法可绘制出多个椭圆线，利用椭圆线生成椭球体表面，进而生成所需的椭球体。国家大剧院计算模型如图2所示。

Figure 2.Computation model

图2.计算模型

在CFD数值模拟中，计算域的横截面根据阻塞率确定，通常要求进口截面的阻塞率不大于5%。另外，计算域入口处到建筑物迎风面的距离一般不小于4h~5h(h为建筑物高度)，建筑物顶面和侧面到各自壁面边界的距离不小于4h；为保证湍流充分发展，计算域出口处到建筑物背风面的距离一般不小于9h~10h。针对上述要求，取计算域大小为300 m × 1000 m × 2400 m，计算模型中心距离入口800 m，距离出口1600 m。

计算域采用非结构四面体网格进行划分，计算模型附近区域网格划分密集，距模型较远区域网格划分稀疏。为减少网格数量，节省计算时间，在流场变化较快的边界层内，额外采用边界层网格进行加密。边界层网格中，第一层网格高度取0.1 m，递增系数取1.1，共划分10层。边界层外第一层网格尺度取2.5 m。经统计，计算域网格总数约70万，计算时间约6小时，计算精度和计算效率均能得到保证。网格划分如图3所示。

Figure 3.Mesh division

图3.网格划分

2.2. 边界条件与湍流模型

在模型计算域的边界上，共设定三类边界条件。在计算域入口边界采用速度入口边界条件，在计算域出口边界采用压力出口边界条件，在计算模型的表面及计算域四周，采用壁面边界条件。三类边界条件包含了计算域的所有边界，使计算域封闭。

在大气边界层内，风速在不同高度按规律变化，风速剖面是描述风速随高度变化的主要方法。目前对风速剖面的描述，比较理想的表示方法有对数律法和指数律法，两者相差不大。由于指数律法比对数律法计算方便，在结构风工程领域，普遍采用指数律法。荷载规范中规定，风速剖面的表达式为：

$v_z=v_{10}{\left(\frac{z}{10}\right)}^{\alpha }$(2)

其中，v₁₀为10 m高度风速平均值，取v₁₀= 8 m/s。α取值与地面粗糙度有关，依据国家大剧院的实际地理位置，粗糙度类别为C类地面粗糙度，取α= 0.22。

速度入口边界条件中的风速剖面、边界处的湍动能k和湍流耗散率ε均可通过用户自定义函数UDF在FLUENT中实现。

本文的数值计算针对不可压缩空气流动，为节省内存空间，采用分离隐式求解器，湍流模型采用定常RNGk-ε模型。速度入口边界和压力出口边界湍流强度取I= 21.5%。边界处的湍流通过给定湍动能k和湍流耗散率ε确定。其表达式分别为：

$k_{\left(z\right)}=\frac{3}{2}{I}^2{v}_{\left(z\right)}^2$(3)

${\epsilon }_{\left(z\right)}=C_{\mu }^{\frac{3}{4}}\frac{k_{\left(z\right)}^{\frac{3}{2}}}{l}$(4)

其中，C_u为湍流模型中的经验常数，近似取0.09；l为湍流长度尺度，取l= 0.07L，L为入口或出口处的水力直径。

2.3. 计算格式与计算内容

采用一阶迎风格式和二阶迎风格式共同计算。在开始阶段时采用一阶迎风格式计算得到一个相对粗糙的解，计算收敛后再用二阶迎风格式计算以提高精度。这样既避免了一阶迎风格式在复杂流场计算方面误差较大的问题，又避免了二阶迎风格式计算时间长、收敛性差的问题。

Figure 4.Wind direction angle

图4.风向角

本次模拟共计算7个风向角，即−90˚、−60˚、−30˚、0˚、30˚、60˚、90˚；考虑周围建筑物，即考虑人民大会堂的影响。风向角分布图如图4所示，其中A、B、C、D为国家大剧院表面上的点。

3. 结果分析

3.1. 风压系数

风压系数是指风在建筑物表面引起的实际压力或吸力与来流风压的比值，其表达式为：

$C_p=\frac{p-p_{\infty }}{\frac{1}{2}\rho {\bar{v}}^2}$(5)

其中，p为测点绝对静压值；p_∞为该测点远前方上游参考高度处绝对静压值；ρ为密度； $\bar{v}$为参考风速，取离地面10 m高度处风速。

图5为考虑周围建筑物的影响，不同风向角时大剧院顶点风压系数变化情况。

Figure 5.Wind pressure coefficient of structure vertex

图5.结构顶点风压系数

结果表明：在任意风向角下，结构顶点的风压系数均为负值；0˚风向角时风压系数绝对值最大，随着风向角变化风压系数绝对值逐渐减小；0˚风向角时，由于狭道效应的影响，大剧院顶部靠近大会堂一侧分离现象明显，使顶部负压增大；−90˚风向角时，由于大会堂的阻挡效应，在大剧院上风向，风速减小，顶部风压系数相应减小；90˚风向角时，由于大会堂前风场回流，同样使得大剧院顶部负压减小。

图6为考虑周围建筑物的影响，不同风向角时大剧院表面A、B、C、D点的风压系数变化情况。其中：图6(a)对应A点的风压系数，图6(b)对应B点的风压系数，图6(c)对应C点的风压系数，图6(d)对应D点的风压系数。

(a) (b)

(c) (d)

Figure 6.Wind pressure coefficient of structure surface

图6.结构表面风压系数

在不同风向角下，各点的方位在迎风面、侧风面、背风面之间转变，在迎风面时由于阻挡作用风压系数为较大正值，在侧风面时由于分离作用风压系数为较大负值，在背风面时风压系数较小，近似为0。周围建筑物对各点风压系数有较大影响：当风向角为−90˚时，周围建筑物位于大剧院上风向，由于建筑物的遮蔽效应，4点风压均为负值，且风压系数相对较小；当风向角为0˚时，A、C点风压均为负值，由于狭道效应C点风压系数绝对值较A点大；当风向角为90˚时，周围建筑物位于大剧院下风向，阻碍大剧院下游风场流动并产生回流现象，故C点出现正风压。

3.2. 速度矢量场与风压分布

对于0˚风向角，图7为结构体系的速度矢量图，图8为风压分布等值线图。

Figure 7.Velocity vector field 1

图7.速度矢量场1

Figure 8.Wind pressure distribution 1

图8.风压分布1

建筑物迎风面产生正压，侧风面和顶面为负压，背风面由于再分离作用产生旋涡。由于周围建筑物的影响，在靠近大会堂一侧的侧面风压和顶部风压分布不同，当风从前方流经大剧院和大会堂之间的空间时，由于狭道效应，流速增加，引起靠近大会堂一侧表面的流动分离现象增加，负压值增大，顶部负压区向人民大会堂一侧移动。

对于−90˚风向角，图9为结构体系的速度矢量图，图10为风压分布等值线图。

Figure 9.Velocity vector field 2

图9.速度矢量场2

Figure 10.Wind pressure distribution 2

图10.风压分布2

在该风向下，周围建筑物位于国家大剧院的上风向，阻挡了大剧院前方的风场流动。当风流经大会堂时，由于正前方两侧流动分离，在大会堂两侧产生负压；当风到达大会堂两侧尾部后，大会堂两侧墙壁突然消失，风向大会堂后方中央流动，并在大会堂后方形成旋涡。由于周围建筑物的遮蔽效应，大剧院前方风速较小，导致大剧院表面风压较小；大剧院位于周围建筑物后方旋涡区，受旋涡影响，在大剧院表面负压占主导作用，但负压值相对较小。整个体系对称分布，但大剧院表面风压分布并不对称，主要因为大剧院距离周围建筑物较近，处于周围建筑物的后方卡门涡街区域，故表面压力分布比较混乱。

对于90˚风向角，图11为结构体系的速度矢量图，图12为风压分布等值线图。

Figure 11.Velocity vector field 3

图11.速度矢量场3

Figure 12.Wind pressure distribution 3

图12.风压分布3

整个大剧院表面，除迎风面外，其他位置受周围建筑物影响均较大。由于周围建筑物的阻挡，大剧院后方部分流体反向回流，导致背风面风压值增大，且为正风压。

3.3. 对比分析

风洞试验在北京大学环境中心国家重点实验室的大气边界层风洞中进行[20]。国家大剧院的缩尺模型比例为1/200，由硬木制成，人民大会堂模型由轻型材料制作。将风洞试验结果与数值模拟结果进行对比分析，选取国家大剧院的顶点风压系数为例，如表1所示。

Table 1.Data comparison

表1.数据对比

风洞试验时，由于空间限制，风洞壁面会对试验结果有较大影响，−90˚风向角时壁面影响最大，结构顶点负风压明显增大。数值模型为理想刚性模型，风洞试验为近似刚性模型，故试验结果通常小于数值模拟结果。但总体而言，二者结果较接近，变化趋势相同，基本满足工程要求。

4. 结论

通过对国家大剧院的数值模拟研究，并将其与风洞试验结果进行对比分析，得出以下结论：

1) 通过对比风洞试验与FLUENT软件平台的模拟结果，结果显示数值模拟计算与风洞试验的结果大体上较为吻合，由于数值模型为理想刚性模型，风洞试验为近似刚性模型，试验结果通常小于数值模拟结果，但总体而言差值不大，且基本满足实际工程的精度要求，表明FLUENT软件平台可以较为准确地模拟建筑表面的风压分布，可以为此类结构的抗风设计提供依据和参考。

2) 大剧院的迎风面由于阻挡作用风压系数为较大正值，在侧风面时由于分离作用风压系数为较大负值，结构顶点风压系数始终为负。且在0˚风向角情况下，建筑物各点风压系数趋于极值，在实际设计过程中应对此工况予以考虑。

3) 大会堂位于大剧院上风向时，由于阻挡效应使风速减小；大会堂两侧产生负压，并在大会堂后方形成旋涡，受旋涡影响大剧院表面负压占主导作用，但负压值相对较小；大剧院距离大会堂较近，处于大会堂的后方卡门涡街区域，导致大剧院表面风压非对称分布。大会堂位于大剧院下风向时，除迎风面外，其他位置受大会堂影响较大；由于大会堂的阻挡作用，大剧院后方部分流体反向回流，导致背风面风压值增大，且为正风压。

4) 基于FLUENT软件平台数值模拟高度可视化的优点，可以直观地得到各区域的速度矢量场与风压分布图形，并对建筑风工程领域的研究分析具有积极的意义。

基金项目

陕西省自然科学基础研究计划(2023-JC-YB-435)。

NOTES

^*第一作者。

参考文献