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

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

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

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

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

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

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

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

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

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

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

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

1 標準壁面函數及附加項

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

對光滑壁面，對數律為

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

對粗糙壁面，對數律為

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

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

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

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

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

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

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

2.1 CFD模型與計算方案

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

2.2 邊界條件

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

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

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

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

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

為說明附加項的重要性，模型常數不變、δB=0.0及3.0時計算結果見圖2。由圖2看出，δB較小時與來流邊界條件相比，地面附近的平均風速偏大;湍動能被低估。平均風速偏大原因為δB較小時壁面第一層網格中心P點處UP偏大。換言之，文獻［15］的修正壁面函數不能提供較小平均風速UP，從而導致計算域來流邊界上的來流邊界條件及計算域內建筑結構前的來流出現差異。δB=0時即為修正的壁面函數式(4)的計算結果。此時γ=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建筑模型常用于對風洞試驗、數值模擬兩種方法有效性驗證［18］。在TJ－2風洞縮尺比1∶50的B類風場條件下，羅攀［17］對TTU建筑模型在不同風向角條件下的繞流進行研究。TTU模型、風向角定義及橫截面測點位置見圖3。本文對60°及90°風向角條件下繞流進行數值模擬。縮尺比1∶50風場對應的空氣動力學粗糙長度為z0=0.001 58 m(實際風場z0=0.08 m)。取參考高度HR=1.0 m、參考風速UR=6 m/s。相同地貌不同縮尺比的模型風場對應的空氣動力學粗糙長度不同。對縮尺比1∶300的B類風場無需附加項，但對縮尺比為1∶50時，則需在標準壁面函數中含附加項。

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

3.3 數值模擬、風洞試驗與現場實測結果對比

定義風壓系數為

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

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

表2 縮尺比1∶50的TTU模型繞流計算域邊界條件及數學意義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模型在計算域內位置及TTU模型表面的網格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 橫截面各測點平均風壓系數的數值模擬結果、風洞試驗結果和現場實測結果比較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模型表面風壓系數的數值模擬結果Fig.6 Numerical results of the pressure coefficients on the TTU surface:wind azimuth at 90 degree and wind azimuth at 60 degree

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

4 結論

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

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

［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，中華人民共和國國家標準:建筑結構荷載規范［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］方平治，顧明，談建國，等.阻塞率對表面風壓系數影響的數值模擬［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］方平治，史軍，王強，等.上海陸家嘴區域建筑群風環境數值模擬研究［J］.建筑結構學報，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－hui3
(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@tongji.edu.cn