楊 易，季長慧，張之遠，胡曉兵,2

(1. 華南理工大學亞熱帶建筑科學國家重點實驗室，廣東，廣州 510641；2. 華中科技大學土木與水利工程學院，湖北，武漢 430074)

計算風工程(computational wind engineering，CWE)中的大渦模擬(large eddy simulation，LES)因其能提供高分辨率的流場時空信息，近年來在風氣候、地形效應(yīng)、污染擴散、空氣動力學效應(yīng)等研究中得到越來越廣泛的應(yīng)用。如Tominaga等[1]運用大渦模擬研究建筑物近場的污染擴散；胡偉成、楊慶山等[2]運用大渦模擬研究三維山丘地形風場及湍流特性；俞怡恬、黃生洪等[3]運用大渦模擬研究龍卷風沖擊高層建筑氣動力效應(yīng)。

大氣邊界層中風場以隨機湍流形式運動，為準確模擬建筑結(jié)構(gòu)的脈動風效應(yīng)，因此需要在數(shù)值模擬中首先對大氣邊界層的湍流風場準確再現(xiàn)。在LES研究中，如何準確、高效地生成一個能夠同時滿足大氣邊界層湍流場的無源性、空間相關(guān)性和脈動風速功率譜(下文簡稱功率譜)特性的湍流入口，并在整個流域內(nèi)保持平衡——即大氣邊界層湍流風場的準確模擬與數(shù)值重構(gòu)，是一個基礎(chǔ)性科學問題，也是當前CWE研究熱點和難題之一[4]。LES數(shù)值模擬建筑結(jié)構(gòu)繞流結(jié)果，在很大程度上取決于入流湍流及平衡態(tài)邊界層的模擬的準確性[5]。

隨機湍流生成方法(random field generation，RFG)作為一種“主動”LES入流湍流生成方法，近年來一直備受國內(nèi)外學者關(guān)注。RFG方法最早在2001年由Smirnov等[6]在Kraichnan[7]的基礎(chǔ)上建議，而后2010年Huang等[8]基于RFG理論上提出離散與合成湍流生成方法(discretizing and synthesizing random flow generation，DSRFG)。2013年Castro和 Paz[9]通過改進DSRFG方法，提出改進的離散與合成湍流生成方法(modified discretizing and synthesizing random flow generation，MDSRFG)；2015年Aboshosha等[10]在DSRFG方法基礎(chǔ)上提出一致性離散湍流合成方法(consistent discrete random flow generation,CDRFG)。作者課題組在2018年基于RFG方法，參考CDRFG方法提出了一種新的窄帶湍流合成法—NSRFG(narrowband synthesis random flow generation)法[11]，基于該方法開展了標準地貌風場大渦模擬研究[12]，并應(yīng)用在超高層建筑風荷載及風振響應(yīng)模擬和參數(shù)優(yōu)化研究中[13]。Feng等[14]利用NSRFG方法研究偏轉(zhuǎn)風對超高層建筑的風效應(yīng)問題，黃銘楓等[15]基于WRF(weather research and forecasting)和NSRFG LES進行臺風多尺度耦合數(shù)值模擬研究。

研究顯示，雖然新建議的NSRFG方法相對此前方法具有計算效率和精度相對更高的優(yōu)點，但其中用來調(diào)整流場空間相關(guān)性的關(guān)鍵參數(shù)—調(diào)諧因子γj的取值規(guī)律，仍需深入探討和進一步明確。筆者團隊在此前研究基礎(chǔ)上，提出改進的INSRFG(improved NSRFG)方法。

1 NSRFG方法

NSRFG方法建議的湍流風場數(shù)學模型[11]如下式所示：式中：ui為 i 方向的速度( i =1,2,3分別代表順風向、橫風向和垂直方向)；x 為 x、y、z方向空間坐標向量； x ?j,n=xj/Lj,n， j=1,2,3 分別代表x、y、z方向； Lj,n=Uav/(γjCjfn)，cj為衰減系數(shù)，γj為調(diào)諧因子； Su,i(fn)為 i 方向的頻率為fn帶寬為Δf的功率譜；φn為在區(qū)間[0,2π]服從均勻分布的隨機相位參數(shù)，即φn～U(0,2π)。kj,n={k1,n,k2,n,k3,n}為服從空間單位球面上均勻分布的隨機向量，以滿足空間分布要求，即 |kn|=1 。此外參數(shù)kj,n需保持湍流場的無源性條件，即滿足下式：

根據(jù)式(2)及 | kn|=1 ，參數(shù)kj,n={k1,n,k2,n,k3,n}進而可以轉(zhuǎn)換成服從空間圓弧曲線上均勻分布的向量。因此，可以利用空間圓曲線的參數(shù)方程計算參數(shù)kj,n的值，如下式所示：

式中：qi,n=pi,n/Li,n，i =1,2,3；θ ～U(0,2π)；系數(shù)An和Bn取值如下：

在NSRFG方法[11]中，參考Aboshosha等[10]的做法，針對我國《建筑結(jié)構(gòu)荷載規(guī)范》中的C類風場，直接對風場衰減系數(shù)和調(diào)諧因子賦值，其模擬得到的湍流場能較好滿足空間相關(guān)性以及功率譜特性條件。

2 NSRFG方法的改進

脈動風速的空間相關(guān)性是檢驗所模擬的湍流風場準確性的關(guān)鍵指標之一。基于湍流合成法的RFG系列方法中，對空間相關(guān)性調(diào)整參數(shù)Li,n的處理是這類方法的關(guān)鍵。

Aboshosha等[10]針對DSRFG方法[8]不能嚴格滿足目標功率譜譜和空間相關(guān)性條件的不足，分別在入流功率譜和脈動風速空間相關(guān)性兩方面進行了修正，進而提出了CDRFG方法。但CDRFG方法中，空間相關(guān)性調(diào)諧因子γj在y和z方向上取值與x方向相同，使得y與z方向的空間相關(guān)性條件得不到保證?？臻g相關(guān)性的不完全滿足顯然將會影響湍流風場數(shù)值模擬結(jié)果的準確性。

NSRFG方法相對CDRFG方法有了明顯改進，但不同地貌類別的湍流風場中衰減系數(shù)cj和調(diào)諧因子γj的適當取值，仍需要大量試算來確定。針對這一問題，本文在CDRFG方法和NSRFG方法的基礎(chǔ)上，同時考慮風速分量u、v、w在三維空間上不同的脈動特性，提出u、v、w三個方向上空間一致相干性的改進建議，并通過理論推導給出數(shù)學表達式，以改進NSRFG方法中衰減系數(shù)和適當調(diào)諧因子取值的難題。

2.1 空間相關(guān)性調(diào)整參數(shù)

基于NSRFG方法，通過式(1)求出空間上任意兩點時程的空間相關(guān)函數(shù)，如下式所示：

式中， Cf=cifd/Uav。其中：d為空間兩點的距離；ci為衰減系數(shù)； Uav為平均速度；γj為調(diào)諧因子。從式(5)可以看出，衰減系數(shù)ci和調(diào)諧因子γj是決定 Lj,n空間相關(guān)性調(diào)整參數(shù)的兩個關(guān)鍵參數(shù)。

下文將通過理論推導，通過無量綱長度尺度β給出調(diào)諧因子γj的數(shù)學表達式。無量綱長度尺度β=ciD/Lu(z) ，其中 Lu(z)是湍流順風向湍流積分尺度，D是調(diào)整相關(guān)性的特征距離[10]，ci是衰減系數(shù)[8]。

2.1.1 調(diào)諧因子γj與空間相關(guān)性R的關(guān)系

給定某一空間相關(guān)性R=0.63的相干函數(shù)擬合過程見圖1，擬合的準則為相干曲線包圍的面積(即寬頻帶中的空間相關(guān)性)相等，即：

圖1 NSRFG方法中給定空間相關(guān)性值R =0.63的相干函數(shù)擬合過程Fig.1 Fitting process of coherence function for a spatial correlation value R =0.63 in NSRFG method

通過以上擬合，可以獲得調(diào)諧因子γj和空間相關(guān)性 Ru1u2的散點圖；然后，對γj和 Ru1u2的散點進行曲線擬合(圖2)，進而可得到 Ru1u2與γj的數(shù)學關(guān)系。

圖2 調(diào)諧因子γ j 與空間相關(guān)性 Ru1u2的關(guān)系Fig.2 Relationship between tuning factor γ j and spatial correlationRu1u2

2.1.2 無量綱長度尺度β與空間相關(guān)性R的關(guān)系

湍流風場順風向兩點脈動風速u1、u2的空間相關(guān)性可通過兩點的協(xié)方差 σu1u2和兩者脈動風速的均方差 σu1、 σu2的比值獲得。順風向脈動風速功率譜采用von Karman譜，空間相關(guān)性 Ru1u2與無量綱長度尺度β可由下式得到：

式中：

同理可得，在y、z方向上，空間相關(guān)性Ru1u2與無量綱長度尺度β如下所示：

通過以上分析，將獲得三個方向上的空間相關(guān)性 Ru1u2與無量綱長度尺度β的分布關(guān)系，如圖3所示；則進一步可推出調(diào)諧因子γj與無量綱長度尺度β的關(guān)系，如圖4所示。

圖3 無量綱尺度β與空間相關(guān)性 Ru1u2的關(guān)系Fig.3 Relationship between dimensionless scaleβ and spatial correlationRu1u2

圖4 調(diào)諧因子γ j 與無量綱尺度β的關(guān)系Fig.4 Relationship between tuning factor γ j anddimensionless scaleβ

2.1.3 調(diào)諧因子γj與無量綱長度尺度β的關(guān)系

以上分別得到空間相關(guān)性 Ru1u2與無量綱長度尺度β以及調(diào)諧因子γj的空間分布關(guān)系。進一步，通過擬合得到關(guān)于調(diào)諧因子γj與無量綱長度尺度β的函數(shù)關(guān)系，如式(10)所示：

式(10)即為本文建議的改進的NSRFG方法(INSRFG方法)中的三維空間相關(guān)性調(diào)整參數(shù)γj的數(shù)學模型。INSRFG方法考慮了大氣邊界層中湍流積分尺度隨高度變化的特點，保證了所生成的脈動風速的空間相關(guān)性特性。

2.2 數(shù)值驗證

INSRFG方法生成的脈動風速場的功率譜和空間相關(guān)性是否與理論邊界層風場模型相符，需要采用數(shù)值模擬方法進行檢驗。為了研究的連續(xù)性和可對比性，數(shù)值驗證采用的數(shù)值模型參數(shù)設(shè)置與文獻[12]保持一致。衰減系數(shù)Cj參考文獻[8]和文獻[12]取值，具體在下文中說明；特征距離D為與建筑高度H相關(guān)的參數(shù)，一般取為0.5H～1.0H，本文根據(jù)作者計算經(jīng)驗取為1.0H(對于本算例為0.6 m)[10]。

2.2.1 脈動風速功率譜

數(shù)值模型中，采用von Karman譜描述脈動風速的統(tǒng)計特性，則空間三維方向上的脈動風速譜采用如下公式[10]：

圖5為3種不同湍流生成方法仿真得到的脈動風速功率譜比較。結(jié)果顯示，NSRFG方法和INSRFG方法二者結(jié)果相近，相對CDRFG方法，所生成的脈動風速功率譜與目標譜更加吻合；而CDRFG方法的數(shù)值模擬結(jié)果幅值波動較大，誤差也相對較大。

圖5 不同LES入口湍流生成方法得到的三維脈動風速時程功率譜比較Fig.5 Comparisons of spectra of 3D fluctuating velocities with different LES inflow turbulence generation methods

2.2.2 空間相關(guān)性

將INSRFG方法模擬得到的湍流場的三維脈動風速空間相關(guān)性與理論目標曲線的比較，如圖6所示?？臻g相關(guān)性目標函數(shù)采用Hemon和Santi[17]提出的公式：

圖6 三維脈動速度分量的空間相關(guān)性比較Fig.6 Comparisons of spatial correlations of 3D fluctuating velocity components

結(jié)果表明，所模擬的湍流場三個方向上的空間相關(guān)性均與目標函數(shù)吻合較好，即INSRFG方法能較好滿足大氣邊界層三維脈動風場對空間相關(guān)性的要求。

3 湍流大氣邊界層風場的LES模擬

為檢驗這一新方法在計算風工程中的適用性，基于ANSYS Fluent軟件平臺采用UDF(userdefined function)功能進行二次開發(fā)，進行與我國《建筑結(jié)構(gòu)荷載規(guī)范》(GB 50009?2012)[18]中定義的四類標準地貌湍流大氣邊界層風場的LES INSRFG模擬。數(shù)值建模中流域尺寸、邊界條件以及網(wǎng)格劃分(圖7)與文獻[12]保持一致(H、B、D為CAARC標準高層建筑1/300縮尺模型尺寸，即H=0.6096 m，B=0.1016 m，D=0.1524 m)。

圖7 數(shù)值模型和網(wǎng)格劃分Fig.7 Numerical model and mesh discretization

INSRFG模型中衰減系數(shù)對于A類地貌取值為c1=10、c2=15、c3=15；其余三種地貌取值均為c1=10、c2=12、c3=12，調(diào)諧因子γj根據(jù)本文建議的式(10)取值。數(shù)值求解格式和風場參數(shù)見表1。

表1 LES數(shù)值求解格式和風場參數(shù)列表Table1 Schemes and parameters in LES simulation

3.1 平均風速和湍流強度剖面

在計算域約x=L/3處(L為計算域順風向長度)設(shè)置一排豎向監(jiān)測點，通過大渦模擬得到各監(jiān)測點的脈動風速時程，進而求出湍流風場的平均風速和湍流強度剖面。圖8為NSRFG方法與INSRFG方法得到的大氣邊界層計算域約x=L/3位置處的湍流風場特性模擬結(jié)果。圖中：Zg為梯度風高度；Ug為梯度風速；Iu為順風向湍流強度；α為風剖面指數(shù)。

圖8 NSRFG方法與INSRFG方法模擬得到的四類標準地貌湍流風場平均風速和湍流強度剖面Fig.8 Mean wind speed and turbulence intensity profiles of four standard wind terrains by NSRFG and INSRFG methods

圖8湍流風場模擬結(jié)果顯示，對于平均速度剖面，模擬得到的計算域約x=L/3位置處的平均風速與理論目標速度剖面二者幾乎完全保持一致，即速度剖面的保持性良好；而對于湍流強度剖面，約x=L/3位置處的湍流強度相對理論目標值略小，這是由于大渦模擬的濾波作用，導致從入口到約x=L/3建筑模型位置的湍流強度發(fā)生衰減所致?；赗FG方法模擬得到的湍流的自保持性問題，有待進一步的研究和解決[4]。相對而言，與NSRFG方法模擬得到的4類標準地貌湍流風場結(jié)果相比，在湍流強度方面，INSRFG方法的模擬精度整體上略高，與規(guī)范值[18]更加吻合。

表2為圖7中H高度處相應(yīng)的順風向脈動速度均方根與目標值的比較。結(jié)果顯示，模擬的中國規(guī)范4類標準地貌下順風向脈動速度分量與目標值的相對誤差在5%以內(nèi)，表明INSRFG方法模擬的湍流風場結(jié)果與規(guī)范值吻合較好。

表2 4類標準地貌順風向脈動風速均方根與目標值比較Table2 Comparison of RMS values of simulated along-wind velocities with the target values for four standard wind terrains

數(shù)值模擬結(jié)果表明，整體上，基于LES INSRFG方法模擬得到的湍流風場，在不經(jīng)任何修正的情況下，可較好滿足大氣邊界層湍流風場特性要求和基本滿足平衡態(tài)要求。

4 結(jié)論

本文基于NSRFG方法，通過對無量綱長度尺度β、調(diào)諧因子γj及空間相關(guān)性R進行分析，建議了調(diào)諧因子γj與無量綱長度尺度β的數(shù)學模型，進而建議了改進的INSRFG方法，總結(jié)如下：

(1)INSRFG方法的優(yōu)點體現(xiàn)在既彌補了CDRFG方法只考慮順風向空間相關(guān)性的局限性，又改進了NSRFG方法衰減系數(shù)和調(diào)諧因子難以準確賦值的不足。

(2)通過數(shù)值模擬仿真和比較分析，驗證了INSRFG方法生成的大氣邊界層湍流風場能較好滿足脈動風速功率譜、空間相關(guān)性等脈動風場基本特性；4類標準地貌LES模擬顯示，能基本滿足計算風工程中對平衡態(tài)大氣邊界層風場的基本要求。

(3)INSRFG湍流合成方法具有數(shù)字表示式簡潔、計算效率高、參數(shù)取值明確等優(yōu)勢，計算精度相對較高，是一種進行建筑結(jié)構(gòu)LES模擬研究具有很好前景的通用LES入口湍流生成方法，可為土木工程數(shù)值風洞技術(shù)的推廣和應(yīng)用提供有價值的參考。