大跨度橋梁脈動風場模擬的插值算法

祝志文， 黃 炎

(1.湖南大學 土木工程學院，長沙 410082；2. 湖南大學 風工程與橋梁工程湖南省重點實驗室，長沙 410082)

大跨度橋梁脈動風場模擬的插值算法

祝志文1，2， 黃 炎1

(1.湖南大學 土木工程學院，長沙 410082；2. 湖南大學 風工程與橋梁工程湖南省重點實驗室，長沙 410082)

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

風場模擬；大跨度橋梁；諧波合成法；插值算法

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

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

1 脈動風場諧波合成理論

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

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

(1)

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

(2)

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

Cohjk(ω)=

(3)

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

(4)

(5)

(6)

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

1.2 譜解法的插值方法

1.3 譜解法的遞歸算法

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

(7)

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

(8)

2 譜解法的應用

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

(9)

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

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

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

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

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

表2 雙軸插值10種工況

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

2.1.1 頻率軸向插值

2.1.2 雙索引頻率軸向插值

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

2.2 合理選點的雙軸插值法

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

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

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

圖7 風場相關系數(shù)對比Fig.7 Correlation coefficient of wind field

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

圖8 風速時程曲線Fig.8 Wind speed samples

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

2.3 風場模擬的遞歸插值法

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

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

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

3 模擬時間對比

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

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

Tab.3 Simulation time of the three methods min

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

4 結(jié) 論

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

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

(3) 譜分解矩陣值在雙索引頻率軸呈線性分布，若雙軸插值法和遞歸插值法在頻率軸向的插值節(jié)點布置相同，那么，兩者模擬效果也相似。在譜分解矩陣高階項高頻部分出現(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 ] 包龍生，劉克同，于玲,等. 大跨度橋梁空間脈動風場的數(shù)值模擬[J]. 沈陽建筑大學(自然科學版)，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 ] 李永樂，周述華，強士中. 大跨度斜拉橋三維脈動風場模擬[J]. 土木工程學報，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 ] 趙林.風場模式數(shù)值模擬與大跨橋梁抖振概率評價[D]. 上海：同濟大學，2003.

[ 9 ] 丁泉順.大跨度橋梁耦合為顫抖振響應的精細化分析[D]. 上海：同濟大學,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

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

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

祝志文 男，博士，教授，博士生導師，1968年生

黃炎 男，博士生，1987年生

U441

A

10.13465/j.cnki.jvs.2017.07.024