基于實(shí)測(cè)風(fēng)特性的臺(tái)風(fēng)模擬研究

楊素珍

(漳州職業(yè)技術(shù)學(xué)院 機(jī)械工程學(xué)院，福建 漳州363000)

臺(tái)風(fēng)往往給人類的生產(chǎn)和生活環(huán)境帶來(lái)巨大破壞[1]，造成大量的經(jīng)濟(jì)損失和人員傷亡，通過對(duì)人們活動(dòng)環(huán)境內(nèi)重要目標(biāo)，如建筑、橋梁、塔架和電線桿等，進(jìn)行臺(tái)風(fēng)風(fēng)載荷作用測(cè)試，研究物體風(fēng)致破壞的機(jī)理，從而找出減小臺(tái)風(fēng)破壞的措施。常規(guī)基于風(fēng)洞試驗(yàn)和現(xiàn)場(chǎng)實(shí)測(cè)的風(fēng)載荷作用檢驗(yàn)手段成本高、耗時(shí)長(zhǎng)；通過模擬臺(tái)風(fēng)的方法仿真實(shí)施目標(biāo)風(fēng)載荷作用[2-3]，可以較好的解決上述問題。因此研究臺(tái)風(fēng)模擬手段意義重大，引起國(guó)內(nèi)外學(xué)者的廣泛關(guān)注[4-9]。

風(fēng)的變化過程可以近似為一種平穩(wěn)隨機(jī)過程，一般可將其分解為穩(wěn)態(tài)平均風(fēng)和動(dòng)態(tài)脈動(dòng)風(fēng)。平均風(fēng)的風(fēng)速值和風(fēng)向不隨時(shí)間變化；脈動(dòng)風(fēng)則隨時(shí)間無(wú)規(guī)律地改變其風(fēng)速和風(fēng)向。脈動(dòng)風(fēng)動(dòng)力成分是強(qiáng)(臺(tái))風(fēng)的重要組成[5]，集中反映其隨機(jī)特性，該特性可以用紊流強(qiáng)度和功率譜密度函數(shù)等來(lái)描述。文獻(xiàn)[4]借助不同坡度三維山體模型風(fēng)洞試驗(yàn)，研究了山體各位置脈動(dòng)風(fēng)速譜的變化規(guī)律。文獻(xiàn)[5]利用臺(tái)風(fēng)“羅莎”實(shí)測(cè)數(shù)據(jù)，研究其高空風(fēng)參數(shù)特性，結(jié)合Kolmogrove理論，提出水平方向和垂直方向脈動(dòng)風(fēng)功率譜密度函數(shù)。上述基于實(shí)際風(fēng)速的方法成本高、耗時(shí)長(zhǎng)。文獻(xiàn)[6] 基于諧波合成法進(jìn)行海上風(fēng)機(jī)的脈動(dòng)風(fēng)速度時(shí)程模擬，考慮了過程順風(fēng)方向和垂直方向計(jì)算點(diǎn)的二維相干性。文獻(xiàn)[7]基于標(biāo)準(zhǔn)達(dá)文普(Davenport)模擬脈動(dòng)風(fēng)的風(fēng)速和風(fēng)壓時(shí)程曲線，并將模擬風(fēng)譜應(yīng)用于巨型框架的減振結(jié)構(gòu)設(shè)計(jì)。文獻(xiàn)[8]提出一種結(jié)合諧波合成法和本征正交分解的脈動(dòng)風(fēng)場(chǎng)數(shù)值模擬方法，來(lái)解決計(jì)算機(jī)內(nèi)存超限問題。文獻(xiàn)[9]采用線性濾波法模擬脈動(dòng)風(fēng)速，生成作用于上海中心大廈結(jié)構(gòu)上的脈動(dòng)風(fēng)載荷時(shí)程，為該大廈的風(fēng)結(jié)構(gòu)設(shè)計(jì)奠定基礎(chǔ)。上述基于脈動(dòng)風(fēng)風(fēng)速時(shí)程模擬方法所用的譜密度函數(shù)是通過大量實(shí)測(cè)數(shù)據(jù)的統(tǒng)計(jì)分析所獲得的經(jīng)驗(yàn)式，但在表示特定地區(qū)風(fēng)時(shí)，由于各地區(qū)風(fēng)移動(dòng)環(huán)境存在較大差異，經(jīng)驗(yàn)式不能完整反映該地區(qū)的實(shí)際風(fēng)特性。

為此，本文提出一種基于特定地區(qū)實(shí)際風(fēng)采集數(shù)據(jù)的風(fēng)模擬方法，該方法利用實(shí)測(cè)臺(tái)風(fēng)風(fēng)速數(shù)據(jù)和最小二乘法對(duì)標(biāo)準(zhǔn)Davenport功率譜密度函數(shù)進(jìn)行修正，再基于一類改進(jìn)型諧波合成法進(jìn)行風(fēng)模擬，在較完整反映特定地區(qū)風(fēng)特性的同時(shí)，保證模擬的快速性和精確性。本文首先分析福建某地區(qū)實(shí)測(cè)臺(tái)風(fēng)特性，然后詳細(xì)闡述修正Davenport模型數(shù)學(xué)建模和基于改進(jìn)型諧波合成法的多點(diǎn)脈動(dòng)風(fēng)模擬算法，最后給出實(shí)驗(yàn)情況。

1 特定地區(qū)實(shí)測(cè)風(fēng)特性分析

為了獲取特定地區(qū)實(shí)際臺(tái)風(fēng)的特性信息，通過專用風(fēng)采集裝置采集特定高度實(shí)際風(fēng)速樣本，以臺(tái)風(fēng)“杜鵑”為例，采集點(diǎn)為福建某城市，取得實(shí)際風(fēng)速樣本10000個(gè)，采樣頻率為0.1s，時(shí)間總長(zhǎng)度為1000s，其反映能量隨頻率分布的頻域峰值譜、功率譜分別如圖1和圖2所示。由圖1可以看出，臺(tái)風(fēng)“杜鵑”的風(fēng)能主要集中在0～3Hz的范圍內(nèi)，超過該范圍基本沒有風(fēng)能的分布，實(shí)際風(fēng)幅值的集中范圍情況與標(biāo)準(zhǔn)Davenport十分相似。

圖1 實(shí)測(cè)風(fēng)頻域幅值譜(截止頻率4Hz)

圖1對(duì)比結(jié)果表明，通過采集裝置采集到的風(fēng)速數(shù)據(jù)是真實(shí)可信的。

圖2 實(shí)測(cè)風(fēng)頻域功率譜

從圖2的功率譜上看，二者明顯偏離過大，吻合度不夠，可見，如果采用標(biāo)準(zhǔn)Davenport進(jìn)行風(fēng)模擬，無(wú)法完全體現(xiàn)采集所在地區(qū)風(fēng)的特點(diǎn)，需要結(jié)合實(shí)測(cè)風(fēng)對(duì)標(biāo)準(zhǔn)Davenport進(jìn)行修正，才能更好地反映該地區(qū)風(fēng)的真實(shí)特性。

2 實(shí)測(cè)風(fēng)修正Davenport模型建模

平穩(wěn)隨機(jī)過程可通過自功率譜密度函數(shù)來(lái)表征其能量在不同頻率上的分布，風(fēng)模擬的關(guān)鍵就是選取能夠反映風(fēng)時(shí)程數(shù)值特征的風(fēng)功率譜密度函數(shù)?；贙olmogrov瑞流理論，風(fēng)工程領(lǐng)域?qū)W者通過對(duì)自然風(fēng)的大量觀測(cè)和實(shí)測(cè)數(shù)據(jù)統(tǒng)計(jì)，已建立了許多滿足瑞流模型的風(fēng)功率譜密度函數(shù)，常用的有Kaimal譜、Simitu譜、ESDU譜和Davenport譜。由于觀測(cè)地區(qū)和數(shù)據(jù)統(tǒng)計(jì)方式的不一致性，不同譜密度函數(shù)所模擬的風(fēng)場(chǎng)頻率結(jié)構(gòu)各有差異，一般無(wú)法與特定地區(qū)的實(shí)際風(fēng)特性完全匹配。

瑞流風(fēng)速主要包括水平、縱向和豎向三個(gè)方向風(fēng)速，其中水平方向產(chǎn)生的風(fēng)壓對(duì)結(jié)構(gòu)的作用最直接，破壞性也最強(qiáng)，目前包括美國(guó)、德國(guó)、加拿大和中國(guó)在內(nèi)的大多數(shù)國(guó)家的風(fēng)載荷規(guī)范中，普遍采用水平脈動(dòng)風(fēng)譜Davenport譜[10]，該譜不隨高度變化，頻譜中高頻諧波分量明顯，參數(shù)配置簡(jiǎn)單，方便進(jìn)行工程應(yīng)用。Davenport譜的水平順向風(fēng)模型為

(1)

式中：x=1200f/U10；f為頻率(Hz)；U10為10m高度處的平均風(fēng)速；k為地面粗糙度系數(shù)(VonKaman常數(shù))。

為使Davenport模型能夠更真實(shí)反映特定地區(qū)的實(shí)際風(fēng)特性，本文提出一種修正Davenport模型，利用特定地區(qū)實(shí)測(cè)數(shù)據(jù)確定該表達(dá)式的待定系數(shù)，讓修正模型中包含特定地區(qū)的客觀變化和行為信息。定義修正Davenport模型的表達(dá)式為

(2)

式中P1、P2、P3、P4為待定系數(shù)。

(3)

式中：F(f)為頻域幅值譜，E(·)表示均值，T為風(fēng)速樣本持續(xù)的時(shí)間。

文中利用實(shí)測(cè)臺(tái)風(fēng)“杜鵑”風(fēng)速數(shù)據(jù)，求取風(fēng)速功率譜，得到求解修正Davenport模型系數(shù)所需的頻率值和功率值樣本，然后基于一類非線性最小二乘法擬合求解得到修正Davenport模型中的待定系數(shù)。

定義待定系數(shù)初值為P0=[P10，P20，P30，P40]，將該初值代入式(2)得到ψ0，則待定系統(tǒng)的求解步驟如下。

(1)對(duì)函數(shù)(2)在P0點(diǎn)進(jìn)行泰勒級(jí)數(shù)展開，并忽略高階項(xiàng)目，得到近似一階表達(dá)式

(4)

式中： ?i=Pi-Pi0，i=1,2,3,4，j=1,2,…,n表示用于最小二乘法擬合數(shù)據(jù)點(diǎn)個(gè)數(shù)序號(hào)；ψj0為序號(hào)j數(shù)據(jù)點(diǎn)初值函數(shù)值。

(2)建立帶求解系數(shù)非線性函數(shù)的近似線性關(guān)系表達(dá)式

(5)

(3)確定求解bi的最小二乘求解準(zhǔn)則為：使n個(gè)數(shù)據(jù)點(diǎn)(xj,yj)與非線性函數(shù)y距離ηj的平方和最小，即求取bi使得

(6)

為最小值。

(4)對(duì)方程(6)求偏導(dǎo)數(shù)，并令其等于零，可得

(7)

再將n(n>4)個(gè)數(shù)據(jù)點(diǎn)代入方程(7)，可建立超定方程組

Φ(X,Y,B)=0

(8)

式中：Y=(y1,y2,…,yn)；B=(b1,b2,b3,b4)；

X=(x11,…,xn1,x12,…,xn2,x13,…,xn3,x14,…,xn4)求解超定方程組中的B值即為最小二乘解。

(5)求得待定系數(shù)為

P=[P1,P2,P3,P4]=B+P0

(9)

對(duì)非線性函數(shù)(2)的近似線性形式(5)，可以采用多元線性擬合方法進(jìn)行求解，降低了計(jì)算復(fù)雜度。由于采用近似線性形式，初步計(jì)算的值存在一定誤差，為提高擬合精度，需將求解后的值重新作為初值，多次進(jìn)行上述求解步驟，直至滿足所需精度要求。

3 改進(jìn)型多點(diǎn)脈動(dòng)風(fēng)諧波合成法

離散傅里葉變換可以用來(lái)描述具有時(shí)間和空間特性的風(fēng)速時(shí)間序列，諧波合成法(譜表示法)的原理是采用一系列具有隨機(jī)頻率的余弦函數(shù)序列來(lái)模擬脈動(dòng)風(fēng)，并考慮結(jié)構(gòu)風(fēng)速的特點(diǎn)。由于同一時(shí)刻作用于結(jié)構(gòu)上的多個(gè)點(diǎn)的風(fēng)載荷不可能完全一樣，因此，在模擬風(fēng)速時(shí)程時(shí)，應(yīng)考慮其相關(guān)性，即采用多個(gè)互相關(guān)的隨機(jī)過程來(lái)描述。

風(fēng)的實(shí)測(cè)數(shù)據(jù)表明，風(fēng)速的脈動(dòng)部分可看作是具有零均值的高斯(Gauss)平穩(wěn)隨機(jī)過程，對(duì)具有不同坐標(biāo)的風(fēng)點(diǎn)(1，2，3，…，n)，其風(fēng)速uj(t)(j=1，2，3，…，n)的譜密度函數(shù)矩陣為

(10)

式中：S11(ω)，S22(ω)，…，Snn(ω)為不同坐標(biāo)點(diǎn)風(fēng)速功率譜；Sij(ω)(i≠j)為i點(diǎn)風(fēng)和j點(diǎn)風(fēng)的相關(guān)功率譜。風(fēng)速功率譜采用式(2)的修正Davenport，以讓模擬風(fēng)包含指定地區(qū)的風(fēng)特性，相關(guān)功率譜采用式(11)表達(dá)式。

(11)

式中：i≠j，Zi、Zj分別為i點(diǎn)風(fēng)和j點(diǎn)風(fēng)的高度值，Vi、Vj分別為i點(diǎn)風(fēng)和j點(diǎn)風(fēng)的平均風(fēng)值。

非負(fù)定矩陣[S(ω)]可以用Cholesky方法分解成如下形式

[S(ω)]=[H(ω)][H*(ω)]T

(12)

式中

(13)

此時(shí)，隨機(jī)過程uj(t)可以用式[14]進(jìn)行模擬(當(dāng)N趨于無(wú)窮大時(shí))

θjm(ωml)+φml]

(14)

式中：N為頻率域內(nèi)的數(shù)據(jù)采集數(shù)目；

φ1l、φ2l、…、φml為相互獨(dú)立的[0，2π]上均勻分布的隨機(jī)相位角序列；Δω=ωu/N，ωu為截?cái)囝l率，其大小通常由功率譜密度矩陣中各項(xiàng)與ω的函數(shù)關(guān)系而定。為避免頻率混淆，模擬的時(shí)間不長(zhǎng)，Δt須滿足式(15)。

Δt≤2π/(2ωu)

(15)

采用諧波合成法脈動(dòng)風(fēng)模擬盡管計(jì)算精度高，但需要在每個(gè)頻率上進(jìn)行大量運(yùn)算，極大地降低了算法的計(jì)算速度。本文通過引入快速傅里葉變換FFT[11]改進(jìn)常規(guī)諧波合成法來(lái)簡(jiǎn)化計(jì)算量，從而提高算法的計(jì)算效率。改進(jìn)后式(14)的表達(dá)式為

(16)

式中：η=1,2,3,…,2N×n-1；j=1,2,3,…,n；λ=1,2,3,…,n-1，為η/2N余數(shù)；Qjm(λΔt)的表達(dá)式為

(17)

式中

(18)

從式(18)可以知道Qjm(λΔt)是Djm(sΔω)的傅里葉變換，因此，可以用FFT進(jìn)行相關(guān)計(jì)算。

4 實(shí)驗(yàn)研究

對(duì)于修正Davenport模型建模，本文以某時(shí)段臺(tái)風(fēng)“杜鵑”的風(fēng)速為樣本，計(jì)算得到相應(yīng)的頻率值和功率值作為最小二乘法求解的數(shù)據(jù)點(diǎn)，選取待定系數(shù)初值P0=[1，1.5，1.5，2]，利用式(4)～式(9)，經(jīng)過多次循環(huán)求解，得到系數(shù)值為

P1=0.4696，P2=1.2721，P3=0.1023，P4=1.8459。

系數(shù)代入式(2)可得到修正Davenport模型。

將所建立的修正前后Davenport模型與實(shí)測(cè)風(fēng)的功率譜進(jìn)行對(duì)比，結(jié)果如圖3所示。

圖3 臺(tái)風(fēng)杜鵑與修正前后Davenport模型功率譜對(duì)比

由圖3可知，模型修正前，臺(tái)風(fēng)杜鵑的功率譜幅值顯然比標(biāo)準(zhǔn)Davenport大，二者的吻合度差，且頻率越大，標(biāo)準(zhǔn)Davenport偏離臺(tái)風(fēng)杜鵑功率譜的趨勢(shì)加大；而采用修正后的Davenport模型，二者重合度得到了較大的提高。

進(jìn)一步以某時(shí)段臺(tái)風(fēng)“蘇迪羅”的風(fēng)速為樣本，修正Davenport模型采用與臺(tái)風(fēng)“杜鵑”實(shí)驗(yàn)相同的系數(shù)值，實(shí)測(cè)風(fēng)與修正前后Davenport模型的功率譜，對(duì)比結(jié)果如圖4所示。

由圖4可以看出，修正后的Davenport模型與修正前相比，同樣具有更好的吻合度，因此，本文基于修正Davenport模型進(jìn)行風(fēng)速模擬，可以更好地體現(xiàn)實(shí)際風(fēng)速樣本采集地的風(fēng)特性。值得一提的是，通過增加風(fēng)速樣本數(shù)量，提高最小二乘法的擬合精度，可以進(jìn)一步提高吻合度。

圖4 臺(tái)風(fēng)蘇迪羅與修正前后Davenport模型功率譜對(duì)比

為了說明本文修正Davenport模型對(duì)臺(tái)風(fēng)的模擬效果，這里以12級(jí)臺(tái)風(fēng)為例，進(jìn)行臺(tái)風(fēng)時(shí)間序列模擬。地面粗糙度取B類，查風(fēng)速等級(jí)表得到Z1=10m高度的平均風(fēng)速為V1=30.9m/s，取風(fēng)速點(diǎn)數(shù)為n=2，第二點(diǎn)高度Z2=7m，根據(jù)指數(shù)風(fēng)剖面平均風(fēng)速沿高度變化的規(guī)律，可求取相關(guān)第二點(diǎn)的平均風(fēng)速值為

(19)

取風(fēng)速時(shí)程點(diǎn)數(shù)N=2048，截?cái)囝l率ωu=10π，利用式(13)和式(16)，風(fēng)模擬結(jié)果如圖5和圖6所示。

從圖5可以看出，臺(tái)風(fēng)隨時(shí)間變化過程具有明顯的脈動(dòng)風(fēng)特性，且風(fēng)速變化跨度大，符合實(shí)際測(cè)量強(qiáng)臺(tái)風(fēng)的變化特點(diǎn)。從圖6可以看出，模擬臺(tái)風(fēng)的功率譜與修正Davenport目標(biāo)譜吻合度高，說明風(fēng)速時(shí)間模擬序列較好地體現(xiàn)了指定地區(qū)實(shí)際臺(tái)風(fēng)的風(fēng)特性。

常規(guī)與改進(jìn)型多點(diǎn)脈動(dòng)風(fēng)諧波合成法的上述臺(tái)風(fēng)時(shí)間序列模擬過程耗時(shí)對(duì)比如表1所示。

圖5 模擬臺(tái)風(fēng)風(fēng)速時(shí)程譜

圖6 模擬臺(tái)風(fēng)功率譜與目標(biāo)譜比較

表1 常規(guī)與改進(jìn)型多點(diǎn)脈動(dòng)風(fēng)諧波合成法耗時(shí)對(duì)比s

由表1可見，基于快速傅里葉變換FFT的改進(jìn)型諧波合成法有效減小了算法的計(jì)算量，整個(gè)模擬過程耗時(shí)僅為常規(guī)諧波合成法的1/8，算法運(yùn)行實(shí)時(shí)性得到了較大的改善。

5 結(jié)束語(yǔ)

臺(tái)風(fēng)模擬是一種快速有效的分析臺(tái)風(fēng)環(huán)境內(nèi)物體風(fēng)致時(shí)程響應(yīng)的手段，通過實(shí)測(cè)風(fēng)數(shù)據(jù)分析表明，采用常規(guī)功率密度函數(shù)進(jìn)行脈動(dòng)風(fēng)模擬的方法，不能完全反映出特定地區(qū)的風(fēng)特性，本文通過利用福建某城市實(shí)測(cè)風(fēng)數(shù)據(jù)和最小二乘法對(duì)標(biāo)準(zhǔn)Davenport模型進(jìn)行再建模，獲得更能表征該地區(qū)風(fēng)特性的修正Davenport模型，在此基礎(chǔ)上，采用多點(diǎn)脈動(dòng)風(fēng)諧波合成法進(jìn)行臺(tái)風(fēng)模擬，并通過FFT方法改進(jìn)其計(jì)算效率，提高了與目標(biāo)譜的吻合度?；诒疚乃岱椒?，通過其他地區(qū)的實(shí)測(cè)風(fēng)數(shù)據(jù)，就可以實(shí)現(xiàn)能反映該地區(qū)風(fēng)特性的臺(tái)風(fēng)模擬，為準(zhǔn)確研究該地區(qū)風(fēng)致破壞奠定良好基礎(chǔ)。