邵建濤，劉 京，趙加寧，李 彪

用于建筑繞流預測的非線性渦粘模型改良

邵建濤1，2，劉 京1，3，趙加寧1，李 彪1

（1.哈爾濱工業(yè)大學市政環(huán)境工程學院，150090哈爾濱；2.華東建筑設(shè)計研究院有限公司，200002上海；3.哈爾濱工業(yè)大學城市水資源與水環(huán)境國家重點實驗室，150090哈爾濱）

為了改良非線性渦粘性模型模擬建筑繞流問題的表現(xiàn)，首先介紹了采用渦粘性模型預測建筑繞流問題的現(xiàn)狀，分析了模擬中存在的問題，然后基于Craft非線性渦粘性模型提出了一種用于預測建筑繞流的改良的非線性渦粘性模型，并利用日本建筑學會提供的風洞實驗數(shù)據(jù)對改良的非線性模型進行分析驗證.結(jié)果表明：改良的非線性渦粘性模型一方面改善了標準k-ε模型建筑前端湍動動能預測過大的問題，預測出了建筑頂部的分離和再附著；另一方面通過增大尾跡區(qū)的湍動動能，改善了渦粘性模型在預測建筑尾跡區(qū)流動中的表現(xiàn).通過改良，非線性渦粘性模型可以較好地預測建筑風環(huán)境.

風環(huán)境；建筑繞流；非線性渦粘性模型；計算流體力學；改良

據(jù)聯(lián)合國人口調(diào)查局統(tǒng)計，2000年時，全球47%的人口居住在城市；預計到2015年，全球城市人口比例將會達52%，且大部分城市化發(fā)生在發(fā)展中國家［1］.伴隨我國經(jīng)濟的發(fā)展，城市化進程的加快，城市人口占總?cè)丝诘谋戎卦絹碓酱螅鞘忻娣e越來越大.城市的熱環(huán)境和風環(huán)境直接影響著城市居民的健康.為此，人們投入了大量的精力研究城市氣候.隨著計算機計算能力的提高和CFD模擬軟件的發(fā)展，在設(shè)計階段，CFD技術(shù)越來越多地用來預測建筑風環(huán)境.建筑風環(huán)境模擬的核心是建筑繞流問題的模擬.多年來，許多學者對建筑繞流問題作了大量的研究［2-6］；然而，建筑繞流問題的模擬長期以來一直是計算風工程領(lǐng)域的難點和熱點.對模擬建筑繞流的難點進行了分析，認為流動的高雷諾數(shù)、建筑前端的撞擊、鈍體尖銳的邊緣和出入口邊界條件是建筑繞流計算中的難點.由于對建筑前端湍流動能預測過大，標準kε模型被認為不適合用來模擬建筑繞流.文獻［7］通過研究發(fā)現(xiàn)修正的k-ε模型（如LK、MMK模型）能夠很好地修正標準k-ε模型在建筑前端高估湍動動能的問題.然而，LK模型和MMK模型存在對建筑再附著長度預測過長的缺點.

由于非線性模擬能夠克服常規(guī)湍流模型的湍流各向同性的缺點，近年來，有研究者嘗試采用非線性湍流模式來模擬建筑繞流問題.Wright等嘗試用非線性湍流模型模擬建筑繞流問題，只有一個非線性的雷諾應(yīng)力模擬得到收斂解，其他模型很難得到收斂解［8］.邵建濤等嘗試利用非線性kε模型對建筑繞流進行了模擬［9］，發(fā)現(xiàn)非線性kε模型對建筑繞流問題不能得到穩(wěn)態(tài)收斂解；而利用非穩(wěn)態(tài)求解能夠得到建筑尾跡區(qū)橫向速度的周期性波動.盡管，利用Craft等提出的非線性kε模型采用URANS方法，可以得到較標準k-ε模型、LK模型、MMK模型更為準確的風環(huán)境預測.但是，URANS方法計算時間較長，很難在存在較大計算區(qū)域的現(xiàn)實工程中廣泛應(yīng)用.本文將嘗試通過風洞實驗數(shù)據(jù)對現(xiàn)有的非線性k-ε模型進行校準和改良，以使其能夠更為快速和準確的預測建筑風環(huán)境.

1 CFD控制方程和湍流模型

1.1 控制方程

不可壓縮流體流動的控制方程為

其中:ui為xi方向的瞬時速度，m／s；xi為空間坐標，m；i，j=1，2，3表示空間坐標系的3個方向；p為瞬時壓力，Pa；ν為粘性系數(shù)，m2／s；ρ為流體密度，kg／m3.

在雷諾平均模型中，瞬時速度ui可以雷諾分解成平均速度和脈動速度兩部分，即

其中:Ui為平均速度，m／s；ui′為脈動速度，m／s．

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

1.2 非線性渦粘性模型

為了更好地在渦粘性模型內(nèi)描述雷諾應(yīng)力，非線性渦粘性模型通常將雷諾應(yīng)力表示為平均變形率張量和平均渦量的高階形式.下式是一個典型的三階非線性渦粘性模型的雷諾應(yīng)力表達式.等號的右端前兩項為線性模型的表達式，等號右端前5項為二階非線性模型的表達式.

式中:νt采用標準k-ε模型同樣的計算方法，即

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

非線性k-ε模型中湍動動能和湍動動能耗散率的控制方程及其經(jīng)驗系數(shù)的取值和標準k-ε模型相應(yīng)的控制方程及其經(jīng)驗系數(shù)的取值相同．Sij和Ωij分別為平均變形率張量和平均渦量，1／s，可由下式計算得到．

非線性渦粘性湍流模型大都十分類似，如Suga模型和Craft模型［10］，Shih模型［11］.各種非線性渦粘性模型的主要區(qū)別在于式（6）中C1到C7的非線性項系數(shù)選擇和式（7）中的cμ不同．

1.3 非線性渦粘性模型的改進

建筑繞流計算的難點在于正確預測建筑周圍的湍動動能分布.由Craft和Shih非線性渦粘性湍流模型的表達式可知，這兩個模型都對標準k-ε模型中的cμ系數(shù)進行了修正，修正了撞擊區(qū)湍動動能預測過大的問題．但根據(jù)前期研究結(jié)果可知，這兩種非線性渦粘性湍流模型對cμ的修正的同時會導致計算得到的尾跡區(qū)湍流粘性過低，穩(wěn)態(tài)計算很難達到收斂［9］．本研究的主要工作是對cμ進行進一步適當?shù)母倪M，以達到改進撞擊區(qū)湍動動能預測過大的問題，同時增大尾跡區(qū)的湍動動能，一方面補償渦脫落造成的湍動增強，另一方面加強計算的穩(wěn)定性．通過對實驗數(shù)據(jù)的分析和大量的數(shù)值試驗，在Craft模型的基礎(chǔ)上，本文提出cμ表達式為

本文改良的非線性渦粘性模型計算雷諾應(yīng)力的經(jīng)驗系數(shù)和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為特征速度的雷諾數(shù)約為24 000．

為了避免計算區(qū)域邊界對模擬結(jié)果產(chǎn)生影響，本算例所采用的計算區(qū)域大小為21b×11b× 11.25b，如圖1所示.由于一階迎風差分的數(shù)值粘性較大，容易造成較大的計算誤差，計算中所有方程中的對流項均采用二階迎風差分.壓力速度耦合采用SIMPLE算法.計算采用立方體網(wǎng)格，如圖2所示.網(wǎng)格的最小間距為0.05H，出現(xiàn)在建筑壁面處，總網(wǎng)格數(shù)約為30萬.計算采用實驗數(shù)據(jù)和平衡湍流假設(shè)設(shè)置入口邊界；其中入口邊界的速度和湍動動能分布，按照風洞實驗數(shù)據(jù)擬合曲線給定；湍動動能耗散率邊界采用假設(shè)平衡態(tài)湍流.出口邊界采用壓力出口，即各物理量沿x方向的梯度為0.建筑表面和地面采用壁面函數(shù)，其中建筑表面按光滑表面處理，地面的粗糙長度z0約為0.005.上表面和側(cè)面采用對稱邊界條件，各物理量在邊界的法向方向均為0，在切線方向梯度為0.

圖1 計算區(qū)域

圖2 計算網(wǎng)格

為了使實驗結(jié)果更加可靠，本研究對網(wǎng)格無關(guān)性進行了檢驗.

2.2 計算結(jié)果

為了便于討論，本文對標準k-ε模型、Craft非線性渦粘模型、本文改良模型3種模型的計算結(jié)果和實驗結(jié)果進行對比討論.

本文首先對預測結(jié)果中的速度場進行討論．圖3為建筑繞流再附著長度示意圖．表1為不同湍流模型的再附著長度模擬結(jié)果．圖4為y／b=0截面的流線圖．從圖中可以看出，除標準k-ε模型外，其他模型由于改善了撞擊區(qū)湍動動能預測過大的問題，都可不同程度地預測出來建筑頂部的渦．和其他研究一樣，標準k-ε模型沒能預測出建筑頂部氣流剝離和再貼附形成的渦旋．

圖3 再附著長度示意圖

表1 不同湍流模型的再附著長度模擬結(jié)果

圖4 y／b=0截面的流線圖

對建筑尾跡區(qū)而言，標準k-ε模型和Craft模型預測的建筑尾跡區(qū)再附著長度過長.本文改良模型預測的尾跡區(qū)在附著長度和風洞實驗值較為接近.

圖5為y／b=0截面上多條豎線上的速度分布．在建筑前，x／b=-0.75處，標準k-ε模型、Craft模型和本文改良模型都能較好地預測建筑前的速度分布，本文改進模型和風洞實驗數(shù)據(jù)更為接近.

圖5 y／b=0截面上多條豎線上的速度分布

而在建筑頂部x／b=0，由于標準k-ε模型沒有預測出建筑頂部的渦，所以標準k-ε模型沒有預測到建筑頂部的回流．而Craft模型和本文改進模型都成功預測出建筑頂部的建筑回流.同時，標準k-ε模型和Craft模型都低估了z=0.2 m處附近的風速，而本文改良模型成功捕捉到了該處的最大風速.在建筑尾跡區(qū)，x／b=1.25處，由于Craft模型預測建筑回流速度過大，標準k-ε模型和本文改良模型預測該處的風速和風洞實驗值較為接近.在再附著點附近，x／b=2處，標準k-ε模型和Craft模型都還能預測出地面附近的回流，這和風洞實驗的結(jié)果不符；而本文改良模型預測結(jié)果和風洞實驗結(jié)果非常接近，沒有出現(xiàn)地面附近的回流，成功預測除了流動進入到邊界層的再發(fā)展區(qū).

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

圖6為y／b=0截面的湍動動能分布云圖.從圖中可以看出，如文獻中指出，標準k-ε模型在建筑前端預測湍動動能過大，在建筑尾跡區(qū)預測湍動過小.通過對計算湍動應(yīng)力計算方法的改進，考慮到湍動的各項異性，Craft模型在一定程度上修正了在標準k-ε模型在建筑前端湍動動能預測過大的影響；然而這種修正也同時造成了尾跡區(qū)湍動動能預測更小.這也導致Craft模型成功預測了建筑頂部的渦，而惡化了建筑尾跡區(qū)渦大小的預測.本文改良模型由于在Craft模型的基礎(chǔ)上改進了湍動應(yīng)力計算方法，同時成功地修正了建筑前端湍動動能預測過大的問題和尾跡區(qū)湍動動能預測過小的問題，提高了湍動動能的預測精度.

圖7為y／b=0截面上多條豎線上的湍動動能分布．從圖中可以看出，在建筑前端，x／b= -0.75處，標準k-ε模型預測建筑前端湍動動能最大；Craft模型改善了標準k-ε模型的預測結(jié)果，但預測的湍動動能較實驗值仍然較大；本文改良模型預測的建筑前端湍動動能值和風洞實驗值最為接近，有較好的預測精度.在建筑頂部，x／b=-0.75處，3種湍流模型都不能很好的預測湍動動能的分布.從圖中可以看出，在建筑屋面附近，標準k-ε模型能夠預測湍動動能的最大值，而Craft模型和本文改良模型都不能預測出該處的湍動動能最大值；在z=0.2 m處附近，本文改良模型能夠很好的預測湍動動能大小，而標準k -ε模型和Craft模型都低估了該處湍動動能的大小.在建筑尾跡區(qū)，x／b=1.25處，Craft模型嚴重低估了該處的湍動動能大小，而標準k-ε模型在一定程度上低估了該處的湍動動能大小，本文改良模型稍微高估了該處的湍動動能大小.在x／b=2處，標準k-ε模型和Craft模型都低估了該處的湍動動能大小，而本文改良模型由于采用了新的湍動應(yīng)力的計算方法，在該處成功準確的預測了該處的湍動動能大小.

圖7 y／b=0截面上多條豎線上的湍動動能分布

2.3 結(jié)果分析

從以上的結(jié)果可以看出，由于本文改良的非線性模型考慮了湍流的各項異性，同時采用了新的湍流應(yīng)力計算方法，本文改良模型一方面修正了標準k-ε模型建筑前端湍動動能預測過大的問題，同時又增大了建筑尾跡區(qū)的湍動動能．因此，本文改良模型一方面預測出了建筑頂部漩渦，另一方面較為準確地預測了建筑尾跡區(qū)的再附著長度．

從建筑繞流的流動特性而言，標準k-ε模型對建筑前端的湍動動能預測過大是由于標準k-ε模型對湍流應(yīng)力的預測具有各項同性特點．很多對標準k-ε模型的修正都可以解決這個問題．但是，大部分修正的k-ε模型都會惡化對建筑尾跡區(qū)的預測，使預測再附著長度更長．事實上，建筑尾跡區(qū)存在著劇烈的非穩(wěn)態(tài)流動，大量的渦脫落使得該部分混合較為劇烈．大部分渦粘性不能預測出該部分的渦脫落，因此對再附著長度的預測都顯得過長．本文改良模型通過增大尾跡區(qū)的湍動動能，加強湍動擴散以求達到和渦脫落相似的混合效果，取得較為準確的尾跡區(qū)的流場．

3 結(jié) 論

1）建筑繞流問題是建筑環(huán)境模擬中的重要問題.基于Craft模型提出了一種用于建筑繞流預測的改良的非線性渦粘性模型，并利用日本建筑學會提供的實驗數(shù)據(jù)對其進行了驗證.

2）改良了模型，一方面修正了標準k-ε模型建筑前端湍動動能預測過大的問題，預測出了建筑頂部的分離和再附著；另一方面通過增大尾跡區(qū)的湍動動能，改善了渦粘性模型在預測建筑尾跡區(qū)流動中的表現(xiàn).改良模型預測的建筑尾跡區(qū)再附著長度和風洞實驗數(shù)據(jù)非常接近.

［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與建筑環(huán)境［M］.朱清宇，譯.北京:中國建筑工業(yè)出版社，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.