大跨度橋梁脈動風(fēng)場模擬的插值算法

祝志文， 黃 炎

(1.湖南大學(xué) 土木工程學(xué)院，長沙 410082；2. 湖南大學(xué) 風(fēng)工程與橋梁工程湖南省重點(diǎn)實(shí)驗(yàn)室，長沙 410082)

大跨度橋梁脈動風(fēng)場模擬的插值算法

祝志文1，2， 黃 炎1

(1.湖南大學(xué) 土木工程學(xué)院，長沙 410082；2. 湖南大學(xué) 風(fēng)工程與橋梁工程湖南省重點(diǎn)實(shí)驗(yàn)室，長沙 410082)

為減少諧波合成法中功率譜矩陣分解的計算量，基于雙索引頻率和譜矩陣特點(diǎn)，提出了譜解矩陣雙軸插值算法和遞歸插值算法。兩種方法均使風(fēng)場模擬計算效率大幅度提高，且得到的模擬風(fēng)場均與傳統(tǒng)諧波合成法在統(tǒng)計意義上相符；研究了插值節(jié)點(diǎn)數(shù)量及分布形式對風(fēng)場模擬的影響，認(rèn)為在譜解矩陣的頻率軸向采取前密后疏分段插值，雙索引頻率軸向采用均勻分布插值節(jié)點(diǎn)的形式，更適合于風(fēng)場模擬。通過對一實(shí)際大跨度橋梁風(fēng)場的模擬，驗(yàn)證了方法的有效性。

風(fēng)場模擬；大跨度橋梁；諧波合成法；插值算法

在時域內(nèi)進(jìn)行大跨度橋梁抖振響應(yīng)分析，首先需要在脈動風(fēng)場離散空間點(diǎn)上，生成所需的具有特定要求的脈動風(fēng)速時程。目前，脈動風(fēng)場數(shù)值模擬方法主要有基于三角級數(shù)的諧波合成法和基于線性濾波技術(shù)的回歸方法等[1]。諧波合成法以功率譜作為權(quán)系數(shù)，與一系列帶隨機(jī)相位的三角級數(shù)的加權(quán)和來逐漸逼近隨機(jī)過程，適用于指定譜特征的平穩(wěn)高斯隨機(jī)過程。該方法有恒幅諧波疊加法和加權(quán)振幅諧波疊加法[2]。加權(quán)振幅諧波疊加利用FFT(Fast Fourier Transform)算法可以明顯減小機(jī)時，但由于頻率均勻分布導(dǎo)致模擬曲線出現(xiàn)周期性。為此，SHINOZUKA等[3]引入雙索引頻率，將頻率微擾量均勻分布在頻率增量內(nèi)，這樣既可采用FFT算法，又可保證模擬曲線的各態(tài)歷經(jīng)性。CAO等針對大跨度橋梁風(fēng)場離散點(diǎn)沿橋向可均勻分布的特點(diǎn)，導(dǎo)出了譜分解的顯示表達(dá)式，減少了部分余弦項(xiàng)，提高了計算效率，但該方法不能考慮三維空間風(fēng)場三個方向之間的相關(guān)性。DING等提出譜分解的三次拉格朗日多項(xiàng)式插值近似建議，但未給出插值節(jié)點(diǎn)的具體布置形式。BAO等[4]引入BP神經(jīng)網(wǎng)絡(luò)方法擬合譜分解函數(shù)曲線，減少了譜分解的計算量。XU等[5]根據(jù)實(shí)際大跨度橋梁風(fēng)速監(jiān)測點(diǎn)少的特點(diǎn)，采用克里金插值，模擬了橋梁抖振分析所需其他未測風(fēng)場離散點(diǎn)的非平穩(wěn)脈動風(fēng)。

本文基于雙索引頻率和譜矩陣特點(diǎn)，提出了譜解矩陣雙軸插值算法和遞歸插值算法。通過對一實(shí)際大跨度橋梁風(fēng)場的模擬，驗(yàn)證了方法的合理性。

1 脈動風(fēng)場諧波合成理論

1.1 傳統(tǒng)諧波合成方法

根據(jù)隨機(jī)過程理論[6-9]，模擬所生成的樣本必須能準(zhǔn)確描述隨機(jī)過程的概率特性。對于平穩(wěn)隨機(jī)過程，脈動風(fēng)速的特性可用功率譜和相關(guān)函數(shù)予以描述，功率譜密度函數(shù)揭示了隨機(jī)振動過程的能量按頻率分布的規(guī)律，相關(guān)函數(shù)則反映各點(diǎn)脈動風(fēng)速之間在時間或空間的相互關(guān)系。三維空間n個風(fēng)場離散點(diǎn)的功率譜矩陣見式(1)。

(1)

在工程實(shí)際中，三維空間點(diǎn)之間的互譜密度通常用實(shí)函數(shù)表達(dá)，因而互譜矩陣S(ω)為實(shí)對稱矩陣。Sjk(ω)是互譜矩陣S(ω)的元素，其中ω為風(fēng)的脈動圓頻率，j,k=1,2,…,n分別表示所模擬的第j個和第k個風(fēng)場離散空間點(diǎn)。當(dāng)j=k時為脈動風(fēng)速的自譜密度，當(dāng)j≠k時為脈動風(fēng)速的互譜密度。互譜密度可由自譜密度和三維空間相干函數(shù)根據(jù)式(2)確定

(2)

其中，脈動風(fēng)速的三維空間相干函數(shù)可表示為

Cohjk(ω)=

(3)

對互譜矩陣S(ω)進(jìn)行Cholesky分解，即S(ω)=H(ω)H(ω)T，H(ω)為下三角矩陣，也是實(shí)數(shù)矩陣，H(ω)T是其轉(zhuǎn)置矩陣。

(4)

(5)

(6)

互譜矩陣S(ω)和譜分解矩陣H(ω)都是雙索引頻率的矩陣，若用四維矩陣表示則為一個n×n×N×n的四維矩陣，其中，第三維對應(yīng)頻率軸，第四維對應(yīng)雙索引頻率軸。顯然，隨著模擬點(diǎn)n的增多，這個四維矩陣呈冪級數(shù)增大，不僅存儲量會非常巨大，而且要進(jìn)行n×N次n×n階矩陣的Cholesky分解，計算量相當(dāng)大，甚至很難實(shí)現(xiàn)。

1.2 譜解法的插值方法

1.3 譜解法的遞歸算法

為了提高譜分解效率，將風(fēng)場模擬式(5)中的|Hjm(ωml)|項(xiàng)寫為矩陣形式

(7)

式(7)的第j行元素就是第j個脈動風(fēng)場離散點(diǎn)模擬所需要的j個譜矩陣分解元素。該矩陣第m列元素僅取式(4)第m個雙索引頻率對應(yīng)譜矩陣分解后的第m列。也就是說，對于每次功率譜矩陣分解，僅有一列為風(fēng)場模擬所需，其余n-1列元素為多余元素。為了減少不必要的譜分解，根據(jù)遞歸算法，對第m個雙索引頻率對應(yīng)的互譜矩陣[S(ωml)]進(jìn)行式(8)形式分塊，使其分塊后的子矩陣[S11]為m×m階方陣。

(8)

2 譜解法的應(yīng)用

在風(fēng)工程領(lǐng)域中，許多學(xué)者通過對自然風(fēng)進(jìn)行長期的觀測，基于大量的實(shí)測數(shù)據(jù)，建立了不同的風(fēng)譜模型，以反映自然風(fēng)脈動的頻率結(jié)構(gòu)。本文采用我國規(guī)范《公路橋梁抗風(fēng)設(shè)計指南》推薦的Simiu譜模擬順風(fēng)向脈動風(fēng)。

(9)

本文針對一座實(shí)際大跨度橋梁開展了風(fēng)場模擬。該橋中跨主纜跨度1 176 m，主纜中心距27 m；中跨加勁梁跨度1 005 m，加勁梁梁寬27.6 m、梁高8.65 m。該橋采用鋼桁加勁梁，吊桿順橋向中心距為14.5 m，相鄰吊桿間設(shè)有兩個桁架節(jié)段。兩側(cè)主塔分別高134 m、67 m。

傳統(tǒng)的大跨度橋梁抖振響應(yīng)分析，脈動風(fēng)荷載作用往往僅考慮主要的風(fēng)荷載作用構(gòu)件加勁梁，忽略了主纜和橋塔上作用的風(fēng)荷載，這樣可能會影響大跨度橋梁抖振響應(yīng)分析的精度，但由于僅在加勁梁上布設(shè)脈動風(fēng)場模擬的空間離散點(diǎn)，因而脈動風(fēng)場模擬的計算量明顯減少，從而對計算資源的要求降低。本文從考慮大跨度橋梁全部主要風(fēng)荷載作用構(gòu)件的角度出發(fā)，風(fēng)場離散空間點(diǎn)按照如下方式布置：沿順橋向中心線布置在主桁上弦節(jié)點(diǎn)處，其編號為1～137；兩側(cè)主塔分別沿高度均勻布置15個離散點(diǎn)，其編號為138～152；上下游主纜風(fēng)場離散點(diǎn)沿主纜布置，其編號為153～374。全橋共374個風(fēng)場離散點(diǎn)，其空間分布見圖1，其中，括號標(biāo)注為橫橋向?qū)ΨQ分布于另一側(cè)主纜上模擬點(diǎn)編號，這樣便于獲得風(fēng)荷載作用的橋梁全部主要構(gòu)件不同空間點(diǎn)的模擬風(fēng)場，目的是為開展全橋全部主要構(gòu)件在風(fēng)荷載作用下響應(yīng)分析。風(fēng)場參數(shù)見表1。

表1 脈動風(fēng)場模擬參數(shù)

圖1 某懸索橋三維風(fēng)場模擬點(diǎn)布置(單位：m)Fig.1 Spacial arrangement of wind field points on a suspension bridge (Unit: m)

為探討功率譜矩陣在頻率和雙索引頻率軸向插值節(jié)點(diǎn)的選取形式對譜分解矩陣插值后的影響，進(jìn)而影響到風(fēng)速樣本模擬。本文在譜分解矩陣的頻率和雙索引頻率軸向分別單獨(dú)進(jìn)行插值，各取5種工況(見表2)分別計算H(ω)與傳統(tǒng)方法譜解矩陣對比，其中，在頻率軸向?yàn)? 024個頻率點(diǎn)，在雙索引頻率軸向?yàn)?74個頻率點(diǎn)，插值步距定義為插值節(jié)點(diǎn)間的頻率點(diǎn)數(shù)量。由此為合理選擇在頻率軸向以及雙索引頻率軸上的插值節(jié)點(diǎn)布置形式提供依據(jù)。

表2 雙軸插值10種工況

2.1 插值節(jié)點(diǎn)分布形式的討論

2.1.1 頻率軸向插值

2.1.2 雙索引頻率軸向插值

由圖3可知，雙索引頻率軸的H(ω)值基本呈線性變化，這是由于雙索引頻率的本質(zhì)是為了延長模擬樣本周期性和保證各態(tài)歷經(jīng)性，而在頻率增量內(nèi)為均勻分布的微小擾動。5種插值工況與傳統(tǒng)方法的譜分解矩陣均保持一致，所以在雙索引頻率軸，插值節(jié)點(diǎn)可以選用均勻分布形式，插值步距可以選用較大間隔。

2.2 合理選點(diǎn)的雙軸插值法

圖2 頻率軸向插值結(jié)果Fig.2 Frequency axis interpolation

圖3 雙索引頻率軸向插值結(jié)果Fig.3 Double-indexing frequency axis interpolation

圖4 雙軸插值法頻率軸對比Fig.4 Frequency direction results

圖5 雙軸插值法雙索引頻率軸對比Fig.5 Double-indexing frequency direction results

如圖6所示，雙軸插值法與傳統(tǒng)方法模擬得到的風(fēng)速樣本整體統(tǒng)計特性較為吻合，所存在的差異主要是由于譜分解矩陣的高階項(xiàng)在高頻頻帶部分趨近于0，從空間意義上說，風(fēng)場空間相關(guān)性的高頻頻帶部分較低頻頻帶部分要弱，相距越遠(yuǎn)的離散點(diǎn)脈動風(fēng)速相關(guān)性越小，而插值造成了在高頻部分的震蕩偏差。但通過后面的功率譜、相關(guān)函數(shù)對比，知其仍可以接受，也可進(jìn)一步加密插值節(jié)點(diǎn)更趨近于傳統(tǒng)方法。

圖6 風(fēng)場統(tǒng)計值Fig.6 Statistic of wind field

進(jìn)一步從整體統(tǒng)計特性上，對比模擬風(fēng)場離散點(diǎn)的相關(guān)系數(shù)矩陣，其為對稱矩陣，主對角線表示各點(diǎn)自相關(guān)系數(shù)為1，為方便對比，將上三角取為傳統(tǒng)方法，下三角取為雙軸插值法(見圖7)。由圖7(a)整體對比可見兩者對稱性很好，且與式(3)表示的任意頻率點(diǎn)上各點(diǎn)相干函數(shù)矩陣波峰位置相契合。進(jìn)一步取主梁(見圖7(b))、主纜(見圖7(c))以及一側(cè)主纜與主梁(見圖7(d))進(jìn)行詳細(xì)對比，在工程中一般認(rèn)為相關(guān)系數(shù)在0.1～0.2以下則兩者基本為不相關(guān)。從圖7(b)可知主梁某離散點(diǎn)相鄰約20點(diǎn)外的相關(guān)系數(shù)下降到了0.1以下，即可認(rèn)為主梁相距超過約145 m則兩點(diǎn)相關(guān)性很小；從圖7(c)可知，主纜存在兩個明顯的波峰，中間波峰為自相關(guān)系數(shù)，另一個波峰為兩側(cè)主纜相鄰點(diǎn)所對應(yīng)的互相關(guān)系數(shù)；由圖7(d)可知除開自相關(guān)系數(shù)的波峰，另一個小的波峰中間高，兩側(cè)低，其正好對應(yīng)主梁與邊跨、中跨主纜的互相關(guān)系數(shù)。綜上可見，雙軸插值法的譜分解矩陣及模擬風(fēng)場樣本在整體統(tǒng)計特性上均與傳統(tǒng)方法基本保持一致。

圖7 風(fēng)場相關(guān)系數(shù)對比Fig.7 Correlation coefficient of wind field

為說明雙軸插值諧波合成法風(fēng)場模擬效果，以加勁梁所在平面為參考基準(zhǔn)面，分別給出了加勁梁跨中第65點(diǎn)、第66點(diǎn)，右側(cè)索塔103 m高處第150點(diǎn)，以及主纜跨中12 m高處第211點(diǎn)(布置見圖1)的風(fēng)速時程，雖然在時域內(nèi)風(fēng)速樣本的統(tǒng)計特性無法看出，但顯然風(fēng)速時程在相鄰的主梁第65點(diǎn)(見圖8(a))、第66點(diǎn)(見圖8(b))及主纜離散點(diǎn)第211點(diǎn)(見圖8(d))的變化趨勢較主塔第150點(diǎn)(見圖8(c))更為相似。進(jìn)一步在頻域內(nèi)對比模擬功率譜函數(shù)與目標(biāo)功率譜函數(shù)，為包絡(luò)住大跨度橋梁風(fēng)場卓越頻率，給出0.01～2 Hz頻帶的功率譜函數(shù)，并采用雙對數(shù)坐標(biāo)表達(dá)形式，其中模擬樣本互功率譜為復(fù)數(shù)形式，圖中幅值為其實(shí)部與虛部的模，由圖9可知，模擬樣本功率譜與目標(biāo)功率譜在整個頻帶均較為一致。由于兩點(diǎn)的相關(guān)性隨距離增大而逐漸衰減，相關(guān)函數(shù)僅給出臨近的互相關(guān)函數(shù)，同時，脈動風(fēng)速經(jīng)過一段時距后，也會與起始狀態(tài)的相關(guān)性逐漸衰減，所以僅給出時滯在200 s以內(nèi)的相關(guān)函數(shù)對比結(jié)果，由圖10可知，模擬相關(guān)函數(shù)與目標(biāo)函數(shù)變化趨勢基本一致，均隨著時間增加，相關(guān)性逐漸降低，約20 s后脈動風(fēng)速樣本與起始狀態(tài)基本沒有了相關(guān)性，即脈動風(fēng)速流經(jīng)約400 m后基本不相關(guān)，這與自然風(fēng)積分尺度在150～450 m相吻合。由此說明了雙軸插值方法模擬的風(fēng)速樣本的可靠性。

圖8 風(fēng)速時程曲線Fig.8 Wind speed samples

圖9 自/互功率譜密度函數(shù)Fig.9 Auto-/cross-spectral density

2.3 風(fēng)場模擬的遞歸插值法

所謂遞歸插值法，是根據(jù)“1.3”所述譜分解矩陣的特點(diǎn)，在雙索引頻率軸向采用遞歸算法。根據(jù)“2.2”所述，在頻率軸向采用分段插值方法。也就是說，在所選的頻率插值節(jié)點(diǎn)上，樣本合成所需的譜分解矩陣與傳統(tǒng)方法完全一致，而頻率插值節(jié)點(diǎn)布置形式根據(jù)“2.1.1”討論采用分段插值，插值后的譜分解矩陣與傳統(tǒng)方法近似相同。由圖6的均值和標(biāo)準(zhǔn)差對比可知，遞歸插值法模擬風(fēng)場效果與雙軸插值方法相同，這是因?yàn)樵陬l率軸向，兩種方法的插值節(jié)點(diǎn)布置相同，而在雙索引頻率軸向，譜分解矩陣值呈線性分布。也進(jìn)一步說明，在高階項(xiàng)高頻部分出現(xiàn)的振蕩偏差均由頻率軸向插值導(dǎo)致，可減小頻率軸向插值步距消除。由于遞歸插值法與雙軸插值法的模擬效果相同，僅給出遞歸插值法的模擬功率譜與目標(biāo)譜對比結(jié)果，由圖11可知，遞歸插值法模擬風(fēng)速樣本功率譜與目標(biāo)譜在整個頻帶上都吻合較好。

圖10 自/互相關(guān)函數(shù)Fig.10 Auto-/cross-correlation function

圖11 自/互功率譜密度函數(shù)Fig.11 Auto-/cross-spectral density

3 模擬時間對比

本文在CPU為3.40 GHz的PC機(jī)上，利用MATLAB編制了諧波合成傳統(tǒng)方法、雙軸插值法和遞歸插值法這三種方法，模擬了圖1所示大跨度懸索橋主梁、主塔、纜索共374個風(fēng)場離散點(diǎn)，其模擬時間對比見表3，雙軸插值法的相對較短，且從存儲量來看，譜矩陣在頻率軸上若采用雙精度浮點(diǎn)類型存儲，傳統(tǒng)方法為374×374×1 024的矩陣，需占用近800 MB內(nèi)存，而雙軸插值法為374×374×183的矩陣，只需150 MB內(nèi)存。同時，由第“2”節(jié)討論可知，雙軸插值法以及遞歸插值法在風(fēng)場模擬統(tǒng)計意義上，與傳統(tǒng)方法基本保持一致，綜上所述，本文建議采用雙軸插值法。

表3 三種方法的模擬時間對比

Tab.3 Simulation time of the three methods min

諧波合成法傳統(tǒng)方法雙軸插值遞歸插值模擬時間3391169

4 結(jié) 論

(1) 本文采用雙軸插值法以及遞歸插值法對三維空間大跨度橋梁脈動風(fēng)場進(jìn)行了模擬，減少了譜分解次數(shù)，提高了譜分解效率。對比傳統(tǒng)方法模擬的風(fēng)速樣本均值、標(biāo)準(zhǔn)差、相關(guān)系數(shù)矩陣、功率譜密度函數(shù)、相關(guān)函數(shù)，驗(yàn)證了這兩種方法的可靠性。

(2) 討論了插值節(jié)點(diǎn)分布形式對譜分解矩陣的影響，得出了在頻率軸采用前密后疏的分段插值，在雙參數(shù)頻率軸向采用較大步距的等間距插值，得到的譜分解矩陣與傳統(tǒng)方法較吻合。

(3) 譜分解矩陣值在雙索引頻率軸呈線性分布，若雙軸插值法和遞歸插值法在頻率軸向的插值節(jié)點(diǎn)布置相同，那么，兩者模擬效果也相似。在譜分解矩陣高階項(xiàng)高頻部分出現(xiàn)的振蕩誤差均來自于頻率軸向插值，可減小插值步距消除。

[ 1 ] CAO Y, XIANG H, ZHOU Y. Simulation of stochastic wind velocity field on long-span bridges[J]. Journal of Engineering Mechanics，2000,126(1):1-6.

[ 2 ] DING Q, ZHU L, XIANG H. An efficient ergodic simulation of multivariate stochastic processes with spectral representation[J]. Probabilistic Engineering Mechanics, 2011,26(2):350-356.

[ 3 ] LIANG J, CHAUDHURI S, SHINOZUKA M. Simulation of nonstationary stochastic processes by spectral representation[J]. Journal of Engineering Mechanics, 2007,133(6):616-627.

[ 4 ] 包龍生，劉克同，于玲,等. 大跨度橋梁空間脈動風(fēng)場的數(shù)值模擬[J]. 沈陽建筑大學(xué)(自然科學(xué)版)，2010，26(2)：238-243. BAO Longsheng, LIU Ketong, YU Ling, et al , Numerical simulation of spatial fluctuating wind field on long-span bridges[J]. Journal of Shenyang Jianzhu University (Natural Sciences), 2010,26(2):238-243.

[ 5 ] XU Y, HU L, KAREEM A. Conditional simulation of nonstationary fluctuating wind speeds for long-span bridges[J]. Journal of Engineering Mechanics, 2014,140(1):61-73.

[ 6 ] DEODATIS G. Simulation of ergodic multivariate stochastic processes[J]. Journal of Engineering Mechanics, 1996,122(8):778-787.

[ 7 ] 李永樂，周述華，強(qiáng)士中. 大跨度斜拉橋三維脈動風(fēng)場模擬[J]. 土木工程學(xué)報，2003，36(10):60-65. LI Yongle, ZHOU Shuhua, QIANG Shizhong. Simulation of three-dimensional fluctuating wind field for large span cable-stayed bridge [J].China Civil Engineering Journal,2003,36(10):60-65.

[ 8 ] 趙林.風(fēng)場模式數(shù)值模擬與大跨橋梁抖振概率評價[D]. 上海：同濟(jì)大學(xué)，2003.

[ 9 ] 丁泉順.大跨度橋梁耦合為顫抖振響應(yīng)的精細(xì)化分析[D]. 上海：同濟(jì)大學(xué),2001.

Interpolation algorithm for fluctuating wind field simulation of long-span bridges

ZHU Zhiwen1，2， HUANG Yan1

(1. College of Civil Engineering, Hunan University, Changsha 410082, China；2. Hunan Provincial Key Laboratory for Wind and Bridge Engineering, Hunan University, Changsha 410082, China)

In order to reduce computing time of power spectrum matrix decomposition in the harmonic superposition method, the biaxial-and recursive interpolation algorithms were respectively proposed based on the characteristics of double-indexing frequency and spectral matrix. The efficiency of the two methods was proved to be significantly improved, and from the statistical point of view, the simulated results were consistent with those using the traditional harmonic superposition method. The effects of the interpolation node number and arrangement form on the simulated results of wind field were also investigated. It was shown that a fine arrangement at front end and a coarse arrangement at rear end in frequency domain, and a uniform node distribution for double-indexing frequency are more suitable to wind field simulation. Through simulating a fluctuating wind field of a long-span suspension bridge, the effectiveness of the proposed methods was verified.

wind field simulation；long-span bridges；harmonic superposition method；interpolation algorithm

國家重點(diǎn)基礎(chǔ)研究發(fā)展計劃(2015CB057701；2015CB057702)；國家自然科學(xué)基金(51278191)；湖南省高校創(chuàng)新平臺開放基金(13K016)

2015-10-23 修改稿收到日期： 2016-03-02

祝志文 男，博士，教授，博士生導(dǎo)師，1968年生

黃炎 男，博士生，1987年生

U441

A

10.13465/j.cnki.jvs.2017.07.024