基于拟周期方法的高层建筑预前模拟

邱顺冬 程 旭

(1.建研科技股份有限公司，北京 100031； 2.西南交通大学，四川 成都 610031)

0 引言

由于设计方法及新建筑材料和技术的创新发展所驱动，新一代高层建筑正在变得越来越苗条、高柔，这就使得高层建筑对诸如风等横向载荷变得更加敏感。此外，为了控制地震产生的惯性力，通常需要降低建筑物的重量，这又进一步有助于高层建筑风荷载及其导致的高层建筑位移的增加。因此，风荷载越来越成为高层建筑横向结构设计主要考虑的因素。目前，高层建筑风荷载的评估主要有风洞试验[1]和CFD(Computational Fluid Dynamics)模拟两种方法。虽然风洞试验便于改变和重复试验条件，然而成本高、受试验条件影响及其所测试验数据有限等缺点严重限制了其发展。另外一方面，凭借其在硬件成本和运行耗能等方面的优势以及流场全貌显示等特点，关于使用CFD模拟来评估高层建筑风荷载的研究变得越来越广泛。

CFD模拟的湍流模型主要可以分为RANS(Reynold-Averaged Navier-Stokes)和LES(Large Eddy Simulation)两种。RANS在平均风特性上表现良好，但其很难预测脉动风成分。而LES方法通过滤波函数将大尺度的涡和小尺度的涡分离开，大尺度的涡直接模拟，小尺度的涡用模型来封闭，这就使得LES从某种程度上可以直接计算得到脉动风成分，故其也正在成为越来越流行的CFD计算方法。而要使用LES进行风荷载的准确评估，入口边界条件则成为一个十分重要的影响因素。目前针对入口边界条件的LES(Large Eddy Simulation)数值模拟方法主要可以分为两大类：人工合成湍流法和预前模拟法。

人工合成湍流法主要是将产生的波动速度叠加目标平均速度。波动速度的产生常用的主要有以下几种方法：谱合成、数字滤波法及拟序结构叠加。为此，国内外很多学者在此基础上进行了改进，如Smirnov et al[2]，Batten et al[3]，Huang et al[4]，Mathey et al[5]。这种方法的优点是产生非均匀湍流,适合复杂的入口。但其很难满足连续性方程，这就导致入流风特性进入计算域内部以后会发生变化，比如在建筑物附近不能达到目的湍流度等，这就给实际应用带来了困难。

预前模拟是在被研究建筑之前通过放置尖劈及粗糙元等钝体，被动模拟出目标风场。此方法最大的优势在于其模拟出的风场满足多点相关性及连续性方程，近似于真实的湍流，对模拟钝体绕流具有十分重要的意义。然而为了模拟出充分发展的湍流，且获得足够的统计数据，往往需要大量的计算网格和计算时间，这也是预前模拟法在工程无法普遍应用的重要原因。由此，为了减少存储量，Spalart[6]提出了周期边界，Lund et al[7]对其进行了改进。Kataoka[8]对Lund et al[7]的方法进行了简化，提出了拟周期边界条件。

本文主要针对拟周期边界条件提出了用于粗糙元设计方法，并且通过大量的LES模拟，验证了公式的正确性，为此方法的应用提供了参考。

1 数值建模及求解设置

1.1 控制方程

大涡模拟(LES)湍流模型控制方程可以表达为：

(1)

(2)

1.2 计算域及边界条件

本模拟的计算域如图1所示，其中计算域尺寸为2 200 m(顺风向)×550 m(横风向)×550 m(竖向)。此外，计算域底部布置间距为100 m的粗糙块，顺风和横风的尺寸均为15 m，而粗糙块高度的设置与风场类别有关。

入口采用Kataoka[8]提出的拟周期边界，其中循环面与入口的距离为800 m。另外，计算域底部以及粗糙块采用无滑移壁面，计算域顶部采用滑移壁面，出口采用压力出口，侧面采用周期边界。注意到，本文采用谐波合成的方法用于增大远地表处的湍流度。

流场的入口速度分量如式(3)所示：

uinlt(y,z,t)=U(z)+φ(θ)×{u(y,z,t)-U(y,z)}recyvinlt(y,z,t)=φ(θ)×{v(y,z,t)-V(y,z)}recy
winlt(y,z,t)=φ(θ)×{w(y,z,t)-W(y,z)}recy

(3)

其中，U,V和W分别为顺流向，横流向和竖直方向的平均速度；u(y,z,t)，v(y,z,t)，w(y,z,t)均为各个方向的瞬时速度；下标的inlt和recy分别为入口和循环面上提取的数值。φ(θ)的表达式如式(4)所示，是由Kataoka[8]提出的阻尼函数。

(4)

其中，θ=y/Ly,Ly为大气边界层高度。而权重数φ(θ)≤1.0，作用是稳定循环后的速度波动，达到目标值。

2 结果与讨论

2.1 LES模拟结果分析

本次模拟的参数主要选取了平均风剖面、湍流度剖面，这两者是通过对大气边界层风场特征参数和规范对比确定的，模拟的目标分别为中国规范指数律风剖面和日本规范湍流度剖面。如图2，图3所示，平均风速剖面和湍流度剖面的模拟结果与所对应的规范值基本吻合，证明了本文采用拟周期方法的适用性。

2.2 粗糙元的设计

在以往基于LES的预前模拟法中并没有给予直接选取粗糙元大小的方法，所以本文在参考Wooding et al[10]风洞试验中所采用的公式的基础上，对风场模拟中的粗糙元进行设计，并提出相应的公式。

Wooding et al[10]在文中推导出的应用于风洞试验的经验公式如式(5)所示:

(5)

其中，δ为湍流剪切层或边界层厚度；κ为卡曼常数；k为粗糙元的高度；λ为粗糙度；φ为关于k和s的函数，φ≈(k/s)β，β为一个常数；U为z=δ处的风速；u*为总剪切速度；C为一个Wooding et al在论文中确定的常数，C≈-2.05。

根据u*=(τ/ρ)1/2，τ为总剪切应力;ρ为流体密度。从阻力系数的计算公式τ=0.5ρU2Cf可以得到总剪切速度的公式：

(6)

其中在Wooding et al的论文中，s=d，s为粗糙元平行于流场的水平维度，d为粗糙元垂直于来流方向的水平宽度，在本文中s与d相同，D为粗糙元之间的间距。

由此可以将式(5)推导成：

(7)

其中，λ=Ar/S，Ar为粗糙元的正面面积，S为每个粗糙元所占有的区域，S=D2，D为粗糙元之间的间距；β为一个常数；κ为卡曼常数，根据Irwin[11]得到κ≈0.35，β≈0.4。Wooding et al和Irwin[11]都提到了，公式的应用范围就是70<δ/kλφ<2 000。

将各个常数代入到式(7)中，可得到式(8)：

(8)

通过数值模拟的结果，本文计算出Cf的数值，采用多项式拟合出Cf和α的关系，Cf=0.3α2-0.03α+0.01。依据规范，不同的地面类型选取不同的α的数值，其中A,B,C,D四类地面所对应的α分别为0.12,0.15,0.22,0.30。根据Hagishima et al[12]得出的结论，λ<0.2时，Cf随着λ的变大而上升。地面粗糙度指数α和Cf在A,B,C,D四类地表情况下，取值都是逐渐变大，所以本文在选取λ时也要在0～0.2的区间逐渐变大。根据王婷婷[13]的LES模拟结果，计算出其模型中式(8)中包含的各个参数，代入式(8)中能够发现比较好的吻合。各参数的数值如表1所示。

表1 各类地表LES模拟得到的参数

本文还通过计算的算例检验这个公式对于粗糙元设计的指导性。本文选取B类地表，对应的地面粗糙度指数α为0.15，大气边界层的厚度为350 m。根据拟合曲线算出的Cf=0.011，对应的λ取0.1，可通过式(6)算出k1.4/d0.4=20.86。流场的正面尺寸为550 m×550 m，粗糙元的间距取为50 m，根据λ=Ar/S，得到kd=250 m2。计算得到d≈14 m，k≈19 m。根据这一数据可以把粗糙元的尺寸和布置情况大致选定，算例所采用的粗糙元尺寸顺流向和横流向的尺寸为15 m，高度取20 m。

对比结果发现，应用公式计算得到的粗糙块的尺寸和间距，LES模拟得到的大气边界层的风场特征参数与规范B类地表能够较好的吻合，如图2，图3所示。这一实例也证实了式(8)在基于LES的预前模拟法粗糙元设计中具有可操作性。

3 结语

本文对B类地表的大气边界层风场进行了LES模拟，通过与规范对比验证模拟结果的正确性，并提出了相应的粗糙块设计方法。

采用拟周期方法，通过布置粗糙块和谐波合成方法来改变计算域入口风特性设置，使得模拟得到的建筑物附近的风速剖面和湍流度剖面与规范值基本符合。进一步通过大量的模拟，提出了基于LES预前模拟法的粗糙元设计式(8)，公式给出了粗糙元的高度和宽度的关系，可以用来快速选取粗糙元的尺寸以及布置粗糙元的位置，并能够得到较好的模拟结果，对今后采用LES湍流模型模拟大气边界层给予一定的指导意义。

参考文献:

[1] Huang G, Chen X. Wind load effects and equivalent static wind loads of tall buildings based on synchronous pressure measurements[J]. Engineering Structures,2007,29(10):2641-2653.

[2] Smirnov A, Shi S, Celik I. Random flow generation technique for large eddy simulations and particle-dynamics modeling[J]. Journal of Fluids Engineering,2001,123(2):359-371.

[3] Batten P, Goldberg U, Chakravarthy S. Interfacing statistical turbulence closures with large-eddy simulation[J]. AIAA journal, 2004,42(3):485-492.

[4] Huang S H, Li Q S, Wu J R. A general inflow turbulence generator for large eddy simulation[J]. Journal of Wind Engineering and Industrial Aerodynamics,2010,98(10-11):600-617.

[5] Mathey F, Cokljat D, Bertoglio J P, et al. Assessment of the vortex method for large eddy simulation inlet conditions[J]. Progress in Computational Fluid Dynamics, An International Journal,2006,6(1-3):58-67.

[6] Spalart P R. Direct simulation of a turbulent boundary layer up toRθ=1 410[J]. Journal of fluid mechanics,1988(187):61-98.

[7] Lund T S, Wu X, Squires K D. Generation of turbulent inflow data for spatially-developing boundary layer simulations[J]. Journal of computational physics,1998,140(2):233-258.

[8] Kataoka H. Numerical simulations of a wind-induced vibrating square cylinder within turbulent boundary layer[J]. Journal of Wind Engineering and Industrial Aerodynamics,2008,96(10-11):1985-1997.

[9] Murakami S. Current status and future trends in computational wind engineering[J]. Journal of wind engineering and industrial aerodynamics,1997(67):3-34.

[10] Wooding R A, Bradley E F, Marshall J K. Drag due to regular arrays of roughness elements of varying geometry[J]. Boundary-Layer Meteorology,1973,5(3):285-308.

[11] Irwin H. Design and use of spires for natural wind simulation. National Research Council, Canada. LTR-LA-233,1979.

[12] Hagishima A, Tanimoto J, Nagayama K, et al. Aerodynamic parameters of regular arrays of rectangular blocks with various geometries[J]. Boundary-Layer Meteorology,2009,132(2):315-337.

[13] 王婷婷.基于FLUENT的大气边界层风场LES模拟[D].北京：北京交通大学结构工程硕士学位论文，2011:43-76.