羅嘯宇， 談?wù)z， 吳曉波， 劉震卿， 涂元剛

(1. 廣東電網(wǎng)有限責(zé)任公司電力科學(xué)研究院， 廣東 廣州 510062；2. 華中科技大學(xué) 土木工程與力學(xué)學(xué)院， 湖北 武漢 430074)

隨著觀測技術(shù)水平的提高，人們逐漸識別出一大部分氣象災(zāi)難的罪魁禍?zhǔn)资窍聯(lián)舯┝?。下?lián)舯┝魇菤饬餮杆傧陆禌_擊地面后，形成的從沖擊中心向四周沿地表方向輻散的一種極具破壞性的強風(fēng)。下?lián)舯┝髟诤娇蘸教旌徒ㄖY(jié)構(gòu)兩個領(lǐng)域?qū)θ祟惿a(chǎn)生活有重大威脅，嚴重情況可能導(dǎo)致災(zāi)難。

目前，學(xué)者主要采用三種研究方法對其展開研究，分別是基于歷史記錄的氣象數(shù)據(jù)的風(fēng)場實測研究、基于下?lián)舯┝鳟a(chǎn)生和發(fā)展原理的氣象理論分析、以及依靠專用的計算機進行實驗室仿真模擬。在實驗方面，Mason等[1]開發(fā)并采用了下?lián)舯┝髂M器，重點研究了下?lián)舯┝黠L(fēng)時程特性。仿真模擬方面，Aboshosha等[2]采用大渦模型研究了下?lián)舯┝黠L(fēng)場特性，著重研究了下?lián)舯┝鞯姆e分尺度、風(fēng)譜特性等。翟偉廉等[3,4]都采用了二維計算域進行數(shù)值模擬，在計算域的邊界條件限制了流場發(fā)展。陳波等[5]采用三維計算域進行數(shù)值模擬，但是只重點研究了風(fēng)場對于平屋面的風(fēng)荷載特性。Aboshosha等[2]進行實際尺寸的數(shù)值模擬但是計算資源消耗過于巨大。因此，本文進行了實驗室縮尺的下?lián)舯┝黠L(fēng)場仿真模擬，以探明下?lián)舯┝黠L(fēng)場特性以及地表粗糙程度對風(fēng)場的影響。

1 下?lián)舯┝黠L(fēng)場模型

目前，下?lián)舯┝鞣抡婺M的風(fēng)場模型主要歸結(jié)為兩種，一種是環(huán)形渦流模型(圖1a)，另一種是沖擊射流模型(圖1b)。沖擊射流模型是將射流垂直噴射向地表，氣流下沉過程中同時也產(chǎn)生環(huán)形渦流，協(xié)同效果下模擬下?lián)舯┝鞯娘L(fēng)場特性效果良好[6,7]，彌補了環(huán)形渦流模型對風(fēng)場內(nèi)核以及外圍區(qū)域模擬效果不理想的缺點。因此，本文選取沖擊射流風(fēng)場模型進行研究。

圖1 下?lián)舯┝黠L(fēng)場模型

2 數(shù)值模型

2.1 控制方程

本文選用大渦模擬的方法進行仿真模擬。如上文所述，仿真模擬中直接模擬計算大尺度渦旋，對于低于網(wǎng)格解析度的部分小尺度渦旋采用亞格子模型考慮其對大尺度渦旋的影響，采用Boussinesq假說以及標(biāo)準(zhǔn)Smagorinsky-Lilly模型計算亞格子應(yīng)力(Sub-Grid Scale，SGS)。過濾亞格子渦旋后的時間依賴Navier-Stokes方程(笛卡爾坐標(biāo)系)如下所示：

(1)

(2)

(3)

(4)

(5)

(6)

為了與光滑地面情況的風(fēng)場作對比，本文選取了浸沒邊界法(Immersed Boundary Method，IBM)來實現(xiàn)對存在建筑物的地面區(qū)域的模擬。這種方法直接在Navier-Stokes方程中添加地面粗糙引起的阻力項，進而可以有效避免建立粗糙元的幾何外形的步驟，公式如下：

(7)

(8)

空間離散采用有限體積法，二階中心差分格式用于對流項與粘性項，二階隱式格式用于非穩(wěn)態(tài)項的時間推進，SIMPLE算法用于壓強速度解耦，見文獻[9]，求解器采用Fluent 6.3.26[10]。

2.2 幾何模型

選取一個完整的三維圓柱體作為計算域，如圖2所示。這使流場演變的過程不會受到限制，模擬結(jié)果更為精確。選取幾何縮尺為1∶1000，取沖擊射流的直徑Dj為1 m。所選用的圓柱體計算域的尺寸如圖2a所示，計算域直徑為8Dj，豎直方向高度為4Dj。其中沖擊射流入口位于豎直高度2Dj處。沖擊射流通道也為圓柱體，延伸至圓柱體計算域上部圓面。根據(jù)下?lián)舯┝鞯娘L(fēng)速特征以及選取的計算域尺寸，本文采用Mason等[11]推薦的速度縮尺1∶3，沖擊射流的速度取為20 m/s。對于粗糙地表條件，動量平衡方程中的阻力項將加入到距離地表10 cm的范圍內(nèi)。

圖2 大渦模擬計算域及其網(wǎng)格分布

2.3 網(wǎng)格分布與邊界條件

使用GAMBIT劃分三維結(jié)構(gòu)化網(wǎng)格，網(wǎng)格數(shù)量為594000。在網(wǎng)格劃分過程中，首先在水平截面中心按照正方形劃分網(wǎng)格，再在圓周方向上劃分出面網(wǎng)格，最后完成體網(wǎng)格的劃分，以確保網(wǎng)格結(jié)構(gòu)質(zhì)量良好。貼地網(wǎng)格高度為0.1 mm，整個計算域最大網(wǎng)格高度為10 cm，最小水平網(wǎng)格尺寸為1 mm，最大水平網(wǎng)格尺寸為10 cm，豎直與水平網(wǎng)格增長率均小于1.05。沖擊射流通道與計算域相交的圓柱面為無剪切應(yīng)力滑移墻面，排除沖擊射流在進入計算域形成下?lián)舯┝髑皩ζ涞挠绊憽A柱體計算域中，底面為壁面邊界條件，圓柱面及上表面為壓力出口。計算時間步長為0.0001 s，總計算時間為50 s。

3 光滑地表流場特性

3.1 數(shù)據(jù)對比

將本文數(shù)值模擬結(jié)果與以往實測風(fēng)場以及仿真模擬結(jié)果以最大風(fēng)速做歸一化處理，繪制風(fēng)速剖面并進行對比，見圖3(圖中，z為高度；ur為徑向速度；urpeak為最大徑向速度)?？梢郧宄乜闯霰緮?shù)值模擬與以往的研究結(jié)果[1,2,4,12,13]吻合良好，由此驗證了本文湍流模型、網(wǎng)格分布、邊界條件以及時間步長等設(shè)置的可靠性。

圖3 風(fēng)速對比

3.2 風(fēng)場隨時間的演變

由于風(fēng)場的對稱性，本文選取一個豎向截面進行研究。圖4(圖中Tn為歸一化時間)為所選取的特征時刻下?lián)舯┝黠L(fēng)場渦量云圖，以此可以清楚地展現(xiàn)出流場演變過程。圖4a中，沖擊射流產(chǎn)生的主環(huán)形渦流開始沖擊地面，地面邊界處已經(jīng)產(chǎn)生次生渦流；圖4b中環(huán)形渦流正沖擊地面，地面邊界處產(chǎn)生次生渦流；圖4c中，主環(huán)形渦流脫落，在徑向上產(chǎn)生位移；圖4d中，主環(huán)形渦流徑向移動；圖4e中，主環(huán)形渦流徑向移動較大，并豎直向上移動；圖4f為風(fēng)場末期，渦流逐漸消散。

圖4 風(fēng)場渦量云圖

3.3 風(fēng)場分析

為了降低坐標(biāo)系對數(shù)據(jù)的影響，將速度分解為軸向速度和徑向速度兩部分，并且僅著重分析變化情況復(fù)雜的徑向速度特征，著眼四個特征位置處的風(fēng)速剖面特征進行分析。在徑向風(fēng)速變化最為復(fù)雜的半徑1Dj～2.5Dj區(qū)域中選取了四個特征位置繪制風(fēng)速剖面圖，如圖5，半徑R分別為1Dj，1.5Dj，2Dj，2.5Dj，而繪制風(fēng)剖面圖所選取的時間點是四個位置處分別達到最大徑向速度的時間點，Tn分別為7.8，8.8，11，12.8對應(yīng)上文中相同時間點的渦量圖?？梢郧宄乜闯霎?dāng)R=1Dj達到最大徑向風(fēng)速時，主環(huán)形渦流中心已經(jīng)過該位置，但渦流并未完全離開該位置。此時剖面處上部接近渦流中心，因此風(fēng)速較低；剖面下部受主環(huán)形渦流沖擊地面的后續(xù)渦流作用，速度較大。在R=1.5Dj處達到最大徑向速度時，剖面正好經(jīng)過渦流中心，這也正是剖面上端風(fēng)速出現(xiàn)負值的原因。剖面下端存在次生渦流，其附近風(fēng)速增強到1.6Vj(Vj為Dj處風(fēng)速)。下?lián)舯┝鳑_擊地面后，會產(chǎn)生沿地面邊界流動的氣流。該氣流受主環(huán)形渦流的負壓力影響，劇烈變動后，形成了與主渦流方向相反的次生渦流。主環(huán)形渦流與次生渦流協(xié)同作用下，近地面風(fēng)速得到較大增強。而R=2Dj處達到最大徑向速度時，主環(huán)形渦流發(fā)生了較明顯的上移，同時次生渦流強度逐漸減弱，因此該處最大風(fēng)速低于R=1.5Dj處。同時，主環(huán)形渦流的上移也導(dǎo)致未被次生渦流增強的位置處風(fēng)速略小，也就表現(xiàn)為圖中風(fēng)速曲線凹陷部分。而R=1Dj位置處，風(fēng)速又得到增強是副環(huán)形渦流沖擊地面的結(jié)果。由于R=2Dj處距風(fēng)場中心較遠，能量耗散，最大風(fēng)速也較低，次生渦流明顯消散，因此近地面風(fēng)速沒被明顯增強。

圖5 最大徑向風(fēng)速剖面

4 地表粗糙影響

圖6為統(tǒng)計了模擬的全部時刻0～0.3Dj高度內(nèi)的風(fēng)速后繪制的最大風(fēng)速包絡(luò)圖。可以清楚地看出地面粗糙對最大徑向風(fēng)速的影響。不加設(shè)地面粗糙時，最大徑向風(fēng)速達到1.8Vj并且發(fā)生在0.25Dj高度處；加設(shè)地面粗糙后，近地面風(fēng)速減小至0.2Vj，并在近地面范圍內(nèi)出現(xiàn)陡坡，然后風(fēng)速隨高度增加迅速增大至1.4Vj，在0.1Dj～0.18Dj之間分布均勻。結(jié)合模型建立的部分，可以分析出速度陡坡是由于地表加設(shè)的粗糙層引起的，其高度與粗糙層高度一致。徑向風(fēng)速達到最大值后，無論地表粗糙與否，徑向風(fēng)速沿高度的衰減率基本一致。綜上，在加大地表粗糙后，最大徑向風(fēng)速減小，并且出現(xiàn)的位置上移。這與Mason等[1]采用方法物理模擬的實驗結(jié)果一致。

圖6 最大徑向速度包絡(luò)圖

圖7是統(tǒng)計豎向0≤Z≤0.2Dj、徑向0.5Dj≤R≤3Dj范圍內(nèi)各點最大徑向速度后繪制的等高線圖。圖7a中出現(xiàn)了兩個峰值，分別位于徑向0.9Dj～1.0Dj和1.2Dj～1.6Dj之間，并且前者峰值分布范圍小于后者。結(jié)合前文渦量圖可以對此做出解釋，兩峰之間即為主環(huán)形渦流沖擊地面的中心，因此較大的徑向風(fēng)速出現(xiàn)在該位置兩側(cè)；后者分布范圍更大的原因正是次生渦流的增強作用以及渦流的豎向移動。加設(shè)粗糙情況下最大徑向風(fēng)速整體上低于光滑情況，這與前文結(jié)論一致，峰值分布在1.1Dj～1.4Dj和0.07Dj～0.18Dj之間，呈狹長形狀。如此不同的分布結(jié)果是兩種地面情況中主環(huán)形渦流的移動路徑不同引起的。

圖7 最大徑向速度等高線

為研究湍流部分的風(fēng)場，選取了湍流強度Iur作為依據(jù)，湍流強度計算方法與Holmes等[14]的方法一致，其計算公式如式(9)。這種方法相當(dāng)于把下?lián)舯┝鞯耐牧鞑糠之?dāng)作在達到最大速度時刻前后的一小段時間內(nèi)分段靜止的過程。

(9)

圖8 湍流強度等高線

5 結(jié) 論

本文利用大渦模擬方法完成了1∶1000縮尺的下?lián)舯┝黠L(fēng)場的CFD(Computational Fluid Dynamics)數(shù)值仿真模擬，研究了下?lián)舯┝鞯难葑冞^程和風(fēng)場特性，得到以下結(jié)論：

(1)下?lián)舯┝鹘孛骘L(fēng)場復(fù)雜多變的主要原因在于環(huán)形渦流與地面沖擊后流場劇烈的變化，地面粗糙會加劇流場變化的劇烈程度以及發(fā)展速度；

(2)下?lián)舯┝鬏S向速度特征較為簡單，主要受沖擊射流的流動路徑影響，在沖擊射流出流區(qū)域至地面范圍內(nèi)軸向速度較大，與沖擊射流的速度方向一致為負向，風(fēng)速大小也相近。正軸向速度主要出現(xiàn)在渦流環(huán)形軸外側(cè)，相對較??；

(3)下?lián)舯┝鲗恿鞑糠謴较蛩俣茸兓謩×遥墉h(huán)形渦流和近地面次生渦流影響較大，主要表現(xiàn)為在近地面迅速增大然后迅速減小，最大徑向風(fēng)速可達1.6Vj左右，最大風(fēng)速的出現(xiàn)位置隨著距離中心點的距離增加有略微上移，地面粗糙程度增大也會導(dǎo)致最大風(fēng)速出現(xiàn)的位置出現(xiàn)較大上移。