大氣邊界層大渦模擬入口湍流生成方法綜述

周 桐，楊慶山，閆渤文，Pham Van Phuc，王京學

(1.北京交通大學土木建筑工程學院，北京 100044；2.結(jié)構(gòu)風工程與城市風環(huán)境北京市重點實驗室，北京 100044；3.重慶大學土木工程學院，重慶 400044；4.清水建設(shè)技術(shù)研究所，日本，東京 135-8530)

現(xiàn)場實測、理論分析和風洞試驗是結(jié)構(gòu)風工程的三種傳統(tǒng)研究方法。隨著計算流體動力學方法(Computational Fluid dynamics，簡稱CFD)的不斷改進和計算機軟硬件水平的飛速發(fā)展，CFD成為結(jié)構(gòu)風工程研究的重要工具，并逐漸發(fā)展成為一個重要分支，被稱為計算風工程(Computational Wind Engineering，簡稱 CWE)。雷諾平均法和大渦模擬是計算風工程的兩種主要求解方法，其中大渦模擬是一種對大于網(wǎng)格尺度的湍流脈動進行直接模擬，而對小于網(wǎng)格尺度的湍流脈動采用亞格子模型來描述的非穩(wěn)態(tài)模擬方法。大渦模擬作為一種兼顧求解精度和求解效率的數(shù)值模擬方法，正逐漸成為結(jié)構(gòu)風效應(yīng)數(shù)值模擬研究的主要選擇[1－3]。城市大氣邊界層研究涉及的尺度通常可劃分為城市尺度(宏觀尺度)、小區(qū)尺度(中尺度)和街區(qū)尺度(微尺度)，而不同尺度上的大氣運動呈現(xiàn)出不同的空氣流動特點[4]，因此需要發(fā)展適用于不同尺度特征的大氣邊界層數(shù)值模擬入口生成方法。例如，針對小區(qū)尺度的大氣微環(huán)境研究，由于同時受到較大尺度氣象因素和小尺度湍流脈動的影響，其入口來流特性應(yīng)具有顯著的多尺度特性，需采用多尺度耦合模擬[5－7]。而本文主要關(guān)注的是近地大氣邊界層內(nèi)的結(jié)構(gòu)風效應(yīng)，需要模擬的主要是小尺度湍流脈動對結(jié)構(gòu)表面的脈動風壓分布和結(jié)構(gòu)周圍的繞流特性[8]，因此真實重現(xiàn)目標風場內(nèi)的小尺度湍流脈動特性是采用大渦模擬準確評估結(jié)構(gòu)風效應(yīng)的重要前提。然而，如何合理定義大渦模擬的入口邊界條件一直都是計算風工程研究領(lǐng)域關(guān)注的重難點問題[9]。理想的結(jié)構(gòu)風效應(yīng)大渦模擬入口邊界條件應(yīng)滿足以下四點要求：1)正確重構(gòu)表征大氣邊界層風場的主要統(tǒng)計量，包括平均風速剖面、湍流強度、脈動風速功率譜和湍流積分尺度等；2)具有與真實風場相接近的流場結(jié)構(gòu)；3)生成的脈動風速具有隨機性；4)適用于并行計算環(huán)境，并且計算成本較低。

預(yù)前模擬法和人工合成法是目前主要的兩類大渦模擬入口湍流生成方法。然而，由于缺乏從結(jié)構(gòu)風工程研究的角度對不同方法的特點以及適用性進行系統(tǒng)探討，研究人員難以合理選擇和正確應(yīng)用大渦模擬入口湍流生成方法開展結(jié)構(gòu)風工程問題的研究。本文首先綜述了不同方法的基本原理，梳理其在結(jié)構(gòu)風工程領(lǐng)域的應(yīng)用和發(fā)展，對比分析不同方法的特點及適用性，為基于大渦模擬方法開展結(jié)構(gòu)風工程問題的研究提供了重要的理論指導(dǎo)。

1 預(yù)前模擬法

在預(yù)前模擬法中，求解的計算域包含兩部分：主計算域和輔助計算域。首先在輔助計算域生成滿足大氣邊界層風場特性的穩(wěn)定來流，然后將其作為主計算域的入口邊界條件。由于輔助計算域內(nèi)生成湍流的具體方式不同，預(yù)前模擬法主要有兩種實現(xiàn)形式，分別是被動模擬法和“回收-變換”法。

1.1 被動模擬法

在建筑結(jié)構(gòu)風洞實驗中，尖劈、粗糙元等裝置被廣泛應(yīng)用于模擬不同地貌類別的大氣邊界層風場。采用風洞試驗的被動模擬技術(shù)，通過建立尖劈、粗糙元的數(shù)值風洞模型，并在入口邊界賦予與風洞試驗相接近的均勻來流，進而一定程度上再現(xiàn)風洞試驗中模擬的大氣邊界層風場，如圖1所示。 被動模擬法的基本思想易于理解，并且技術(shù)難度較低，因而被研究人員應(yīng)用于生成結(jié)構(gòu)風效應(yīng)大渦模擬的目標風場[10－12]。然而，這種方法在實際應(yīng)用時主要存在以下兩點局限性：1)工作效率低。不同的目標風場對應(yīng)于不同的風場布置，因此在模擬不同目標風場條件下的結(jié)構(gòu)風效應(yīng)前，需要調(diào)整數(shù)值風洞中尖劈、粗糙元的幾何形狀和幾何位置，并重新劃分網(wǎng)格。2)計算成本高。計算域內(nèi)的網(wǎng)格分辨率很大程度上決定了采用被動模擬法生成的大氣邊界層風場是否與風洞試驗結(jié)果相一致。

圖1 被動模擬法數(shù)值模型Fig.1 Numerical model of the passive simulation method

為解決上述不足，Enoki和Ishihara[13]提出了在計算域中建立“隱性粗糙元”的方法，即通過在動量方程中添加阻力源項，以近似地模擬粗糙元對來流產(chǎn)生的拖曳力效應(yīng)。

式中：Cf=CD,ui/(1-γ0)2，CD,ui為阻力系數(shù)，建議取值為0.4[14]；γ0=V0/Vgrid，γ0為體積占有率，V0為粗糙元體積，Vgrid為包含粗糙元的網(wǎng)格體積；l0=V0/A0，l0為特征長度，A0為單個粗糙元的迎風面積。

在此基礎(chǔ)上，Liu等[15]將其延伸應(yīng)用于光滑三維小山和二維山脊周圍湍流場的大渦模擬。這種方法避免了在輔助計算域內(nèi)建立數(shù)值粗糙元模型，當目標風場對應(yīng)的粗糙元的幾何尺寸、排列形式發(fā)生改變時，可通過修改相關(guān)程序代碼快速地進行調(diào)整，有效地提高了工作效率。然而，該方法難以描述尖劈對于風場的影響。

1.2 “回收-變換”法

基于湍流邊界層的相似變換理論，“回收-變換”法將提取得到的湍流邊界層下游的平均統(tǒng)計量和脈動統(tǒng)計量進行調(diào)節(jié)后作為計算域的入口邊界條件。這種方法最初應(yīng)用于低雷諾數(shù)平板湍流邊界層[16]，并不適用于土木工程結(jié)構(gòu)所處的高雷諾數(shù)粗面湍流邊界層。為適應(yīng)結(jié)構(gòu)數(shù)值抗風研究的需求，在Lund等[17]方法的基礎(chǔ)上，Nozawa和Tamura[18]在入口平面和循環(huán)站間建立數(shù)值粗糙元并合理地調(diào)整速度縮尺函數(shù)，實現(xiàn)了粗糙壁面湍流邊界層的模擬，然后將其作為低矮建筑結(jié)構(gòu)風效應(yīng)大渦模擬的入口邊界條件。通過引入邊界層厚度沿順流向不發(fā)生變化的假設(shè)，Kataoka[19]簡化了Lund等[17]提出的擬周期邊界條件，將下游循環(huán)站的脈動風速分量縮尺后與入口平面的平均風速分量疊加：

式中：u、v、w分別為順風向、橫風向和豎向的瞬時風速；分別為順風向、橫風向和豎向的平均風速；下標inlt表示入口平面；下標recy表示循環(huán)站；φ(θ)為權(quán)函數(shù)。

為進一步模擬具有較高湍流強度的大氣邊界層風場，基于 Kataoka[19]提出的簡化方法，朱偉亮和楊慶山[20]通過建立數(shù)值粗糙元和擾流桿模型增大了近地面和邊界層高處的湍流強度，如圖2所示。

圖2 “回收-變換”法數(shù)值模型Fig.2 Numerical model of the recycling-rescaling method

由于主計算域和輔助計算域的粗糙度存在顯著差異，來流的湍流強度在主計算域內(nèi)會發(fā)生一定程度的衰減，尤其是邊界層高處。王婷婷和楊慶山[21]通過在一定高度以上添加滿足正態(tài)分布的隨機數(shù)，提高了邊界層高處的湍流強度，但其風速時程呈現(xiàn)出明顯的周期性，因此生成的流場特性與真實大氣邊界層風場不符。Aboshosha等[22]采用隨機傅里葉模式生成分形表面，以模擬任意粗糙壁面的湍流邊界層。Li等[23]提出了基于脈動速度比例系數(shù)λ(z,t)的主動控制方法，以生成滿足目標湍流強度特性的邊界層。

式中：Iu,obj(y,z,t)為目標湍流強度；Iu,re(y,z,t)為瞬時湍流強度；為橫向平均值；為循環(huán)站處的平均速度；N為模擬時間步數(shù)；Δt為時間步長。

提高生成目標湍流的效率并避免虛假低頻脈動的產(chǎn)生對于“回收-變換”法在結(jié)構(gòu)風工程領(lǐng)域的應(yīng)用具有重要意義。Liu和Pletcher[24]通過動態(tài)變換循環(huán)站的位置，使循環(huán)站處的湍流充分發(fā)展，有效減少了生成目標湍流的時間；Jewkes等[25]采用鏡像方法有效地避免了入口邊界和循環(huán)站間的誤差反饋和累積，縮短了計算域長度，從而降低了計算成本；Stevens等[26]通過實時將輔助計算域出口區(qū)域內(nèi)的數(shù)據(jù)傳輸?shù)街饔嬎阌?，提高了求解效率。Spille-Kohoff和 Kaltenbach[27]采用在計算域入口邊界與循環(huán)站間添加源項的方法，對流場進行有效的反饋控制，降低了人工速度場長時間相關(guān)形成虛假低頻脈動的可能性；Morgan等[28]將動態(tài)變換的橫向反射應(yīng)用于循環(huán)站，從而避免虛假低頻脈動的產(chǎn)生。

2 人工合成法

人工合成法是基于嚴格的數(shù)理推導(dǎo)構(gòu)造滿足大氣邊界層風場特性(平均風速剖面、湍流強度、脈動風速功率譜等)的入口邊界條件。根據(jù)生成湍流脈動的具體方式不同，人工合成法主要分為四類，分別是譜合成法、本征正交分解重構(gòu)法、數(shù)字濾波法、渦方法。

2.1 譜合成法

基于對頻譜空間和波譜空間的隨機傅里葉變換，譜合成法可分為傳統(tǒng)的諧波合成法和傅里葉合成法兩類。

2.1.1 諧波合成法

諧波合成法是利用脈動風速的目標功率譜密度和空間相關(guān)性在頻域內(nèi)構(gòu)造含有高斯隨機系數(shù)項的三角級數(shù)，脈動風速樣本可以表達為一系列的正弦和余弦函數(shù)。

式中：Δω為頻率增量，Δω=(ωmax-ωmin)/N1；N1為采樣頻率點數(shù)，理論上N1→∞；Hjm(ωml)為目標功率譜S(ωml)的Cholesky分解矩陣H(ωml)中的元素；φml為均勻分布于[0,2π]的隨機相位角。

Rice[29]首先提出了諧波合成法的基本思想，但是其僅適用于一維單變量平穩(wěn)高斯隨機過程，無法模擬大氣邊界層風場中不同位置處的脈動風速時程。經(jīng)過Shinozuka等[30－32]的改進，將諧波合成法發(fā)展到多變量、非平穩(wěn)高斯過程的模擬，進而適用于生成真實的隨機脈動風場。基于頻率雙索引的概念，Deodatis[33]提高了諧波合成法的求解精度，模擬了各態(tài)歷經(jīng)的多變量平穩(wěn)高斯隨機過程。李正農(nóng)等[34]采用非均勻圓頻率間隔的方法有效地解決了以往諧波合成法模擬的周期性問題。上述方法改善了模擬精度，但是卻降低了模擬效率。采用諧波合成法生成大渦模擬脈動入口邊界時，通常需要模擬入口平面多點的脈動風速時程，而隨著模擬點數(shù)的增加會導(dǎo)致其計算效率的顯著降低[35]，因此提高其求解效率至關(guān)重要。由于諧波合成法求解過程涉及對復(fù)雜矩陣進行大量的 Cholesky分解以及三角級數(shù)疊加，因此Cholesky分解與三角級數(shù)的疊加是制約其計算效率的主要因素。Yang[36－37]采用快速傅里葉變換技術(shù)簡化了三角級數(shù)的疊加，進而大幅度提升了諧波合成法的模擬效率。孫瑛等[38]基于本征特征正交分解技術(shù)，建立能夠近似代表實際風場頻譜特性的少量階時間主坐標的互譜密度矩陣，在保證計算精度的前提下提高了計算效率。Ding等[39]將脈動風速的互譜密度矩陣簡化為實對稱矩陣，利用三次Lagrange插值方法減少了Cholesky分解的次數(shù)，進而縮短了計算時長。羅俊杰和韓大建[40]引入并改進基于矩陣分塊的“遞歸”算法以提高Cholesky分解效率，并采用矩陣乘法取代了效率低下的疊加算法。李春祥和劉晨哲[41]利用徑向基神經(jīng)網(wǎng)絡(luò)對Cholesky分解后的互譜密度矩陣進行插值，在保證諧波合成法模擬精度的同時，顯著地提高了模擬效率。Huang等[42]采用相位分解，將演化功率譜密度矩陣/功率譜密度矩陣轉(zhuǎn)換為實模矩陣，進而提高了Cholesky分解效率?；陔p索引頻率和譜分解矩陣的特點，祝志文和黃炎[43]提出在雙索引頻率軸向采取均勻分布插值，譜分解矩陣的頻率軸向采取前密后疏分段插值能夠大幅提高隨機脈動風場的模擬效率。對于Cholesky分解矩陣的非對角線元素，采用三次Lagrange插值方法的模擬值在某些插值區(qū)間存在較大波動。因此，陶天友和王浩[44]提出了基于Hermite插值方法的諧波合成法，以兼顧模擬精度和計算效率。

諧波合成法具有嚴格的數(shù)學理論基礎(chǔ)，并且模擬的脈動風速功率譜與脈動風速的自/互相關(guān)函數(shù)均能夠與目標值吻合較好。然而，由于其生成的脈動風速時程無法嚴格滿足連續(xù)性條件，因此會在計算域內(nèi)產(chǎn)生虛假的壓力脈動并可能導(dǎo)致計算過程的發(fā)散。同時，由于該方法無法獨立生成不同空間點的脈動風速時程，且生成的風速時程數(shù)據(jù)需要預(yù)先儲存，因此不適用于并行計算。

2.1.2 傅里葉合成法

傅里葉合成法是基于三維能量譜在波數(shù)域內(nèi)生成與目標功率譜相一致的、空間相關(guān)的湍流場。

Kraichnan和 Robert[45]首先提出了傅里葉合成法的基本思想，但其僅適用于生成均勻、各向同性的湍流場。以Kraichnan和Robert[45]提出的算法為基礎(chǔ)，Smirnov等[46]引入了湍流長度尺度和時間尺度，并通過比例和正交變換生成了非均勻各向異性的湍流場。然而，采用上述方法僅能生成滿足高斯譜形式的湍流場，并不適用于結(jié)構(gòu)抗風研究。然而，通過改變的分布形式可以模擬具有不同能量譜特性的脈動速度場。當各向同性地分布在半徑為k0的空間球面或平面圓表面時，對應(yīng)脈動速度場的能量譜分別為利用E1(k)和E3(k)僅在k0處不為零的特性，以構(gòu)造滿足任意能量譜的脈動速度場：

式中：km為波長值；M為能量譜離散的區(qū)間樣本數(shù)?；谏鲜隼碚?，Huang等[47]提出了 DSRFG(discretizing and synthesizing random flow generation)方法，其脈動速度場的表達形式為：

在結(jié)構(gòu)風工程領(lǐng)域，DSRFG方法主要具有以下幾點優(yōu)勢：1)嚴格滿足連續(xù)性條件；2)生成的湍流場基本滿足脈動風速功率譜；3)能夠獨立生成不同空間點的脈動風速時程，適用于并行計算。然而，DSRFG方法中的Ls取值并未考慮其與相干函數(shù)的內(nèi)在聯(lián)系。基于數(shù)學推導(dǎo)，Castro和 Paz[48]指出采用DSRFG方法生成的隨機脈動風場的湍流強度與離散譜的區(qū)間樣本數(shù)相關(guān)。當離散譜的區(qū)間樣本數(shù)達到一定數(shù)量時，湍流強度的模擬值能夠較好地吻合目標值。在 DSRFG方法的基礎(chǔ)上，Castro和Paz[48]通過時間尺度參數(shù)τ0的引入來考慮脈動風速的時間相關(guān)性。

脈動風的空間相關(guān)性對于準確評估柔性結(jié)構(gòu)(高層建筑、高聳結(jié)構(gòu)、大跨度橋梁等)的風致耦合作用十分重要，Aboshosha等[49]在 DSRFG方法的基礎(chǔ)上，提出了 CDRFG(consistent discretizing random flow generation)方法。這種方法修正了脈動風速譜能量在頻率上的分布，并且通過建立空間尺度參數(shù)與頻率的關(guān)系，比較合理地考慮了脈動風的豎向相關(guān)性。

式中：C為相干函數(shù)衰減系數(shù)；D=0.5~1.0h；h為建筑高度。類似于DSRFG方法，為使生成的脈動速度場滿足連續(xù)性條件，須滿足下列方程：

由于DSRFG方法、CDRFG方法需要耗費較多的計算資源并且脈動風速的表達形式相對復(fù)雜，因此，Yu等[50]在DSRFG方法、CDRFG方法的基礎(chǔ)上，通過窄帶過程的模擬和疊加，并考慮三個方向上脈動風的空間相關(guān)性提出了表達形式更為簡潔、嚴格滿足連續(xù)性條件和脈動風速功率譜的 NSRFG(narrow-band synthesizing random flow generation)方法。

式中：ui(i=1,2,3)為順風向、橫風向和豎向的脈動風速；N3為功率譜離散數(shù)目；Si(fn)(i=1,2,3)分別為順風向、橫風向和豎向的脈動風速功率譜；fn=( 2n-1)Δf2，Δf為帶寬；φn~U(0,2π)；為使生成的脈動速度場滿足連續(xù)性條件，須滿足下列方程：

從空間解析幾何的角度分析，上述方程為一空間圓曲線，因此ki,n(i=1,2,3)可由式(25)求解：

2.2 本征正交分解重構(gòu)法

本征正交分解法是一種隨機脈動風場的高效描述方法[51－52]，僅需通過少量階與時間主坐標和空間本征模態(tài)的級數(shù)組合來重構(gòu)隨機脈動風場的主要信息[53－55]。

式中：a(n)(t)為時間主坐標；ψ(x)為本征模態(tài)；n為所選本征模態(tài)階數(shù)。

基于熱線風速儀或粒子圖像測速儀采集得到平面內(nèi)各點的脈動風速時程的本征正交分解，結(jié)合隨機模擬技術(shù)，理論上可以構(gòu)造出滿足風洞流場特性的入口邊界條件。熱線風速儀具有較高時間分辨率，但由于熱線風速儀同步測量的空間點數(shù)量有限，因此其空間分辨率較低。對此，Druault等[56]采用線性隨機估計填補了測量位置間的空缺信息，進而生成了大渦模擬的入口邊界條件。然而，線性隨機估計的準確性很大程度上取決于空間點的距離。相比于熱線風速儀，粒子圖像測速儀具有較高空間分辨率。因此，Perret等[57－58]通過對三維粒子圖像測速儀的測量數(shù)據(jù)進行本征正交分解，進而生成了空間發(fā)展平板混合層大渦模擬的入口邊界條件。然而，由于其采樣頻率相對較低，因此需要引入隨機時間序列來提高時間分辨率。

采用本征正交分解重構(gòu)法生成的入口湍流需要經(jīng)過較長的發(fā)展距離才能演化成為真實湍流，Johansson和Andersson[59]基于本征正交分解法準確定位具有最多能量的模態(tài)，并通過納維-斯托克斯方程的伽遼金映射，將低能量、小尺度的模態(tài)疊加到高能量、大尺度的模態(tài)上以平衡耗散水平，實現(xiàn)了速度分量間更加真實的能量分布。

2.3 數(shù)字濾波法

數(shù)字濾波法的基本思想是利用數(shù)字濾波器將離散隨機序列轉(zhuǎn)變?yōu)闈M足指定時空相關(guān)性的隨機序列。

式中：N5為濾波器的長度()，Δx為網(wǎng)格尺寸，l為湍流積分尺度；rm為滿足零均值、單位方差的隨機序列；bj為濾波系數(shù)，其表達式為。

隨機序列的時間相關(guān)性同樣可采用指數(shù)函數(shù)來定義：

式中：C為常數(shù)(π/4)；T為拉格朗日時間尺度，T=l/Ucon，Ucon為平均對流速度。

基于數(shù)字濾波法生成的瞬時速度場為：

式中：Ui為平均速度；aij為雷諾應(yīng)力Rij的Cholesky分解[17]。

Klein等[60]首先提出了數(shù)字濾波法的基本思想，指出采用該方法能夠再現(xiàn)目標湍流場的二階統(tǒng)計量和自相關(guān)性，但是其計算效率受限于網(wǎng)格形式。Kempf等[61]基于并行算法程序的改進，使其在復(fù)雜的非結(jié)構(gòu)化網(wǎng)格形式下仍然具有較高的計算效率。Xie和 Castro[62]提出了簡化的濾波方法，在每個時間步內(nèi)僅需要進行二維濾波，因此進一步地提高了該方法的計算效率。同時，該方法是易于并行的，但其生成的湍流場無法嚴格滿足零散度條件，因此會在計算域內(nèi)引起顯著的壓力脈動，進而無法準確評估結(jié)構(gòu)上的瞬時風荷載(極值)。為避免虛假的壓力脈動，Kim等[63]改進了不可壓縮流動求解器中的速度-耦合程序，Daniels等[64]將其應(yīng)用于評估高層建筑物的極值風荷載和表面脈動壓力。采用數(shù)字濾波法生成的湍流場通常無法與納維-斯托克斯方程相容，導(dǎo)致其入口的湍流特性無法在計算域內(nèi)得到良好的保持，然而，結(jié)構(gòu)風工程領(lǐng)域主要關(guān)注建筑物附近的來流湍流特性。為滿足結(jié)構(gòu)數(shù)值抗風研究的需要，Lamberti等[65]在Kim等[63]方法的基礎(chǔ)上，以建筑物區(qū)域的湍流特性為目標，采用梯度優(yōu)化算法來確定入口湍流生成程序中的輸入?yún)?shù)。該方法的基本步驟可簡要概括為以下兩步：1)基于貝塞爾曲線擬合目標湍流特征統(tǒng)計量(如雷諾應(yīng)力、湍流積分尺度等)剖面；2)基于標量化技術(shù)確定目標函數(shù)(以雷諾應(yīng)力為例)。

式中：yp為貝塞爾曲線控制點p的高度；γu、γv和γw為權(quán)重因子；下標i、b和exp分別表示入口平面、建筑物和風洞實驗的雷諾應(yīng)力值。

2.4 渦方法

基于渦量的拉格朗日描述，二維渦方法(VM)在入口平面構(gòu)造隨機分布的渦以產(chǎn)生湍流脈動。

式中：N6為入口平面的渦數(shù)量；Γi為渦環(huán)通量，為入口平面的面積；k為湍流動能；x為每個渦的二維坐標；xi為每個渦中心的二維坐標；z為順流向的單位矢量；σi為渦的特征尺度。

基于湍流混合長假設(shè)，Mathey等[66]將局部渦的特征尺度與入口平面的湍流動能和湍流耗散率相聯(lián)系()，并采用簡化的線性運動學模型來描述順流向的湍流脈動。在結(jié)構(gòu)風工程領(lǐng)域，湍動能和耗散率通常采用下述公式進行計算：

在二維渦方法的基礎(chǔ)上，Jarrin等[67－68]提出了SEM(synthetic eddy method)，采用形函數(shù)來定義具有時空相關(guān)性的三維渦旋相干結(jié)構(gòu)，進而生成脈動速度場。

式中：Ui為順流向的平均速度；N7為入口平面的渦數(shù)量；aij為雷諾應(yīng)力Rij的Cholesky分解；為獨立隨機變量，在(-1,1)上服從均勻分布；xk為渦的空間位置坐標；fσ(x-xk)為xk位置處渦的速度分布函數(shù)。

式中：f為形函數(shù)，常見的包括三角函數(shù)、高斯函數(shù)等；σ為渦的特征尺度，由湍動能、耗散率、網(wǎng)格尺度共同決定。

式中：κ=0.41；Δ=m ax(Δx,Δy,Δz)；δ為邊界層厚度。

基于泰勒湍流凍結(jié)假設(shè)以定義脈動速度場的時間相關(guān)性，在每個迭代時間步內(nèi)的第k個渦的空間位置為：

式中，Uc為平均速度。

基于上述合成渦方法能夠生成滿足部分大氣邊界層風場特性的湍流場，如平均風速剖面、湍流強度、湍流積分尺度、時空相關(guān)性。同時，其各點風速時程的生成過程獨立，因此適用于并行計算。然而，由于其生成渦的特征尺度單一，因此不符合脈動風速功率譜特性。針對這一問題，Luo等[69]提出了 MSSEM (multi-scale synthetic eddy method)，其基本思想是合成具有不同頻譜能量的多尺度渦以構(gòu)造滿足任意脈動風速功率譜特性的湍流場。

式中：U(x)為順流向的平均速度；M為渦的特征尺度的數(shù)量；Nm為m階尺度渦的數(shù)量，其定義式為為a的m階分量，其中a為雷諾應(yīng)力的 Cholesky分解，其定義式為為隨機3×3矩陣，在(-1,1)上服從均勻分布；為位置處的速度分布函數(shù)。

基于不同尺度渦之間的獨立性假設(shè)，方向i上的脈動速度功率譜可表示為：

式中：n為頻率；為與m階尺度渦相關(guān)的脈動速度功率譜。

式中，F(xiàn)i為形函數(shù)f的傅里葉變換。

相應(yīng)地，無量綱化的脈動速度功率譜為：

結(jié)合式(42)~式(45)可得：

Luo等[70]采用MSSEM生成了滿足大氣邊界層風場特性的脈動入口邊界條件，進而運用大渦模擬比較準確地評估了高層建筑上的風壓分布和風荷載特性。

3 結(jié)論

生成滿足大氣邊界層風場特性的入口邊界條件是合理運用大渦模擬開展結(jié)構(gòu)風工程研究的關(guān)鍵瓶頸問題。本文對預(yù)前模擬法和人工合成法這兩類大渦模擬入口湍流生成方法的基本原理和發(fā)展歷程進行綜述，并從結(jié)構(gòu)風工程研究的角度對不同方法的特點及適用性進行對比分析，結(jié)果表明：

(1)采用預(yù)前模擬法生成的入口邊界條件具有相對真實的流場結(jié)構(gòu)與時空相關(guān)性。同時，由于其滿足流體運動的基本方程(連續(xù)性方程、納維-斯托克斯方程)，因此入口湍流特性在計算域內(nèi)能夠得到良好的保持，并且數(shù)值計算易于收斂。然而，預(yù)前模擬法無法直接定義目標湍流特性(湍流強度、湍流積分尺度和脈動風速功率譜)，調(diào)試過程相對復(fù)雜。當網(wǎng)格質(zhì)量和網(wǎng)格分辨率足夠高時，運用被動模擬法能夠比較真實地再現(xiàn)大氣邊界層風洞的流場特性。然而，被動模擬法需要在主要計算域外增加額外的預(yù)前模擬計算域，網(wǎng)格量增加較多，從而降低了工作效率、提高了計算成本，另一方面，要求主要計算域邊界形狀規(guī)則，這兩點在一定程度上限制了其在結(jié)構(gòu)風工程領(lǐng)域研究中的應(yīng)用?；凇盎厥?變換”法生成的入口邊界條件滿足大氣邊界層的平均風速剖面，其局限性主要體現(xiàn)在邊界層的湍流強度較低，與近地脈動風場的高湍流度特征不符。為適應(yīng)結(jié)構(gòu)數(shù)值抗風研究的需求，通常采用前置粗糙元法來增大近地面的湍流度，但邊界層高處的湍流度仍然低于目標值。在“回收-變換”法的基礎(chǔ)上，結(jié)合基于脈動速度比例系數(shù)的主動控制方法，在一定程度上能夠生成滿足目標湍流強度特性的大氣邊界層風場。

(2)采用人工合成法生成的入口邊界條件滿足大氣邊界層風場特性與時空相關(guān)性。同時，該方法的計算效率較高、適用范圍較廣。由于人工合成湍流無法滿足納維-斯托克斯方程，需要在計算域入口邊界與計算模型之間預(yù)留較長的發(fā)展距離來形成目標湍流，而在此過程中，湍流特性可能因數(shù)值模型、網(wǎng)格尺寸、時間步長等因素發(fā)生改變?；谥C波合成法能夠生成滿足任意目標譜特性和空間相關(guān)性的脈動入口邊界，但由于其無法與連續(xù)性方程相容，會在計算域內(nèi)引起顯著的壓力脈動，進而無法準確評估結(jié)構(gòu)的瞬時風荷載。譜合成法中的DSRFG、CDRFG和NSRFG方法嚴格滿足連續(xù)性條件，入口湍流的空間相關(guān)性可通過空間尺度參數(shù)調(diào)整，其中NSRFG方法的并行計算效率相對較高。本征正交分解重構(gòu)法的應(yīng)用很大程度上取決于是否具有類似大氣邊界層湍流風場的實驗數(shù)據(jù)，因此其在結(jié)構(gòu)數(shù)值抗風研究中存在較大的局限性。數(shù)字濾波法生成的入口邊界無法滿足納維-斯托克斯方程，因此入口湍流特性無法在計算域內(nèi)得到良好的保持。相比于譜合成法、本征正交分解重構(gòu)法和數(shù)字濾波法，渦合成法生成的入口湍流能更快地發(fā)展成為真實湍流。其中，MSSEM 生成的湍流場能較好地滿足脈動風速功率譜，進而更適用于結(jié)構(gòu)風工程研究。結(jié)合大渦模擬入口湍流生成方法在結(jié)構(gòu)風工程領(lǐng)域的研究現(xiàn)狀，對后續(xù)研究作出以下幾點展望：

(1)在“回收-變換”法的基礎(chǔ)上，進一步發(fā)展多目標湍流特性(湍流強度、湍流積分尺度等)的主動控制技術(shù)，提高“回收-變換”法的調(diào)試效率。

(2)人工合成法生成的入口湍流特性可能無法在計算域內(nèi)較好地保持，結(jié)合相關(guān)優(yōu)化算法對人工合成湍流程序中的輸入?yún)?shù)進行適當調(diào)整，使來流在感興趣的研究區(qū)域具有目標湍流特性。

(3)在滿足大氣邊界層主要統(tǒng)計特性(平均風速剖面、湍流強度和脈動風速功率譜)的前提下，運用人工合成法生成比較合理的大氣邊界層湍流結(jié)構(gòu)，即旋渦尺度隨高度的增加而增加，并且隨地面粗糙度的增加而減小。

(4)來流湍流特性會顯著影響結(jié)構(gòu)表面的風荷載特性，并且不同結(jié)構(gòu)(高層建筑、低矮建筑、大跨度屋蓋結(jié)構(gòu)等)風荷載特性具有顯著差異，因此應(yīng)進一步對比各種入口湍流生成方法在不同結(jié)構(gòu)風效應(yīng)大渦模擬研究中的適用性。