下擊暴流非平穩脈動風速數值模擬

李錦華，吳春鵬，陳水生

(華東交通大學 土木建筑學院，南昌 330013)

雷暴天氣會形成局部強下沉氣流猛烈沖擊地面，產生突發性、破壞性強風即下擊暴流。所謂下擊暴流，即能在地面產生17.9 m/s以上輻散風的強烈下沉氣流[1]。Chay等[2]從本質上詳細闡述過下擊暴流風速的形成方式。在雷暴環境中尺度較小的微下擊暴流發生頻率較高，在近地面產生的最大風速達75 m/s[3]。破壞性極大，屬近地面災害性強風。建筑物主體及圍護結構均會產生嚴重損傷、破壞，甚至垮塌[4-5]。

通過現場觀測、實驗模擬、理論分析及數值仿真等方法對下擊暴流展開廣泛研究。對下擊暴流數次實測發現，下擊暴流風速場中最大風速出現在近地面，而距地面較高處風速較小，與傳統大氣邊界層風速場分布明顯不同[3]。Holmes等[6]通過碰撞射流理論研究下擊暴流行進過程中某一固定位置產生的平均風速提出隨時間變化的平均風速模型。Wood等[7]通過下沉氣流的風洞實驗數據提出與計算流體動力學分析模型吻合的風速場平均風速模型。Savory等[8]用Holmes模型研究下擊暴流作用的輸電塔破壞。由于該模型未考慮風速場的隨機波動成分，會低估下擊暴流風荷載作用的結構動力響應。下擊暴流風速具有較強非平穩性。瞿偉廉等[9]將下擊暴流分解為確定的時變平均風速與調制的非平穩脈動風速進行數值模擬研究。李春祥等[10]采用Deodatis的均勻調制非平穩隨機場模擬方法模擬下擊暴流非平穩脈動風速。張文福等[11]通過時變函數均勻調制基于AR模型的平穩脈動風速生成具有空間相關性的下擊暴流非平穩脈動風速。據進化譜理論[12]，非平穩脈動風速模擬應通過與時間頻率有關的非均勻調制函數對功率譜進行調制獲得非平穩風速的時變功率譜即進化譜，再通過非平穩隨機過程模擬方法進行模擬。為簡化，通常將非均勻調制函數設為僅與時間有關的均勻調制函數，將非平穩隨機過程模擬轉化為均勻調制函數對平穩隨機過程的調制，可避免據時變功率譜模擬非平穩隨機過程難度。

目前，對非平穩脈動風速模擬通常將時變平均風速函數設為平穩脈動風速的均勻調制函數無相關理論依據。為此，李錦華等[13]將非平穩脈動風速離散成若干段在足夠短時間Δt內近似為平穩脈動風速的短時時間序列，嚴格推導出與時間、頻率有關的非均勻調制函數及相應非平穩脈動風速進化譜。

本文據進化譜理論，通過下擊暴流時變平均風速模型及與時間頻率有關的非均勻調制函數對功率譜進行調制獲得下擊暴流非平穩脈動風速時變功率譜；將平穩隨機過程AR模型考慮時變特征建立非平穩隨機過程的TAR時變模型；通過時變功率譜、TAR時變模型實現下擊暴流非平穩脈動風速的有效模擬。

1 下擊暴流時變平均風速

(1)

式中：f(t)為時間函數，最大值為1；V(z)為最大平均風速的豎向分布函數。

三種最大平均風速豎向風剖面[14]為

Vicroy豎向分布模型

V(z)=1.22[e-0.15z/zmax-e-3.2175z/zmax]Vmax

(2)

Wood豎向分布模型

(3)

Oseguera 、Bowles豎向分布模型

(4)

式中：Vmax為豎向分布風速中最大風速；zmax為最大風速所處高度位置；erf為誤差函數；δ為最大風速一半所處高度位置；r為距離下擊暴流風場中心徑向位置；R為下擊暴流風場輻射半徑；z*為大氣邊界層外特征高度；ε為大氣邊界層內特征高度；λ為尺度因子。

在碰撞射流理論中徑向風速為軸對稱，在某高度處徑向風速[6]可表示為

(5)

式中：Vr,max為風速場中某高度處最大風速；rmax為最大風速點與下擊暴流中心水平距離；r為觀測點與下擊暴流中心距離；Rr為徑向長度比例系數；Π為強度系數，表明雷暴強度隨時間變化。

考慮下擊暴流自身在移動，故觀測點與雷暴中心距離rt為隨時間變化，觀測點平均風速可表示為徑向風速Vr與下擊暴流平移速度Vo的矢量合成，即

Vp(t)=Vr(rt)+Vο

(6)

時間函數f(t)可定義為空間點任意時刻平均風速與最大平均風速比值，表示為豎向風剖面隨時間變化的函數，即

f(t)=Vp(t)/max|Vp(t)|

(7)

2 下擊暴流非平穩脈動風速

與大氣邊界層風速不同，下擊暴流的平均風速具有明顯的時變特征，其脈動風速具有較強的非平穩特征。基于非平穩隨機過程的進化譜理論，下擊暴流的非平穩脈動風速的模擬，可首先通過一調制函數對某一功率譜進行調制獲得時變功率譜，然后通過非平穩隨機過程模擬算法進行模擬。該模擬過程不僅需要獲得有效的調制函數，而且還需建立非平穩隨機過程的模擬算法。

2.1 非均勻調制函數

據非平穩隨機過程數值模擬理論，非平穩脈動風速的數值模擬關鍵是獲得時變功率譜。時變功率譜即進化譜的獲取可通過對現場實測非平穩風速進行時變譜估計或據進化譜理論通過調制函數對功率譜進行調制獲得。現場實測對大多數值模擬研究而言并不可行，而更方便可行的方法為用功率譜調制。為獲得有效調制函數，李錦華等[13]將非平穩脈動風速離散成若干段在足夠短時間Δt內可近似為平穩脈動風速的短時時間序列，基于Kaimal譜推導出與時間、頻率有關的非均勻調制函數，簡稱Kaimal調制函數。

(8)

Davenport調制函數

(9)

Harris調制函數

(10)

Simiu調制函數

(11)

2.2 進化譜

據進化譜理論[12]，零均值非平穩隨機過程f(t)可表示為

(12)

式中：A(ω,t)為非均勻調制函數；Z(ω)為正交過程，且滿足

E[dZ(ω)]=0

(13)

E[dZ*(ω)dZ(ω′)]=J(ω)δ(ω-ω′)dωdω′

(14)

式中：*表示共軛；δ為狄拉克函數。

非平穩隨機過程均值為

(15)

相關函數為

Rff(t,τ)=E[f*(t)f(t+τ)]=

τ)e-iωteiω′(t+τ)E[dZ*(ω)dZ(ω′)]=

(16)

τ=0時

(17)

因此，進化譜S(ω,t)可通過時頻函數A(ω,t)對功率譜非均勻調制表示[16]，即

S(ω,t)=|A(ω,t)|2J(ω)

(18)

據式(17)、(18)可建立非平穩隨機過程相關函數與時變功率譜間關系為

Rff(t,τ)=E[f*(t)f(t+τ)]=

(19)

2.3 非平穩脈動風速TAR時變模型建立

在平穩隨機過程模擬中，較成熟方法主要有譜表示及線性濾波。為進一步建立非平穩隨機過程模擬方法，Liang等[16]詳細推導考慮時變特征的非平穩隨機過程譜表示方法。在平穩隨機過程模擬中譜表示方法模擬精度較高，但模擬效率較低。對大型工程的風場模擬往往需在譜表示方法基礎上改進模擬效率。將譜表示成考慮時變特征模擬非平穩隨機過程，模擬效率更難滿足實際工程應用。線性濾波法因計算量小、速度快而廣泛用于隨機過程模擬，但其主要針對平穩高斯隨機過程模擬。李錦華等[17]考慮隨機過程的非高斯特性，實現基于線形濾波法中AR、ARMA模型非高斯隨機過程模擬。本節將繼續對線性濾波法中的AR模型考慮時變特征，建立非平穩隨機過程的TAR時變模型。

設TAR時變模型階數為p，非平穩隨機過程模擬樣本點數為N，采樣時間間隔為Δt，基于TAR(p) 時變模型的非平穩隨機過程模擬公式可表示為

(20)

在t=t′時刻，有

(21)

確定Ai(t′)可將式(21)兩邊同時右乘f(t′-jΔt)，并取數學期望，有

iΔt)f(t′-jΔt)]+E[L(t′)w(t′)f(t′-jΔt)]

(22)

L(t′)Rwf(t′,jΔt)，即

jΔt-iΔt)+L(t′)Rwf(t′,jΔt)

(23)

對右二項，Rwf(t′,jΔt)可理解為w(t′)與f(t′-jΔt)的互相關函數。w(t)為該系統輸入，而f(t)為系統輸出。當前輸出只依賴于當前與過去的輸入，與將來輸入無關，因此w(t′)與f(t′-jΔt)互相獨立，故

Rwf(t′,jΔt)=0

(24)

將式(24)代入式(23)，可確定Ai(t′)，即

(25)

展開式為

確定L(t′)可將式(21)兩邊同時右乘w(t′)，并取數學期望，有

(26)

當前輸入w(t′)與過去輸出f(t′-iΔt)互相獨立，則

Rfw(t′-iΔt,iΔt)=0

(27)

將式(27)代入式(26)，得

Rfw(t′,0)=L(t′)Rww(t′,0)

(28)

據白噪聲特性，有Rww(t′,0)=1，則

Rfw(t′,0)=L(t′)

(29)

將式(21)兩邊同時右乘f(t′)，并取數學期望，有

-iΔt)+L(t′)Rwf(t′,0)

(30)

因

Rwf(t′,0)=E[w(t′)f(t′)]=

E[f(t′)w(t′)]=Rfw(t′,0)

(31)

將式(29)、(31)代入式(30)，有

(32)

令A0(t′)=-1，則

(33)

因此，時變模型系數Ai(t′)、L(t′)(t′=1,2,3,…)，可分別據式(25)、(33)確定。

3 數值模擬

3.1 下擊暴流時變平均風速模擬

圖1 下擊暴流徑向風速與平移速度矢量合成

位于同高度的6個下擊暴流風速觀測點P1～P6見圖1。初始時刻，觀測點與下擊暴流中心距離均為r0=3 500 m，方位角為θ01=0°、θ02=15°、θ03=30°、θ04=45°、θ05=60°、θ06=90°。考慮下擊暴流自身的移動，則某高度觀測點與雷暴中心距離rt隨時間變化，觀測點平均風速可表示為徑向風速Vr與下擊暴流平移速度Vo的矢量合成:

Vp(t)=

(34)

下擊暴流平均風速豎向分布模型V(z)采用Vicroy模型，下擊暴流時變平均風速可表示為

(35)

考慮下擊暴流豎向分布風速中最大風速Vmax=80 m/s，所處高度位置zmax=67 m；風速場中某高度處徑向最大風速Vr,max取47 m/s，與下擊暴流中心水平距離rmax取1 000 m，徑向長度比例系數Rr取700 m，雷暴強度∏隨時間變化[18]可表示為

(36)

下擊暴流平移速度取Vo=8 m/s。

圖2 下擊暴流時變平均風速

各觀測點時變平均風速見圖2，對徑向距離相同各觀測點，時變平均風速波動隨方位角的增大而減弱。由式(8)～式(11)知，非平穩脈動風速功率譜非均勻調制函數與時變平均風速有關，因此時變平均風速波動性將影響非平穩脈動風速功率譜的時變性。考慮觀測點P1的時變平均風速具有較強波動性，本文基于建立的TAR時變模型進一步對該位置處豎向不同點非平穩脈動風速進行模擬。

3.2 基于TAR時變模型的下擊暴流非平穩脈動風速模擬

現采用建立的TAR時變模型來模擬下擊暴流行進路線上初始距離雷暴中心3 500 m處的豎向(Z1=15 m，Z2=25 m，Z3=45 m)三個不同點的非平穩脈動風速。模擬采用的上限截止圓頻率取為4 rad/s，模擬的時間長度為800 s，TAR時變模型階數P=15。下擊暴流非平穩脈動風速的時變功率譜采用Kaimal非均勻調制函數對Kaimal譜進行調制。三個不同點處的下擊暴流Kaimal非均勻調制函數見圖3。

圖3 下擊暴流Kaimal非均勻調制函數

圖4 下擊暴流風速時程與非平穩脈動成分

圖5 下擊暴流非平穩脈動風速時變功率譜

不同點間的相干函數采用Davenport相干函數模型，即

(37)

式中：zi，zj分別為第i，j點豎向空間坐標；Cz為空間豎向兩點間衰減系數；Ui，Uj分別為第i，j點平均風速，取時變平均風速算術平均值。

基于TAR模型模擬的下擊暴流風速時程及脈動風速成分見圖4。由圖4看出，脈動風速振幅與時變平均風速大小有關，即時變平均風速越大脈動風速振幅越大，與實際風場特性吻合，非均勻調制函數的有效性得以驗證。為說明模擬下擊暴流脈動風速具有目標非平穩特征，本文對模擬100組下擊暴流脈動風速樣本進行時變功率譜估計見圖5。由圖5看出，樣本估計功率譜時變特征明顯，此與目標時變譜的時變特征吻合。

圖6 單樣本下擊暴流非平穩脈動風速統計平均估計譜與目標譜對比

圖7 單樣本下擊暴流非平穩脈動風速統計平均自相關函數與目標自相關函數對比

圖8 單樣本下擊暴流非平穩脈動風速統計平均自相關函數

4 結 論

由于計算簡化的需要，模擬非平穩脈動風速通常將調制函數設為僅與時間有關的函數，使非平穩脈動風速模擬轉化為時間函數對平穩脈動風速調制，但無相關理論依據。本文通過建立有效功率譜非均勻調制函數及非平穩脈動風速模擬方法，對下擊暴流非平穩脈動風速進行數值模擬。并通過模擬的下擊暴流非平穩脈動風速相關性、功率譜分析驗證該模擬的有效性。

[1] 李春祥,李錦華,于志強.輸電塔線體系抗風設計理論與發展[J].振動與沖擊, 2009, 28(10):15-25.

LI Chun-xiang,LI Jin-hua,YU Zhi-qiang.A review of wind-resistant design theories of transmission tower-line systems[J].Journal of Vibration and Shock, 2009, 28(10):15-25.

[2] Chay M T, Letchford C W.Pressure distributions on a cube in a simulated thunderstorm downburst-part a: stationary downburst observations[J].Journal of Wind Engineering and Industrial Aerodynamics, 2002, 90(7):711-732.

[3] Letchford C W, Mans C, Chay M T.Thunderstorms-their importance in wind engineering (a case for the next generation wind tunnel)[J].Journal of Wind Engineering and Industrial Aerodynamics, 2002, 90(12/15):1415-1433.

[4] Holmes J D.Modeling of extreme thunderstorm winds for wind loading of structures and risk assessment[C].Proceedings of 10th International Conference on Wind Engineering, Copenhagen, Denmark, 1999: 1409-1415.

[5] Sengupta A, Haan F L, Sarkar P P, et al.Transient loads on buildings and tornado winds in microburst[J].Journal of Wind Engineering and Industrial Aerodynamics, 2008, 96(10/11): 2173-2187.

[6] Holmes J D, Oliver S E.An empirical model of a downburst[J].Engineering Structures, 2000, 22(9):1167-1172.

[7] Wood G S, Kwok K C S, Motteram N A, et al.Physical and numerical modelling of thunderstorm downdraft[J].Journal of Wind Engineering and Industrial Aerodynamics, 2001, 89(6):535-552.

[8] Savory E, Gerard A R P, Zeinoddini M,et al.Modelling of tornado and microburst-induced wind loading and failure of a lattice transmission tower[J].Engineering Structures, 2001, 23(4):365-75.

[9] 瞿偉廉,王錦文.下擊暴流風荷載的數值模擬[J].武漢理工大學學報, 2008, 30(2):70-74.

QU Wei-lian,WANG Jin-wen.Numerical simulation of downburst wind loads[J].Journai of Wuhan University of Technology, 2008, 30(2):70-74.

[10] 李春祥,劉晨哲,申建紅,等.土木工程下擊暴流風速數值模擬的研究[J].振動與沖擊, 2010, 29(12): 49-54.

LI Chun-xiang,LIU Chen-zhe,SHEN Jian-hong,et al.Numerical simulations of downburst wind speeds in civil engineering[J].Journal of Vibration and Shock, 2010, 29(12): 49-54.

[11] 張文福,謝丹,劉迎春,等.下擊暴流空間相關性風場模擬[J].振動與沖擊,2013,32(10):12-16.

ZHANG Wen-fu,XIE Dan,LIU Ying-chun,et al.Simulation of downburst wind fueld with spatial correlation[J].Journal of Vibration and Shock, 2013,32(10):12-16.

[12] Priestley M B.Power spectral analysis of non-stationary random processes[J].Journal of Sound and Vibration, 1967, 6(1):86-97.

[13] 李錦華,李春祥,申建紅.非平穩脈動風速的數值模擬[J].振動與沖擊, 2009, 28(1):18-23.

LI Jin-hua,LI Chun-xiang,SHEN Jian-hong.Numerical simulation of non-stationary fluctuating wind velocity[J].Journal of Vibration and Shock, 2009, 28(1):18-23.

[14] Chen L, Letchford C W.A deterministic-stochastic hybrid model of downbursts and its impact on a cantilever structure[J].Engineering Structures, 2004, 26(5):619-629.

[15] Vicroy D D.Assessment of microburst models for downdraft estimation[J].Journal of Aircraft, 1992, 29(6):1043-1048.

[16] Liang J, Chaudhuri S R, Shinozuka M.Simulation of non-stationary stochastic processes by spectral representation[J].Journal of Engineering Mechanics ASCE, 2007, 133 (6): 616-627.

[17] Li Jin-hua, Li Chun-xiang.Simulation of non-gaussian stochastic process with target power spectral density and lower-order moments[J].Journal of Engineering Mechanics ASCE, 2012, 138(5):391-404.

[18] Chay M T, Albermani F, Wilson R.Numerical and analytical simulation of downburst wind loads[J].Engineering Structures, 2006, 28 (2): 240-254.