阻塞率对表面风压系数影响的数值模拟

基金项目：上海市博士后基金科研资助计划项目（07R214162）

作者简介：方平治（1974），男，河南平舆人，副研究员，工学博士，Email：freedomfpz@163.com。

摘要：以同济大学TJ2风洞模拟的缩尺比为1∶50的B类风场为目标，给出了基于标准kε湍流模型的数值风场。作为数值风场的一个应用，对TTU模型的定常绕流进行模拟，并与风洞试验以及现场实测结果进行比较，重点探讨了阻塞率对表面风压系数的影响。结果表明：对于TTU模型的屋顶和背风面，阻塞率增大可导致风压系数绝对值有增大趋势；对于TTU模型的迎风面，由计算域侧面引起的阻塞率增大，导致风压系数有减小趋势，而由计算域顶面引起的阻塞率增大，则导致风压系数有增大趋势。

关键词：TTU模型；表面风压系数；壁面函数；阻塞率；数值模拟；kε湍流模型

中图分类号：TU312 文献标志码：A

0 引 言

风洞试验是风工程研究的重要手段，然而存在周期长、耗费大、不能够提供足够多的流场信息等缺点，因此基于计算流体动力学的数值风洞得到日益广泛的应用。尽管计算机技术迅猛发展，但计算量仍是数值风洞在工程应用中需要考虑的一个重要问题，其中计算域是影响因素之一。数值计算中通常对计算域进行优化，以得到和计算域无关的数值计算结果，同时计算域尽可能小。近年来，在欧洲和日本分别有研究小组对计算流体动力学在工程中的应用进行研究，给出了一些建议，并对计算域问题进行了探讨[12]。

为了有效评价计算域，必须对风场进行高精度模拟。在计算域内再现大气边界层是数值风洞的一个基本要求。数值方法的误差不仅来源于湍流模型以及钝体绕流的复杂性，数值模拟大气边界层的不足也是一个重要因素。围绕数值模拟大气边界层问题，许多学者进行了研究。基于标准湍动能湍动能耗散率（kε）湍流模型，Richards等[3]研究了大气边界层的数值模拟方法，给出了和湍流模型相容的来流边界条件，包括平均风速、湍动能和湍动能耗散率。基于文献[3]中的研究成果，Yang等[4]重新推导了来流边界条件，给出了和现场实测以及风洞试验结果相吻合的湍动能来流边界条件。Blocken等[5]对数值模拟大气边界层存在的壁面函数问题进行了阐述。本文中笔者对壁面函数问题进行了研究[67]，从非均匀粗糙度入手，在标准壁面函数[8]基础上，增加一个附加项，克服了壁面函数问题对模拟大气边界层的不利影响；以同济大学TJ2风洞模拟的缩尺比为1∶300的B类和D类风场为来流边界条件，建立了基于kε系列湍流模型的数值风场，并对附加项的有效性进行了验证。这些数值风场满足现行《建筑结构荷载规范》（GB 50009—2001）[9]中关于平均风速和湍流强度的规定。

以同济大学TJ2风洞模拟的缩尺比为1∶50的B类风场为目标，本文中首先建立了基于标准kε湍流模型的数值风场。作为数值风场的一个应用，对TTU模型的定常绕流进行了计算，并与风洞试验和现场实测结果进行了比较；在此基础上，重点探讨了阻塞率对表面风压系数的影响。此外，还对出流面和入流面位置对表面风压系数的影响进行了计算。

1 基于标准kε湍流模型数值风场

1.1 标准壁面函数的附加项

对于大气边界层，其粗糙度是由大小不等、分布不均匀的粗糙元形成的。在数值模拟大气边界层研究中，为了考虑这种非均匀性，在标准壁面函数基础上，给出了一个附加项，即

式中：U为壁面附近的平均风速；u*为摩擦速度；κ为冯·卡门常数，κ=0.42；z为到壁面的垂直距离；Ks为物理粗糙高度；C2为常数；δB为附加项，δB=ln（α）/κ，α为参数。

参数α一般大于1，δB为正值，因此，具有附加项的标准壁面函数可以提供更小的平均风速，这对于具有较大粗糙度的地貌类型和小缩尺比的风场尤其重要。当α=1，δB=0时，式（1）即为标准壁面函数。附加项克服了单独依靠由近壁面网格提供粗糙度的限制，即克服了壁面函数问题。

1.2 数值风场

采用标准kε湍流模型，以同济大学TJ2风洞模拟的缩尺比为1∶50的B类风场为目标，本文中首先给出数值风场。数值风场参考高度HG=1.0 m处的参考风速UG=6.0 m·s-1。风场的粗糙长度与对应的湍流模型常数以及边界条件见表1，2，其中P为风压。墙边界条件采用具有附加项的标准壁面函数。附加项由FLUENT 6.3的用户自定义函数（User Defined Function，UDF）编程实现，并通过用户自定义壁面函数（Userdefined Wall Functions，UDWF）代替标准壁面函数。附加项δB=5.5，给出α=10，因此，z′0=z0/α=0.000 158 m，其中，z0为空气动力学粗糙长度。由于利用了附加项，满足FLUENT 6.3关于Ks方案中，zP=0.004 m。采用基于GreenGauss Node Based的求解器；速度压力耦合方式为SIMPLEC；对流项的离散格式为QUICK；压力插值方式为PRESTO；其他为FLUENT的默认设置。由于MMK kε湍流模型和标准kε湍流模型的差别仅在于湍流粘性系数，上述湍流模型常数和边界条件同样适用于MMK kε湍流模型。

2 理想计算域条件下TTU模型定常绕流的数值模拟2.1 TTU低矮建筑标准模型简介

TTU模型是Texas Tech University提出的一种标准模型[10]，该模型处于B类地貌类型中，实际尺寸为9.1 m×13.7 m×4.0 m。TTU研究小组在该模型表面上布置了许多压力测点，通过现场实测得到大量压力数据，这些数据被很多研究小组用来检验风洞试验结果和数值模拟结果的有效性。本文中对TTU模型的定常绕流进行计算，选择的工况为90°和60°两个风向角，并将数值计算结果与现场实测以及风洞试验[11]结果进行比较，以验证数值风场以及数值模拟结果的可靠性。TTU模型、风向角β的定义和模型壁面测点位置如图1所示。

2.2 数值模拟结果

2.2.1 理想计算域

风洞试验中，通常要求阻塞率足够小，以减小风洞的侧面和顶面对建筑物绕流产生的不利影响[1214]，上述原则在数值计算中同样适用。本文中以90°风向角工况为例，从无限大计算域出发，通过数值试验寻求该理想计算域。经过数值试验，理想计算域在流向x、横向y和竖向z三个方向的尺寸分别为：-6H≤x≤20H，-10H≤y≤10H，0≤z≤7H，H为参考高度，其中，H=0.08 m为缩尺模型的高度，相应的阻塞率为2.5%。图2（a）中给出了该理想计算域和坐标系，其中坐标原点位于TTU建筑模型底面中心。采用非均匀结构化网格，建筑物壁面附近的最小网格尺寸控制在0.075H（注意和计算域紧贴地面第1层网格的垂直尺寸0.008 m，即0.1H的区别）。图2（b）中给出了建筑物壁面上的均匀网格，流向x、横向y和竖向z三个方向上的面网格单元数分别为24，34，10。整个计算域的网格数为70×104。

以入流面位置的优化为例，图3中给出了y=0面上x=-6H（理想计算域的入流面），-5.5H，-5H，-4H，-3H处的平均风速、湍动能剖面和来流边界条件的比较。由图3可见：在x=-3H处，无论是平均风剖面，还是湍动

能剖面，均受到建筑物的干扰，而在x=-5H处，平均风剖面和湍动能剖面基本不受建筑物的影响。因此，可将入流面放置在x≤-5H位置处，本文中取x=-6H。计算域出流面位置、顶面位置以及侧面位置的优化过程均采取相同的思路。

2.2.2 90°风向角的计算结果及其分析

利用上述理想的计算域和标准kε湍流模型，对TTU模型的定常绕流进行计算。在迎风面屋檐附近，由于绕流复杂，图2中的均匀网格已出现不足。图4（a）中给出了建筑物壁面上的非均匀网格，壁面附近最小网格尺寸为0.062 5H，在流向x、横向y和竖向z三个方向上的网格单元数分别为26，34，13，对应的整个计算域的网格数为75×104。 图4（b）中给出了模型横截面上11个测点（参考图1）平均风压系数的数值计算、风洞试验以及现场实测结果的比较。测点i的平均风压系数CPi的定义为

式中：Pi为测点i的平均压力；PH为参考高度H处的静压；UH为参考高度H处的来流平均风速；ρ为空气密度，ρ=1.225 kg·m-3。

由图4可见：①风洞试验结果再现了现场实测结果，说明风洞试验结果是可靠的；②在迎风面屋檐附近，数值计算结果与风洞试验结果的差异较大，而在其他测点处，两者基本吻合。由于标准kε湍流模型自身的不足，导致湍动能在迎风面屋檐处堆积，从而导致风压系数过高估计，本文中采用MMK kε湍流模型进行计算，计算结果如图4（b）所示。由此可见，MMK kε湍流模型能够给出和风洞试验较为接近的计算结果。

2.2.3 60°风向角的计算结果及其分析

在风向角为60°时，屋顶出现三角翼涡，模型周围的绕流更为复杂。图5（a）中给出了建筑物壁面上的非均匀网格，壁面附近的最小网格尺寸为0.062 5H，对应的整个计算域的网格数为94×104。由于网格和风向出现倾斜交叉，数值计算误差（数值粘性）会增大，因此，需要采用更高阶的对流项离散格式。图5（b）中给出了MMK kε湍流模型对流项离散格式分别采用QUICK（MMK/Secondorder）和MUSCL（MMK/Thirdorder）的计算结果（共10个测点）。由图5可见：由于采用了三阶的对

流项离散格式，迎风面屋顶的数值计算精度得到明显提高。当风向角为240°时，角点（参考图1）的风洞试验结果为CP=-0.83；现场实测结果为CP=-0.88；而数值计算结果为CP=-0.85（根据对称性，由本文风向角为60°时的计算结果得出）。

3 计算域

在数值计算中，合适的计算域可以保证在满足计算精度要求的同时大幅度减少网格数量。本文中的理想计算域较大，在数值计算中不实用。下面将进一步优化计算域，寻求在工程应用中合适的计算域。以90°风向角为例，在理想计算域基础上，即-6H≤x≤20H，-10H≤y≤10H，0≤z≤7H，依次将出流面的位置变为10H、侧面的位置变为6H、入流面的位置变为5H以及顶面的位置变为4H，与之对应的计算结果分别以MMK_D10H，MMK_D10H_S6H，MMK_D10H_S6H_U5H，MMK_D10H_S6H_U5H_T4H表示，如图6所示。最终的计算域为-5H≤x≤10H，-6H≤y≤6H，0≤z≤4H，对应的整个计算域的网格数由75×104降到36×104，阻塞率由2.5%升到7%。由图6可见：①尽管计算域和网格数量大幅度减小，上述计算结果还是可以接受的；②计算域的不足会对计算结果造成一定影响。具体而言，对于TTU模型的屋顶和背风面，阻塞率增大可导致风压系数绝对值有增大趋势；对于TTU模型的迎风面，由计算域侧面引起的阻塞率增大导致风压系数有减小趋势，而由计

算域顶面引起的阻塞率增大导致风压系数有增大趋势。另外，入流面和出流面位置对计算结果也会产生影响。对于TTU模型屋顶和背风面，入流面和出流面距建筑物越近，风压系数绝对值越大；对于TTU模型迎风面，出流面位置距建筑物越近，风压系数越小，而入流面位置距建筑物越近，则风压系数越大。

4 结 语

以同济大学TJ2风洞模拟的缩尺比为1∶50的B类风场为来流边界条件，本文中首先建立了基于标准kε湍流模型的缩尺比为1∶50三维数值风场，这些数值风场满足现行的《建筑结构荷载规范》关于平均风速和湍流强度的规定；在此基础上，对TTU模型的定常绕流进行了计算，重点对阻塞率问题进行了探讨。主要结论如下：

（1）在建立小缩尺比的数值风场（相对于缩尺比1∶300），如本文中缩尺比为1∶50的B类风场，附加项是必要的；小缩尺比数值风场在低矮建筑模型的数值模拟中具有重要意义。

（2）当计算域为-5H≤x≤10H，-6H≤y≤6H，0≤z≤4H时，仍可以得到和风洞试验较为接近的数值结果；此时，与理想计算域相比，网格数由75×104下降到36×104，阻塞率由2.5%上升到7%，这个合适的计算域在计算风工程中具有重要的参考意义。

（3）对于TTU模型的屋顶和背风面，阻塞率增大可导致风压系数绝对值有增大趋势；对于TTU模型的迎风面，由计算域侧面引起的阻塞率增大，导致风压系数有减小趋势，而由计算域顶面引起的阻塞率增大，导致风压系数有增大趋势。另外，入流面和出流面位置对计算结果也会产生影响。对于TTU模型屋顶和背风面，入流面和出流面距建筑物越近，风压系数绝对值越大；对于TTU模型迎风面，出流面位置距建筑物越近，风压系数越小，而入流面位置距建筑物越近，风压系数越大。

参考文献：

[1] TOMINAGA Y，MOCHIDA A，YOSHIE R，et al.AIJ Guidelines for Practical Applications of CFD to Pedestrian Wind Environment Around Buildings[J].Journal of Wind Engineering and Industrial Aerodynamics，2008，96（10/11）：17491761.

[2]FRANKE J.Recommendations of the COST Action C14 on the Use of CFD in Predicting Pedestrian Wind Environment[C]//MURAKAMI S，MATSUMOTO M，TAMURA Y.Proceedings of the Fourth International Symposium on Computational Wind Engineering.Yokohama：CWE，2006：529532.

[3]RICHARDS P J，HOXEY R P.Appropriate Boundary Conditions for Computational Wind Engineering Models Using the Model[J].Journal of Wind Engineering and Industrial Aerodynamics，1993，46/47：145153.

[4]YANG Y，GU M，CHEN S，et al.New Inflow Boundary Conditions for Modelling the Neutral Equilibrium Atmospheric Boundary Layer in Computational Wind Engineering[J].Journal of Wind Engineering and Industrial Aerodynamics，2009，97（2）：8895.

[5]BLOCKEN B，STATHOPOULOS T，CARMELIET J.CFD Simulation of the Atmospheric Boundary Layer：Wall Function Problems[J].Atmospheric Environment，2007，41（2）：238252.

[6]FANG P Z，GU M，TAN J G，et al.Modeling the Neutral Atmospheric Boundary Layer Based on the Standard kε Turbulent Model：Modified Wall Function[C]//APCWE.Proceedings of the 7th AsiaPacific Conference on Wind Engineering.Taipei：APCWE，2009：812.

[7]方平治，顾 明，谈建国.计算风工程中基于kε系列湍流模型的数值风场[J].水动力学研究与进展：A辑，2010，25（4）：475483.

[8]LAUNDER B E，SPALDING D B.The Numerical Computation of Turbulent Flows[J].Computer Methods in Applied Mechanics and Engineering， 1974，3（2）：269289.

[9]GB 50009—2001，建筑结构荷载规范[S].

GB 50009—2001，Load Code for the Design of Building Structures[S].

[10]LEVITAN M L，MEHTA K C，VANN W P，et al.Field Measurements of Pressures on Texas Tech Building[J].Journal of Wind Engineering and Industrial Aerodynamics，1991，38（2/3）：227234.

[11]罗 攀.基于标准模型的风洞试验研究[D].上海：同济大学，2004.

[12]李加武，张宏杰，韩万水.斜拉桥风致响应的雷诺数效应[J].中国公路学报，2009，22（2）：4247.

[13]李加武，崔 欣，张宏杰，等.粗糙度对雷诺数效应的影响[J].长安大学学报：自然科学版，2009，29（2）：5659，64.

[14]汤 华，王松帆.基于风洞试验的高层建筑风荷载研究及抗风设计[J].建筑科学与工程学报，2012，29（2）：96101.