楊素珍

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

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

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

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

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

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

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

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

圖2 實測風(fēng)頻域功率譜

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

2 實測風(fēng)修正Davenport模型建模

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

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

(1)

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

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

(2)

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

(3)

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

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

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

(1)對函數(shù)(2)在P0點進行泰勒級數(shù)展開，并忽略高階項目，得到近似一階表達式

(4)

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

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

(5)

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

(6)

為最小值。

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

(7)

再將n(n>4)個數(shù)據(jù)點代入方程(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)

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

3 改進型多點脈動風(fēng)諧波合成法

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

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

(10)

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

(11)

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

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

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

(12)

式中

(13)

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

θjm(ωml)+φml]

(14)

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

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

Δt≤2π/(2ωu)

(15)

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

(16)

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

(17)

式中

(18)

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

4 實驗研究

對于修正Davenport模型建模，本文以某時段臺風(fēng)“杜鵑”的風(fēng)速為樣本，計算得到相應(yīng)的頻率值和功率值作為最小二乘法求解的數(shù)據(jù)點，選取待定系數(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模型與實測風(fēng)的功率譜進行對比，結(jié)果如圖3所示。

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

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

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

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

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

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

(19)

取風(fēng)速時程點數(shù)N=2048，截斷頻率ωu=10π，利用式(13)和式(16)，風(fēng)模擬結(jié)果如圖5和圖6所示。

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

常規(guī)與改進型多點脈動風(fēng)諧波合成法的上述臺風(fēng)時間序列模擬過程耗時對比如表1所示。

圖5 模擬臺風(fēng)風(fēng)速時程譜

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

表1 常規(guī)與改進型多點脈動風(fēng)諧波合成法耗時對比s

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

5 結(jié)束語

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