用于建筑绕流预测的非线性涡粘模型改良

邵建涛，刘 京，赵加宁，李 彪

用于建筑绕流预测的非线性涡粘模型改良

邵建涛1，2，刘 京1，3，赵加宁1，李 彪1

（1.哈尔滨工业大学市政环境工程学院，150090哈尔滨；2.华东建筑设计研究院有限公司，200002上海；3.哈尔滨工业大学城市水资源与水环境国家重点实验室，150090哈尔滨）

为了改良非线性涡粘性模型模拟建筑绕流问题的表现，首先介绍了采用涡粘性模型预测建筑绕流问题的现状，分析了模拟中存在的问题，然后基于Craft非线性涡粘性模型提出了一种用于预测建筑绕流的改良的非线性涡粘性模型，并利用日本建筑学会提供的风洞实验数据对改良的非线性模型进行分析验证.结果表明：改良的非线性涡粘性模型一方面改善了标准k-ε模型建筑前端湍动动能预测过大的问题，预测出了建筑顶部的分离和再附着；另一方面通过增大尾迹区的湍动动能，改善了涡粘性模型在预测建筑尾迹区流动中的表现.通过改良，非线性涡粘性模型可以较好地预测建筑风环境.

风环境；建筑绕流；非线性涡粘性模型；计算流体力学；改良

据联合国人口调查局统计，2000年时，全球47%的人口居住在城市；预计到2015年，全球城市人口比例将会达52%，且大部分城市化发生在发展中国家［1］.伴随我国经济的发展，城市化进程的加快，城市人口占总人口的比重越来越大，城市面积越来越大.城市的热环境和风环境直接影响着城市居民的健康.为此，人们投入了大量的精力研究城市气候.随着计算机计算能力的提高和CFD模拟软件的发展，在设计阶段，CFD技术越来越多地用来预测建筑风环境.建筑风环境模拟的核心是建筑绕流问题的模拟.多年来，许多学者对建筑绕流问题作了大量的研究［2-6］；然而，建筑绕流问题的模拟长期以来一直是计算风工程领域的难点和热点.对模拟建筑绕流的难点进行了分析，认为流动的高雷诺数、建筑前端的撞击、钝体尖锐的边缘和出入口边界条件是建筑绕流计算中的难点.由于对建筑前端湍流动能预测过大，标准kε模型被认为不适合用来模拟建筑绕流.文献［7］通过研究发现修正的k-ε模型（如LK、MMK模型）能够很好地修正标准k-ε模型在建筑前端高估湍动动能的问题.然而，LK模型和MMK模型存在对建筑再附着长度预测过长的缺点.

由于非线性模拟能够克服常规湍流模型的湍流各向同性的缺点，近年来，有研究者尝试采用非线性湍流模式来模拟建筑绕流问题.Wright等尝试用非线性湍流模型模拟建筑绕流问题，只有一个非线性的雷诺应力模拟得到收敛解，其他模型很难得到收敛解［8］.邵建涛等尝试利用非线性kε模型对建筑绕流进行了模拟［9］，发现非线性kε模型对建筑绕流问题不能得到稳态收敛解；而利用非稳态求解能够得到建筑尾迹区横向速度的周期性波动.尽管，利用Craft等提出的非线性kε模型采用URANS方法，可以得到较标准k-ε模型、LK模型、MMK模型更为准确的风环境预测.但是，URANS方法计算时间较长，很难在存在较大计算区域的现实工程中广泛应用.本文将尝试通过风洞实验数据对现有的非线性k-ε模型进行校准和改良，以使其能够更为快速和准确的预测建筑风环境.

1 CFD控制方程和湍流模型

1.1 控制方程

不可压缩流体流动的控制方程为

其中:ui为xi方向的瞬时速度，m／s；xi为空间坐标，m；i，j=1，2，3表示空间坐标系的3个方向；p为瞬时压力，Pa；ν为粘性系数，m2／s；ρ为流体密度，kg／m3.

在雷诺平均模型中，瞬时速度ui可以雷诺分解成平均速度和脉动速度两部分，即

其中:Ui为平均速度，m／s；ui′为脉动速度，m／s．

把式（3）分别代入式（1）、（2）后，可以得到

1.2 非线性涡粘性模型

为了更好地在涡粘性模型内描述雷诺应力，非线性涡粘性模型通常将雷诺应力表示为平均变形率张量和平均涡量的高阶形式.下式是一个典型的三阶非线性涡粘性模型的雷诺应力表达式.等号的右端前两项为线性模型的表达式，等号右端前5项为二阶非线性模型的表达式.

式中:νt采用标准k-ε模型同样的计算方法，即

其中:k为湍动动能，m／s；ε为湍动动能耗散率，m2／s3.

非线性k-ε模型中湍动动能和湍动动能耗散率的控制方程及其经验系数的取值和标准k-ε模型相应的控制方程及其经验系数的取值相同．Sij和Ωij分别为平均变形率张量和平均涡量，1／s，可由下式计算得到．

非线性涡粘性湍流模型大都十分类似，如Suga模型和Craft模型［10］，Shih模型［11］.各种非线性涡粘性模型的主要区别在于式（6）中C1到C7的非线性项系数选择和式（7）中的cμ不同．

1.3 非线性涡粘性模型的改进

建筑绕流计算的难点在于正确预测建筑周围的湍动动能分布.由Craft和Shih非线性涡粘性湍流模型的表达式可知，这两个模型都对标准k-ε模型中的cμ系数进行了修正，修正了撞击区湍动动能预测过大的问题．但根据前期研究结果可知，这两种非线性涡粘性湍流模型对cμ的修正的同时会导致计算得到的尾迹区湍流粘性过低，稳态计算很难达到收敛［9］．本研究的主要工作是对cμ进行进一步适当的改进，以达到改进撞击区湍动动能预测过大的问题，同时增大尾迹区的湍动动能，一方面补偿涡脱落造成的湍动增强，另一方面加强计算的稳定性．通过对实验数据的分析和大量的数值试验，在Craft模型的基础上，本文提出cμ表达式为

本文改良的非线性涡粘性模型计算雷诺应力的经验系数和Craft非线性模型相同.

本模型中C1～C7分别为-0.1、-0.1、

2 算例分析

2.1 算例概况

本文将采用文献中最为常见的日本建筑学会标准算例中的1∶1∶2型建筑为例对本文改良的非线性涡粘性模型进行验证，该建筑如图1所示［12］；建筑物的长度和宽度b均为0.08 m，高度H为0.16 m.以建筑高度H为特征长度和建筑高度处风速UH为特征速度的雷诺数约为24 000．

为了避免计算区域边界对模拟结果产生影响，本算例所采用的计算区域大小为21b×11b× 11.25b，如图1所示.由于一阶迎风差分的数值粘性较大，容易造成较大的计算误差，计算中所有方程中的对流项均采用二阶迎风差分.压力速度耦合采用SIMPLE算法.计算采用立方体网格，如图2所示.网格的最小间距为0.05H，出现在建筑壁面处，总网格数约为30万.计算采用实验数据和平衡湍流假设设置入口边界；其中入口边界的速度和湍动动能分布，按照风洞实验数据拟合曲线给定；湍动动能耗散率边界采用假设平衡态湍流.出口边界采用压力出口，即各物理量沿x方向的梯度为0.建筑表面和地面采用壁面函数，其中建筑表面按光滑表面处理，地面的粗糙长度z0约为0.005.上表面和侧面采用对称边界条件，各物理量在边界的法向方向均为0，在切线方向梯度为0.

图1 计算区域

图2 计算网格

为了使实验结果更加可靠，本研究对网格无关性进行了检验.

2.2 计算结果

为了便于讨论，本文对标准k-ε模型、Craft非线性涡粘模型、本文改良模型3种模型的计算结果和实验结果进行对比讨论.

本文首先对预测结果中的速度场进行讨论．图3为建筑绕流再附着长度示意图．表1为不同湍流模型的再附着长度模拟结果．图4为y／b=0截面的流线图．从图中可以看出，除标准k-ε模型外，其他模型由于改善了撞击区湍动动能预测过大的问题，都可不同程度地预测出来建筑顶部的涡．和其他研究一样，标准k-ε模型没能预测出建筑顶部气流剥离和再贴附形成的涡旋．

图3 再附着长度示意图

表1 不同湍流模型的再附着长度模拟结果

图4 y／b=0截面的流线图

对建筑尾迹区而言，标准k-ε模型和Craft模型预测的建筑尾迹区再附着长度过长.本文改良模型预测的尾迹区在附着长度和风洞实验值较为接近.

图5为y／b=0截面上多条竖线上的速度分布．在建筑前，x／b=-0.75处，标准k-ε模型、Craft模型和本文改良模型都能较好地预测建筑前的速度分布，本文改进模型和风洞实验数据更为接近.

图5 y／b=0截面上多条竖线上的速度分布

而在建筑顶部x／b=0，由于标准k-ε模型没有预测出建筑顶部的涡，所以标准k-ε模型没有预测到建筑顶部的回流．而Craft模型和本文改进模型都成功预测出建筑顶部的建筑回流.同时，标准k-ε模型和Craft模型都低估了z=0.2 m处附近的风速，而本文改良模型成功捕捉到了该处的最大风速.在建筑尾迹区，x／b=1.25处，由于Craft模型预测建筑回流速度过大，标准k-ε模型和本文改良模型预测该处的风速和风洞实验值较为接近.在再附着点附近，x／b=2处，标准k-ε模型和Craft模型都还能预测出地面附近的回流，这和风洞实验的结果不符；而本文改良模型预测结果和风洞实验结果非常接近，没有出现地面附近的回流，成功预测除了流动进入到边界层的再发展区.

图6 y／b=0截面的湍动动能分布云图（m2／s2）

图6为y／b=0截面的湍动动能分布云图.从图中可以看出，如文献中指出，标准k-ε模型在建筑前端预测湍动动能过大，在建筑尾迹区预测湍动过小.通过对计算湍动应力计算方法的改进，考虑到湍动的各项异性，Craft模型在一定程度上修正了在标准k-ε模型在建筑前端湍动动能预测过大的影响；然而这种修正也同时造成了尾迹区湍动动能预测更小.这也导致Craft模型成功预测了建筑顶部的涡，而恶化了建筑尾迹区涡大小的预测.本文改良模型由于在Craft模型的基础上改进了湍动应力计算方法，同时成功地修正了建筑前端湍动动能预测过大的问题和尾迹区湍动动能预测过小的问题，提高了湍动动能的预测精度.

图7为y／b=0截面上多条竖线上的湍动动能分布．从图中可以看出，在建筑前端，x／b= -0.75处，标准k-ε模型预测建筑前端湍动动能最大；Craft模型改善了标准k-ε模型的预测结果，但预测的湍动动能较实验值仍然较大；本文改良模型预测的建筑前端湍动动能值和风洞实验值最为接近，有较好的预测精度.在建筑顶部，x／b=-0.75处，3种湍流模型都不能很好的预测湍动动能的分布.从图中可以看出，在建筑屋面附近，标准k-ε模型能够预测湍动动能的最大值，而Craft模型和本文改良模型都不能预测出该处的湍动动能最大值；在z=0.2 m处附近，本文改良模型能够很好的预测湍动动能大小，而标准k -ε模型和Craft模型都低估了该处湍动动能的大小.在建筑尾迹区，x／b=1.25处，Craft模型严重低估了该处的湍动动能大小，而标准k-ε模型在一定程度上低估了该处的湍动动能大小，本文改良模型稍微高估了该处的湍动动能大小.在x／b=2处，标准k-ε模型和Craft模型都低估了该处的湍动动能大小，而本文改良模型由于采用了新的湍动应力的计算方法，在该处成功准确的预测了该处的湍动动能大小.

图7 y／b=0截面上多条竖线上的湍动动能分布

2.3 结果分析

从以上的结果可以看出，由于本文改良的非线性模型考虑了湍流的各项异性，同时采用了新的湍流应力计算方法，本文改良模型一方面修正了标准k-ε模型建筑前端湍动动能预测过大的问题，同时又增大了建筑尾迹区的湍动动能．因此，本文改良模型一方面预测出了建筑顶部漩涡，另一方面较为准确地预测了建筑尾迹区的再附着长度．

从建筑绕流的流动特性而言，标准k-ε模型对建筑前端的湍动动能预测过大是由于标准k-ε模型对湍流应力的预测具有各项同性特点．很多对标准k-ε模型的修正都可以解决这个问题．但是，大部分修正的k-ε模型都会恶化对建筑尾迹区的预测，使预测再附着长度更长．事实上，建筑尾迹区存在着剧烈的非稳态流动，大量的涡脱落使得该部分混合较为剧烈．大部分涡粘性不能预测出该部分的涡脱落，因此对再附着长度的预测都显得过长．本文改良模型通过增大尾迹区的湍动动能，加强湍动扩散以求达到和涡脱落相似的混合效果，取得较为准确的尾迹区的流场．

3 结 论

1）建筑绕流问题是建筑环境模拟中的重要问题.基于Craft模型提出了一种用于建筑绕流预测的改良的非线性涡粘性模型，并利用日本建筑学会提供的实验数据对其进行了验证.

2）改良了模型，一方面修正了标准k-ε模型建筑前端湍动动能预测过大的问题，预测出了建筑顶部的分离和再附着；另一方面通过增大尾迹区的湍动动能，改善了涡粘性模型在预测建筑尾迹区流动中的表现.改良模型预测的建筑尾迹区再附着长度和风洞实验数据非常接近.

［1］United Nations，Department of Economic and Social Affairs.World urbanization prospects:the 2009 revision［EB／OL］.［2009-03-01］.http:／／esa.un.org／unpd／wup／CD-ROM-2009／WUP2009-F03-Urban-Population. xls.

［2］STATHOPOULOST，BASKARANB.Computer simulation of wind environmental conditions around buildings［J］.Engineering Structures，1996，18（11）:876-885.

［3］TSUCHIYA M，MURAKAMI S，MOCHIDA A，et al. Development of a new k-ε model for flow and pressure fields aroundbluffbody［J］.JournalofWind Engineering and Industrial Aerodynamics，1997，67-68:169-182.

［4］MURAKAMI S，OOKA R，MOCHIDA A，et al.CFD analysis of wind climate from human scale to urban scale［J］.Journal of Wind Engineering and Industrial Aerodynamics，1999，81（1／2／3）:57-81.

［5］WANG X，MCNAMARA K F.Evaluation of CFD simulation using RANS turbulence models for building effects on pollution dispersion［J］.Environmental Fluid Mechanics，2006，6（2）:181-202.

［6］村上周三.CFD与建筑环境［M］.朱清宇，译.北京:中国建筑工业出版社，2007.

［7］MURAKAMI S.Current status and future trends in computational wind engineering［J］.Journal of Wind Engineering and Industrial Aerodynamics，1997，67-68:3-34.

［8］WRIGHT N G，EASOM G J.Non-linear turbulence model results for flow over a building at full-scale［J］. Applied Mathematical Modelling，2003，27（12）: 1013-1033.

［9］SHAO Jiantao，LIU Jing，ZHAO Jianing.Evaluation of various non-linear k-ε models for predicting wind flow around an isolated high-rise building within the surface boundary layer［J］.Building and Environment，2012，57（11）:145-155.

［10］CRAFT T J，LAUNDER B E，SUGA K.Development and application of a cubic eddy-viscosity model of turbulence［J］.International Journal of Heat and Fluid Flow，1996，17（2）:108-115.

［11］SHIH T H，ZHU J，LUMLEY J L.A realizable reynolds stress algebraicequationmodel［R］.Cleveland: NASA，1993.

［12］TOMINAGA Y，MOCHIDA A，MURAKAMI S，et al. Comparison of various revised k-ε models and LES applied to flow around a high-rise building model with 1∶1∶2 shape placed within the surface boundary layer［J］.Journal of Wind Engineering and Industrial Aerodynamics，2008，96（4）:389-411.

（编辑 魏希柱）

Improvement of the non-linear eddy viscosity model applied to predicting wind flow around building

SHAO Jiantao1，2，LIU Jing1，3，ZHAO Jianing1，LI Biao1
（1.School of Municipal and Environmental Engineering，Harbin Institute of Technology，150090 Harbin，China；2.East China Architectural Design＆Research Institute Co.Ltd.，200002 Shanghai，China；3.State Key Laboratory of Urban Water Resource and Environment，Harbin Institute of Technology，150090 Harbin，China）

The aim of this paper is to improve the performance of the non-linear eddy viscosity model for simulating the wind flow around the building.Firstly，the state of art of the predicting wind flow around building using RANS model was introduced，and the problems in the simulation were analyzed.Then an improved non-linear eddy viscosity for predicting the wind flow around buildings was proposed based on Craft model.The improved non-linear eddy viscosity was validated and analyzed through the wind tunnel data provided by AIJ.The results showed that the proposed non-linear eddy viscosity improved the overestimation of turbulent kinetic energy in impingement region by the standard k-ε model，and predicted better results in the wake region behind buildings simultaneously through strengthening the eddy viscosity in the wake region.After the improvement，the non-linear eddy viscosity model can predict the wind environment around buildings better.

wind environment；flow around building；non-linear eddy viscosity model；computational fluid dynamics；improvement

TU111.19+3

A

0367-6234（2014）04-0050-07

2013-03-30.

国家自然科学基金资助项目（40505025）；

城市气象科学研究基金资助项目（UMRF201004）.

邵建涛（1983—），男，博士，高级工程师；

刘 京（1972—），男，教授，博士生导师；

赵加宁（1956—），女，教授，博士生导师.

邵建涛，shaojiantao＠gmail.com.