周 桐，楊慶山，閆渤文，Pham Van Phuc，王京學(xué)

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

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

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

1 預(yù)前模擬法

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

1.1 被動(dòng)模擬法

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

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

為解決上述不足，Enoki和Ishihara[13]提出了在計(jì)算域中建立“隱性粗糙元”的方法，即通過(guò)在動(dòng)量方程中添加阻力源項(xiàng)，以近似地模擬粗糙元對(duì)來(lái)流產(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為特征長(zhǎng)度，A0為單個(gè)粗糙元的迎風(fēng)面積。

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

1.2 “回收-變換”法

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

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

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

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

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

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

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

2 人工合成法

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

2.1 譜合成法

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

2.1.1 諧波合成法

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

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

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

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

2.1.2 傅里葉合成法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.3 數(shù)字濾波法

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

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

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

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

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

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

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

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

2.4 渦方法

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

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

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

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

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

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

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

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

式中，Uc為平均速度。

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

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

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

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

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

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

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

Luo等[70]采用MSSEM生成了滿足大氣邊界層風(fēng)場(chǎng)特性的脈動(dòng)入口邊界條件，進(jìn)而運(yùn)用大渦模擬比較準(zhǔn)確地評(píng)估了高層建筑上的風(fēng)壓分布和風(fēng)荷載特性。

3 結(jié)論

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

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

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

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

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

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

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