方平治，顧明，談建國，韓志惠

(1.中國氣象局上海臺風研究所，上海 200030; 2.同濟大學土木工程防災(zāi)國家重點實驗室，上海 200092;3.上海市氣象科學研究所，上海 200030)

數(shù)值模擬大氣邊界層中解決壁面函數(shù)問題方法研究

方平治1，顧明2，談建國3，韓志惠3

(1.中國氣象局上海臺風研究所，上海 200030; 2.同濟大學土木工程防災(zāi)國家重點實驗室，上海 200092;3.上海市氣象科學研究所，上海 200030)

在雷諾平均N－S方程湍流模型框架內(nèi)，壁面函數(shù)常用來模型化壁面附近的低雷諾數(shù)流動。探討基于標準湍流模型數(shù)值模擬大氣邊界層中出現(xiàn)的壁面函數(shù)問題。在已有標準壁面函數(shù)基礎(chǔ)上，通過增加一個附加項來模型化地表上面建筑結(jié)構(gòu)等粗糙元由于大小不一、錯亂分布對地表附近空氣流動產(chǎn)生的附加影響。通過模擬縮尺比為1∶300的具有較大空氣動力學粗糙長度的中性大氣邊界層，以及縮尺比為1∶50的TTU低矮建筑模型在中性大氣邊界層內(nèi)的繞流，對附加項的有效性和使用場合進行評估和說明。結(jié)果表明:附加項對于解決壁面函數(shù)問題，即在計算域內(nèi)保持來流邊界條件是必要的。

大氣邊界層;標準壁面函數(shù);壁面函數(shù)問題;物理粗糙高度;空氣動力學粗糙長度

重現(xiàn)中性大氣邊界層為計算風工程的基本要求［1］。數(shù)值模擬方法精度不僅受鈍體繞流復雜性影響，數(shù)值模擬建筑結(jié)構(gòu)所在大氣邊界層的不足也是重要因素。在無建筑結(jié)構(gòu)的計算域中，源于中性大氣邊界層的來流邊界條件應(yīng)在整個計算域流向方向保持一致。Richards等［2］發(fā)現(xiàn)來流邊界條件與建筑結(jié)構(gòu)前來流間存在差異。基于標準k－ε湍流模型，Richards等［3］提出數(shù)值模擬中性大氣邊界層的4個經(jīng)典假設(shè); Richards等［4］對數(shù)值模擬中性大氣邊界層的基本理論進行更普適研究，并論述邊界條件、湍流模型及壁面函數(shù)三者間相互關(guān)系。

基于雷諾平均N－S方程的湍流模型通常采用兩種方法處理壁面附近的低雷諾數(shù)流動，即壁面函數(shù)法與近壁面模型法。前者采用標準壁面函數(shù)［5］;后者主要增加壁面附近的網(wǎng)格精度，但會增大計算量。Craft等［6－7］提出解析壁面函數(shù)。與前述兩方法不同，解析壁面函數(shù)假定壁面附近的湍流粘性發(fā)生變化，通過直接積分控制方程求解壁面附近的低雷諾數(shù)流動。Miles等［8］通過直接模擬設(shè)置在建筑結(jié)構(gòu)前的粗糙元獲得來流邊界條件。以上各種方法中，壁面函數(shù)法因在保證精度的同時兼具經(jīng)濟性、有效性，廣泛用于工程中。然而，壁面函數(shù)問題，即在計算域地面附近，來流邊界條件與建筑結(jié)構(gòu)前來流間存在差異，因在模擬具有較大空氣動力學粗糙長度的大氣邊界層中經(jīng)常出現(xiàn)而使大氣邊界層的模擬精度受到影響［9］。

數(shù)值模擬中性大氣邊界層時，相容的來流邊界條件函數(shù)表達式同樣具有重要作用。為滿足標準k－ε湍流模型的控制方程，文獻［3］的湍動能來流邊界條件的函數(shù)表達式為常數(shù)。而據(jù)文獻［10－11］研究結(jié)果，湍動能在地面附近沿高度變化并存在峰值。Yang等［12］對相容性重新研究，給出與實際觀測結(jié)果相吻合的湍動能來流邊界條件函數(shù)表達式，并用于其它基于雷諾平均N－S方程的k－ε湍流模型中［13］。

基于標準k－ε湍流模型，本文給出在模擬中性大氣邊界層中解決壁面函數(shù)問題方法。該方法強調(diào)在已有標準壁面函數(shù)基礎(chǔ)上導出的附加項的重要作用。本文方法及所有數(shù)值計算在商用CFD軟件Ansys 13中實現(xiàn)。本文在回顧CFD中常用標準壁面函數(shù)基礎(chǔ)上，導出附加項用于模型化地表面粗糙元因大小不一、錯亂分布對地表附近空氣流動產(chǎn)生的附加影響;通過分別模擬縮尺比為1∶300的具有較大空氣動力學粗糙長度的中性大氣邊界層、縮尺比為1∶50的TTU低矮建筑模型在中性大氣邊界層內(nèi)的繞流，對附加項的有效性及使用場合進行評估、說明，并對本文方法進行評價及總結(jié)。

1 標準壁面函數(shù)及附加項

在計算風工程中，通常采用基于雷諾平均N－S方程的湍流模型及壁面函數(shù)方法模擬壁面附近的低雷諾數(shù)流動。ρu*z/μ≥100時，壁面附近的流動常用對數(shù)律模擬［14］，其中z為至壁面距離，u*為摩擦風速，ρ為空氣密度，μ為空氣動力粘性系數(shù)。

對光滑壁面，對數(shù)律為

式中:U為流向平均風速;κ=0.42為馮·卡門常數(shù)［3］; C1=5.5。

對粗糙壁面，對數(shù)律為

式中:Ks為物理粗糙高度，表示粗糙度;C2為常數(shù)。

據(jù)式(1)，標準壁面函數(shù)［5］為

式中:U*為無量綱風速;kP為壁面第一層網(wǎng)格中心P點處湍動能;UP為P點風速;z*為無量綱距離;zP為P點至壁面距離;Cμ為模型常數(shù);E=9.79。

對由大小相等、分布均勻的砂粒型粗糙元形成的粗糙壁面，在標準壁面函數(shù)基礎(chǔ)上，據(jù)式(2)，修正的壁面函數(shù)［15］為

在商用CFD軟件如Ansys 13中，要求Ks≤zP。然而工程應(yīng)用中，中性大氣邊界層對應(yīng)的Ks通常大于數(shù)值計算的zP，從而導致壁面函數(shù)問題［9］。在標準壁面函數(shù)基礎(chǔ)上，通過引入新的附加項δB，可一定程度上解決該問題。γ=1時式(9)退化為式(4)。附加項δB也可認為是ΔB的一部分，故δB無新意，但式(9)提供了解決壁面函數(shù)問題的思路。

2 具有較大空氣動力學粗糙長度的中性大氣邊界層模擬

據(jù)我國建筑結(jié)構(gòu)荷載規(guī)范［16］，同濟大學在TJ－2大氣邊界層風洞中分別對縮尺比1∶300的A、B、C、D類四類風場進行模擬［17］。因C、D類風場具有較大空氣動力學粗糙長度z0，在數(shù)值模擬中需采用本文方法，即以具有最大空氣動力學粗糙長度的D類風場為例說明其應(yīng)用。模型風場的z0=6.756×10－3m(對應(yīng)實際風場z0=2.0 m)，取參考高度zR=1.0 m，對應(yīng)參考風速為UR=5 m/s。

2.1 CFD模型與計算方案

縮尺計算域的大小為12×1.5 m(Lx×Lz)，其中x表示流向方向，z表示豎向方向，見圖1。采用非均勻結(jié)構(gòu)化網(wǎng)格方案劃分計算域網(wǎng)格，地面附近沿高度方向網(wǎng)格最小尺寸為0.01 m(由此給出zP=0.005 m)，流向方向網(wǎng)格尺寸為0.03 m，網(wǎng)格總數(shù)為24 000。速度－壓力耦合方式為SIMPLEC，對流項的離散格式為QUICK，壓力插值方式為PRESTO，所有計算方案的數(shù)值實現(xiàn)方法可參考Ansys 13的幫助文件。

2.2 邊界條件

圖1為入流面、出流面、頂面及地面位置處所用邊界條件，該邊界條件名稱源于Ansys 13。各邊界條件數(shù)學意義見表1。平均風速及湍動能來流邊界條件源于同濟大學TJ－2風洞模擬縮尺比1∶300的D類風場。湍動能k=0.5σu2，其中σu為流向脈動風速均方根值。計算中設(shè)局部地區(qū)湍流處于平衡狀態(tài)，即ε=Ck?U/?z，湍動能沿高度變化［12］。由于模型風場的z0=6.756 ×10－3m，即使據(jù)Ks=20z0［9］，該模型風場的物理粗糙高度Ks也遠大于zP=0.005 m。因此，需在標準壁面函數(shù)后添加附加項δB解決此矛盾。上述方法采用Ansys 13的用戶自定義壁面函數(shù)實現(xiàn)。為強調(diào)附加項δB的重要性，用戶輸入值Ks=0.005 m，該值為地面附近網(wǎng)格能允許的最大值。

圖1 O－XZ坐標系、縮尺計算域及邊界條件Fig.1 Coordinate system O－XZ，scaled computational domain and boundary conditions

表1 具有最大空氣動力學粗糙長度風場的邊界條件及數(shù)學意義Tab.1 Boundary conditions and their corresponding mathematical implications in simulating the neutral ABL with the largest aerodynamic roughness length

2.3 湍流模型常數(shù)及附加項δB對地面附近流動影響

在標準k－ε湍流模型中共5個湍流模型常數(shù)，即Cμ，σk，σε，C1ε，C2ε;其中Cμ為最重要的模型常數(shù)，其余均與Cμ有直接或間接關(guān)系。各模型常數(shù)的計算方法見文獻［3，12］。據(jù)C=u/k，Cμ由來流邊界條件決定。由于標準k－ε湍流模型中Cμ為常數(shù)，而k沿高度變化，因此須尋找一高度對應(yīng)的k值獲得合適的Cμ及其它湍流模型常數(shù)。圖2的計算結(jié)果表明，模型常數(shù)Cμ=0.23，σε=0.77，σk=0.59，C1ε=1.44，C2ε= 1.92及δB=5.5時，來流邊界條件可獲得較好模擬。由于在計算域流向方向來流邊界條件均保持較好，用于比較的計算結(jié)果來源于出流邊界位置處。

為說明附加項的重要性，模型常數(shù)不變、δB=0.0及3.0時計算結(jié)果見圖2。由圖2看出，δB較小時與來流邊界條件相比，地面附近的平均風速偏大;湍動能被低估。平均風速偏大原因為δB較小時壁面第一層網(wǎng)格中心P點處UP偏大。換言之，文獻［15］的修正壁面函數(shù)不能提供較小平均風速UP，從而導致計算域來流邊界上的來流邊界條件及計算域內(nèi)建筑結(jié)構(gòu)前的來流出現(xiàn)差異。δB=0時即為修正的壁面函數(shù)式(4)的計算結(jié)果。此時γ=1?？梢?，為模擬D類風場或為模擬具有較大空氣動力學粗糙長度的風場，須采用附加項。

圖2 附加項δB對地面附近流動影響Fig.2 Effects of the additional term δB on the flow near the land surface

31 ∶50縮尺比TTU低矮建筑模型繞流模擬

3.1 TTU建筑模型與縮尺比風場簡介

TTU建筑模型常用于對風洞試驗、數(shù)值模擬兩種方法有效性驗證［18］。在TJ－2風洞縮尺比1∶50的B類風場條件下，羅攀［17］對TTU建筑模型在不同風向角條件下的繞流進行研究。TTU模型、風向角定義及橫截面測點位置見圖3。本文對60°及90°風向角條件下繞流進行數(shù)值模擬。縮尺比1∶50風場對應(yīng)的空氣動力學粗糙長度為z0=0.001 58 m(實際風場z0=0.08 m)。取參考高度HR=1.0 m、參考風速UR=6 m/s。相同地貌不同縮尺比的模型風場對應(yīng)的空氣動力學粗糙長度不同。對縮尺比1∶300的B類風場無需附加項，但對縮尺比為1∶50時，則需在標準壁面函數(shù)中含附加項。

圖3 TTU模型、風向角定義及橫截面上測點位置Fig.3 TTU model，definition of the wind azimuth and positions of the measurement points on the edge of the cross section

3.2 CFD模型、計算方案及邊界條件

文獻［19］詳細考察了計算域?qū)?shù)值計算結(jié)果影響。本文計算域在流向、橫向及垂直方向尺寸分別為－6H≤x≤20H，－10H≤y≤10H，0≤z≤7H，見圖4 (a)，其中H為縮尺TTU模型高度。地面附近第一層網(wǎng)格尺度0.008 m，對應(yīng)zP=0.004 m，見圖4(b)。計算域的入流面、出流面、頂面、側(cè)面、地面及TTU模型表面邊界條件(名稱源于Ansys 13)見表2。標準k－ε湍流模型常數(shù)為Cμ=0.13，σε=1.02，σk=0.78，C1ε= 1.44及C2ε=1.92。為較好模擬該模型風場，須在標準壁面函數(shù)中用附加項δB=4.5。用戶輸入量Ks=0.004 m，為zP=0.004 m允許最大值。

3.3 數(shù)值模擬、風洞試驗與現(xiàn)場實測結(jié)果對比

定義風壓系數(shù)為

式中:Pi為測點平均壓力;PH為TTU模型屋頂高度處靜壓;UH為未受擾動的來流平均風速;ρ=1.225 kg/m3為空氣密度。

圖5為模型橫截面各測點平均風壓系數(shù)的數(shù)值模擬結(jié)果、風洞試驗結(jié)果及現(xiàn)場實測結(jié)果比較。圖5(a)為在迎風面屋檐處，采用二階QUICK格式標準k－ε湍流模型(SKE)計算結(jié)果。由于降低了湍動能在屋檐處的堆積［20］，MMK k－ε湍流模型(MMK)結(jié)果較滿意。需要說明的是，SKE及MMK兩湍流模型的控制方程及模型常數(shù)均相同。對圖5(b)中繞流更復雜工況，因數(shù)值粘性降低，采用三階MUSCL格式的MMK k－ε湍流模型結(jié)果亦較滿意。由圖5看出，數(shù)值模擬結(jié)果與風洞試驗結(jié)果更接近，可能在數(shù)值計算中，由于采用附加項，風場模擬的得好，故風洞試驗結(jié)果更接近數(shù)值計算結(jié)果。

表2 縮尺比1∶50的TTU模型繞流計算域邊界條件及數(shù)學意義Tab.2 Boundary conditions and their corresponding mathematical implications in simulating the flow around the TTU model in a neutral wind field with the scale of 1∶50

圖4 風向角為90°時計算域、O－XYZ坐標系、TTU模型在計算域內(nèi)位置及TTU模型表面的網(wǎng)格Fig.4 Scaled computational domain for the wind azimuth at 90 degree，O－XYZ coordinate system and location of the TTU model and the mesh on the TTU model surface

圖5 橫截面各測點平均風壓系數(shù)的數(shù)值模擬結(jié)果、風洞試驗結(jié)果和現(xiàn)場實測結(jié)果比較Fig.5 Comparisons of the mean pressure coefficients at the measurement points among the numerical，wind tunnel and real－site results:wind azimuth at 90 degree and wind azimuth at 60 degree

圖6 TTU模型表面風壓系數(shù)的數(shù)值模擬結(jié)果Fig.6 Numerical results of the pressure coefficients on the TTU surface:wind azimuth at 90 degree and wind azimuth at 60 degree

需要指出的是，本文數(shù)值風場參考高度HR=1.0 m對應(yīng)參考風速為UR=6 m/s，該風速相對較低。據(jù)計算流體動力學中關(guān)于無量綱尺度z*=ρu*zP/μ定義，若忽略雷諾數(shù)影響，在z*保持不變條件下，低風速可采用較大網(wǎng)格尺度，從而降低計算量。對較復雜的建筑群數(shù)值模擬，可考慮低風速從而降低計算量［21］。不同工況條件下TTU模型表面風壓系數(shù)分布見圖6。由圖6可見，對90°風向角工況，在迎風面屋檐及墻角處，由于分離流、風壓系數(shù)變化較劇烈，且風壓系數(shù)沿橫截面對稱分布;對60°風向角工況，由于出現(xiàn)三角翼渦，加劇了漩渦脫落的復雜性，其風壓系數(shù)變化更復雜。

4 結(jié)論

(1)本文對基于標準k－ε湍流模型數(shù)值模擬大氣邊界層中的壁面函數(shù)問題進行研究。商用CFD軟件中通常要求Ks≤zP;而對模擬具有較大空氣動力學粗糙長度的大氣邊界層，Ks通常大于zP。該矛盾為出現(xiàn)壁面函數(shù)問題的根本原因，從而導致計算域地面附近的來流邊界條件及建筑結(jié)構(gòu)前的來流出現(xiàn)差異。通過分析，用附加項δB模型化地表上面建筑結(jié)構(gòu)等形成的粗糙元因大小不一、錯亂分布對地表附近空氣流動產(chǎn)生的附加影響。

(2)通過模擬縮尺比1:300具有較大空氣動力學粗糙長度二維中性大氣邊界層及縮尺比1:50的TTU低矮建筑模型在三維中性大氣邊界層內(nèi)繞流，對附加項的有效性及使用場合進行評估、說明。結(jié)果表明，由于能在計算域的流向方向保持來流邊界條件，附加項能有效解決壁面函數(shù)問題。本文參考風速較低，TTU建筑模型計算表明不影響計算結(jié)果可靠性;且對較復雜建筑群數(shù)值模擬，低風速可降低計算量。

［1］Jensen M.The model law for phenomena in the natural wind［J］.Igenioren(International Edition)，1958，2:121－128.

［2］Richards P J，Younis B A.Comments on“prediction of the wind－generated pressure distribution around buildings”by E H Mathews［J］.Journal of Wind Engineering and Industrial Aerodynamics，1990，34(1):107－110.

［3］Richards P J，Hoxey R P.Appropriate boundary conditions for computational wind engineering models using thek－εmodel［J］.JournalofWindEngineeringandIndustrial Aerodynamics，1993，46/47:145－153.

［4］Richards P J，Norris S E.Appropriate boundary conditions forcomputational wind engineering models revisited［J］.Journal of Wind Engineering and Industrial Aerodynamics，2011，99: 257－266.

［5］Launder B E，Spalding D B.The numerical computation of turbulent flows［J］.Computer Methods in Applied Mechanics and Engineering，1974，3(2):269－289.

［6］Suga K，Craft T J，Iacovides H.An analytical wall－function for turbulent flows and heat transfer over rough walls［J］. International Journal of Heat and Fluid Flow，2006，27(5): 852－866.

［7］Craft T J，Gerasimov A V，Iacovides H，et al.Progress in the generalization of wall－function treatments［J］.International Journal of Heat and Fluid Flow，2002，23(2):148－160.

［8］Miles S，Westbury P.Practical tools for wind engineering in the built environment［J］.QNET－CFD Network Newsletter，2003，21:11－14.

［9］Blocken B，Stathopoulos T，Carmeliet J.CFD simulation of the atmospheric boundary layer:wall function problems［J］. Atmospheric Environment，2007，41(2):238－252.

［10］ESDU.Engineering sciences data item 85020［S］.EngineeringSciencesDataUnit，251－259RegentStreet，London，1985.

［11］Hargreaves D M，Wright N G.On the use of the k－εmodel in commercial CFD software to model the neutral atmospheric boundary layer［J］.Journal of Wind Engineering and Industrial Aerodynamics，2007，95(5):355－369.

［12］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):88－95.

［13］Yang W，Quan Y，Jin X Y，et al.Influences of equilibrium atmosphere boundary layer and turbulence parameter on wind loads of low－rise buildings［J］.Journal of Wind Engineering and Industrial Aerodynamics，2008，96(10/11):2080－2092.

［14］普朗特，等.著.郭永懷，等.譯.流體力學概論［M］.北京:科學出版社，1981.

［15］Cebeci T，Bradshaw P.Momentum transfer in boundary layers［M］.New York:Hemisphere Publishing Corporation，1977.

［16］GB50009－2001，中華人民共和國國家標準:建筑結(jié)構(gòu)荷載規(guī)范［S］.

［17］羅攀.基于標準模型的風洞試驗研究［D］.上海:同濟大學，2004.

［18］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):227－234.

［19］方平治，顧明，談建國，等.阻塞率對表面風壓系數(shù)影響的數(shù)值模擬［J］.建筑科學與工程學報，2013，30(3): 101－106.

FANG Ping－zhi，GU Ming，TAN Jian－guo，et al.Numerical simulation of effect of blockage ratio on facade pressure coefficient［J］.JournalofArchitectureandCivil Engineering，2013，30(3):101－106.

［20］Murakami S，Mochida A，Kondo K，et al.Development of new k－ε model for flow and pressure fields around bluff body［J］.Journal of Wind Engineering and Industrial Aerodynamics，1997，66/67:169－182.

［21］方平治，史軍，王強，等.上海陸家嘴區(qū)域建筑群風環(huán)境數(shù)值模擬研究［J］.建筑結(jié)構(gòu)學報，2013，34(9):104－111.

FANG Ping－zhi，SHI Jun，WANG Qiang，et al.Numerical study on the wind environment among tall buildings in Lujiazui zone［J］.Journal of Building Strucutures，2013，34(9): 104－111.

Method to solve the wall function problem in simulation of atmospheric boundary layer

FANG Ping－zhi1，GU Ming2，TAN Jian－guo3，HAN Zhi－h(huán)ui3
(1.Shanghai Typhoon Institute of China Meteorological Administration，Shanghai 200030，China;
2.State Key Laboratory for Disaster Reduction in Civil Engineering，Tongji University，Shanghai 200092，China; 3.Shanghai Institute of Meteorological Science，Shanghai 200030，China)

Wall function is preferred to model the low Reynolds－number flow near wall based on the Reynoldsaveraged Navier－Stokes turbulent models.Then wall function problem in simulating the atmospheric boundary layer based on the standard turbulent model was investigated.An extra term which considers the extra effects induced by non－uniform and irregular distribution of rough elements such as various structures on the land surface was proposed and appended to the widely accepted standard wall function.The effectiveness and application situation of the proposed term were demonstrated by simulating a neutral wind field with the scale of 1:300 featured by larger aerodynamic roughness length. The flow around the TTU model in a neutral wind field with the scale of 1:50 was then simulated.The extra term is shown to be necessary to solve the wall function problem，i.e.，the preservation inlet flow boundary conditions in computation domain in both cases

atmospheric boundary layer;standard wall function;wall function problem;physical roughness height; aerodynamic roughness length

O353.1

A

10.13465/j.cnki.jvs.2015.02.015

2013－11－08修改稿收到日期:2014－02－14

方平治男，博士，副研究員，1974年生

顧明男，博士，教授，博士生導師，1957年生

郵箱:migngu@#edu.cn