邱順冬 程 旭

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

0 引言

由于設計方法及新建筑材料和技術的創(chuàng)新發(fā)展所驅動，新一代高層建筑正在變得越來越苗條、高柔，這就使得高層建筑對諸如風等橫向載荷變得更加敏感。此外，為了控制地震產生的慣性力，通常需要降低建筑物的重量，這又進一步有助于高層建筑風荷載及其導致的高層建筑位移的增加。因此，風荷載越來越成為高層建筑橫向結構設計主要考慮的因素。目前，高層建筑風荷載的評估主要有風洞試驗[1]和CFD(Computational Fluid Dynamics)模擬兩種方法。雖然風洞試驗便于改變和重復試驗條件，然而成本高、受試驗條件影響及其所測試驗數(shù)據(jù)有限等缺點嚴重限制了其發(fā)展。另外一方面，憑借其在硬件成本和運行耗能等方面的優(yōu)勢以及流場全貌顯示等特點，關于使用CFD模擬來評估高層建筑風荷載的研究變得越來越廣泛。

CFD模擬的湍流模型主要可以分為RANS(Reynold-Averaged Navier-Stokes)和LES(Large Eddy Simulation)兩種。RANS在平均風特性上表現(xiàn)良好，但其很難預測脈動風成分。而LES方法通過濾波函數(shù)將大尺度的渦和小尺度的渦分離開，大尺度的渦直接模擬，小尺度的渦用模型來封閉，這就使得LES從某種程度上可以直接計算得到脈動風成分，故其也正在成為越來越流行的CFD計算方法。而要使用LES進行風荷載的準確評估，入口邊界條件則成為一個十分重要的影響因素。目前針對入口邊界條件的LES(Large Eddy Simulation)數(shù)值模擬方法主要可以分為兩大類：人工合成湍流法和預前模擬法。

人工合成湍流法主要是將產生的波動速度疊加目標平均速度。波動速度的產生常用的主要有以下幾種方法：譜合成、數(shù)字濾波法及擬序結構疊加。為此，國內外很多學者在此基礎上進行了改進，如Smirnov et al[2]，Batten et al[3]，Huang et al[4]，Mathey et al[5]。這種方法的優(yōu)點是產生非均勻湍流,適合復雜的入口。但其很難滿足連續(xù)性方程，這就導致入流風特性進入計算域內部以后會發(fā)生變化，比如在建筑物附近不能達到目的湍流度等，這就給實際應用帶來了困難。

預前模擬是在被研究建筑之前通過放置尖劈及粗糙元等鈍體，被動模擬出目標風場。此方法最大的優(yōu)勢在于其模擬出的風場滿足多點相關性及連續(xù)性方程，近似于真實的湍流，對模擬鈍體繞流具有十分重要的意義。然而為了模擬出充分發(fā)展的湍流，且獲得足夠的統(tǒng)計數(shù)據(jù)，往往需要大量的計算網(wǎng)格和計算時間，這也是預前模擬法在工程無法普遍應用的重要原因。由此，為了減少存儲量，Spalart[6]提出了周期邊界，Lund et al[7]對其進行了改進。Kataoka[8]對Lund et al[7]的方法進行了簡化，提出了擬周期邊界條件。

本文主要針對擬周期邊界條件提出了用于粗糙元設計方法，并且通過大量的LES模擬，驗證了公式的正確性，為此方法的應用提供了參考。

1 數(shù)值建模及求解設置

1.1 控制方程

大渦模擬(LES)湍流模型控制方程可以表達為：

(1)

(2)

1.2 計算域及邊界條件

本模擬的計算域如圖1所示，其中計算域尺寸為2 200 m(順風向)×550 m(橫風向)×550 m(豎向)。此外，計算域底部布置間距為100 m的粗糙塊，順風和橫風的尺寸均為15 m，而粗糙塊高度的設置與風場類別有關。

入口采用Kataoka[8]提出的擬周期邊界，其中循環(huán)面與入口的距離為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分別為入口和循環(huán)面上提取的數(shù)值。φ(θ)的表達式如式(4)所示，是由Kataoka[8]提出的阻尼函數(shù)。

(4)

其中，θ=y/Ly,Ly為大氣邊界層高度。而權重數(shù)φ(θ)≤1.0，作用是穩(wěn)定循環(huán)后的速度波動，達到目標值。

2 結果與討論

2.1 LES模擬結果分析

本次模擬的參數(shù)主要選取了平均風剖面、湍流度剖面，這兩者是通過對大氣邊界層風場特征參數(shù)和規(guī)范對比確定的，模擬的目標分別為中國規(guī)范指數(shù)律風剖面和日本規(guī)范湍流度剖面。如圖2，圖3所示，平均風速剖面和湍流度剖面的模擬結果與所對應的規(guī)范值基本吻合，證明了本文采用擬周期方法的適用性。

2.2 粗糙元的設計

在以往基于LES的預前模擬法中并沒有給予直接選取粗糙元大小的方法，所以本文在參考Wooding et al[10]風洞試驗中所采用的公式的基礎上，對風場模擬中的粗糙元進行設計，并提出相應的公式。

Wooding et al[10]在文中推導出的應用于風洞試驗的經驗公式如式(5)所示:

(5)

其中，δ為湍流剪切層或邊界層厚度；κ為卡曼常數(shù)；k為粗糙元的高度；λ為粗糙度；φ為關于k和s的函數(shù)，φ≈(k/s)β，β為一個常數(shù)；U為z=δ處的風速；u*為總剪切速度；C為一個Wooding et al在論文中確定的常數(shù)，C≈-2.05。

根據(jù)u*=(τ/ρ)1/2，τ為總剪切應力;ρ為流體密度。從阻力系數(shù)的計算公式τ=0.5ρU2Cf可以得到總剪切速度的公式：

(6)

其中在Wooding et al的論文中，s=d，s為粗糙元平行于流場的水平維度，d為粗糙元垂直于來流方向的水平寬度，在本文中s與d相同，D為粗糙元之間的間距。

由此可以將式(5)推導成：

(7)

其中，λ=Ar/S，Ar為粗糙元的正面面積，S為每個粗糙元所占有的區(qū)域，S=D2，D為粗糙元之間的間距；β為一個常數(shù)；κ為卡曼常數(shù)，根據(jù)Irwin[11]得到κ≈0.35，β≈0.4。Wooding et al和Irwin[11]都提到了，公式的應用范圍就是70<δ/kλφ<2 000。

將各個常數(shù)代入到式(7)中，可得到式(8)：

(8)

通過數(shù)值模擬的結果，本文計算出Cf的數(shù)值，采用多項式擬合出Cf和α的關系，Cf=0.3α2-0.03α+0.01。依據(jù)規(guī)范，不同的地面類型選取不同的α的數(shù)值，其中A,B,C,D四類地面所對應的α分別為0.12,0.15,0.22,0.30。根據(jù)Hagishima et al[12]得出的結論，λ<0.2時，Cf隨著λ的變大而上升。地面粗糙度指數(shù)α和Cf在A,B,C,D四類地表情況下，取值都是逐漸變大，所以本文在選取λ時也要在0～0.2的區(qū)間逐漸變大。根據(jù)王婷婷[13]的LES模擬結果，計算出其模型中式(8)中包含的各個參數(shù)，代入式(8)中能夠發(fā)現(xiàn)比較好的吻合。各參數(shù)的數(shù)值如表1所示。

表1 各類地表LES模擬得到的參數(shù)

本文還通過計算的算例檢驗這個公式對于粗糙元設計的指導性。本文選取B類地表，對應的地面粗糙度指數(shù)α為0.15，大氣邊界層的厚度為350 m。根據(jù)擬合曲線算出的Cf=0.011，對應的λ取0.1，可通過式(6)算出k1.4/d0.4=20.86。流場的正面尺寸為550 m×550 m，粗糙元的間距取為50 m，根據(jù)λ=Ar/S，得到kd=250 m2。計算得到d≈14 m，k≈19 m。根據(jù)這一數(shù)據(jù)可以把粗糙元的尺寸和布置情況大致選定，算例所采用的粗糙元尺寸順流向和橫流向的尺寸為15 m，高度取20 m。

對比結果發(fā)現(xiàn)，應用公式計算得到的粗糙塊的尺寸和間距，LES模擬得到的大氣邊界層的風場特征參數(shù)與規(guī)范B類地表能夠較好的吻合，如圖2，圖3所示。這一實例也證實了式(8)在基于LES的預前模擬法粗糙元設計中具有可操作性。

3 結語

本文對B類地表的大氣邊界層風場進行了LES模擬，通過與規(guī)范對比驗證模擬結果的正確性，并提出了相應的粗糙塊設計方法。

采用擬周期方法，通過布置粗糙塊和諧波合成方法來改變計算域入口風特性設置，使得模擬得到的建筑物附近的風速剖面和湍流度剖面與規(guī)范值基本符合。進一步通過大量的模擬，提出了基于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.