鄭恒斌，顏全勝，2，王衛(wèi)鋒，2，吳 杰

1)華南理工大學土木與交通學院，廣州510640;2)華南理工大學亞熱帶建筑科學國家重點實驗室，廣州510640

在實際大氣邊界層湍流中，脈動風速不僅是時間函數，且隨空間位置(x，y，z)變化，是一個單變量四維隨機場［1］.通常情況下，當不考慮x、y、z方向的風速互相關時，可將脈動風場看成一個一維多變量的隨機向量過程［2］.對于隨機風場的模擬方法許多學者都進行了研究，其中由Rice提出并經Shinozuka等［3-4］發(fā)展的基于原型譜表示的諧波合成法，因具有適用范圍廣、精度高等優(yōu)點，成為脈動風場最主要的模擬方法，代表了脈動風場隨機模擬的主流方向.對傳統(tǒng)的諧波合成方法，不少學者也對其計算效率不高，占用計算機資源較大進行了改進［2，5-13］.其中由 Chen［8］較早采用并由胡亮等［9-13］發(fā)展并不斷完善的基于本征正交分解 (proper orthogonal decomposition，POD)的譜表示法是近年研究熱點.

本研究采用該法進行脈動風場的模擬，其模擬思路為在功率譜矩陣的分解中通過POD法代替?zhèn)鹘y(tǒng)的Cholesky分解法，結合模態(tài)截斷技術在滿足計算精度和隨機過程統(tǒng)計特性的前提下，用少數能量較大的低階“風吹動模態(tài)”［14］的線性組合來近似表達整個隨機過程，因此相對其他以Cholesky分解法為核心算法的模擬方法而言，可大大提高計算效率.根據分解對象的不同，該法可以分為基于相干函數矩陣的POD型譜表示法［9］和基于互功率譜密度矩陣的POD型譜表示法［10］，前者通過引入雙索引頻率，實現(xiàn)了各態(tài)歷經脈動風場的模擬［11］，后者尚不具備各態(tài)歷經性，不便于應用，尚需拓展完善.本研究以脈動風場POD型譜表示法為基礎，根據橋梁主梁處脈動風場的特點，對互功率譜密度矩陣進行簡化，推導了大跨度橋梁主梁脈動風場模擬的實值簡化公式，引入雙索引頻率和快速傅里葉變換 (fast Fourier transform，F(xiàn)FT)加速算法，編制了相應的FORTRAN程序，實現(xiàn)了大跨度橋梁主梁處基于互譜密度矩陣POD型譜表示法的各態(tài)歷經脈動風場的快速模擬.本方法同樣適用于大跨度橋梁三維脈動風場的模擬.

1 脈動風場模擬簡化公式

將需要模擬的脈動風場近似看成為零均值、平穩(wěn)的一維 n變量隨機向量過程，即 V0(t)=［v01(t)，v02(t)，…，v0n(t)］T，其中 v0j(t)(j－ 1，…，n)表示各分量過程，上標“0”表示需要模擬的目標值，其互譜密度矩陣可表示為

其中，S0p(ω)和S0q(ω)分別為各分量過程的自功率譜密度(假設為雙邊譜);γpq(ω)為變量之間的相干函數.一般情況下，大跨度橋梁位置所處的地表狀態(tài)基本不變，地表起伏較小，橋面的高差也較小.針對這一條件，將橋面脈動風場模擬的各點近似認為處于同一高度和相同的地表狀態(tài)參數下，從而具有相同的自功率譜密度函數值，即式 (1)中S011(ω)=S022(ω)= … =S0nn(ω)=S0(ω).因此，式(1)可簡化為

實際工程中，脈動風的功率譜密度值和空間各點的相關函數值均為實數，因此式(3)是一個實數域的對稱矩陣，且經證明為半正定的Hermite矩陣［15］.通過采用矩陣的正交分解法可得到n個互異的特征值η1(ω)≥η2(ω)≥…≥ηn(ω)，對應的特征向量φ1(ω)，φ2(ω)，…，φn(ω)全為實數.對特征值和特征向量進行正交化處理后，原功率譜密度矩陣可展開為

引入模態(tài)截斷技術后，式 (4)可進一步近似為

其中，Ns為POD的模態(tài)截斷階數，且Ns?n，其值需要通過誤差分析得到.根據推導［4，10］，可得到脈動風場基于互譜密度矩陣的POD型譜表示法的模擬公式為

其中，Vj(t)為j點處t時刻的脈動風速;φj，k(ωl)為第k階特征向量的第j行元素值;M為頻率分割數;Δω =(ωup－ω0)/M為頻率步長，ωup和ω0分別為需要模擬的上下限截止頻率;ωl=l×Δω;θ(k)l為獨立在［0，2π］區(qū)間的均勻分布的隨機相位角.本研究采用以能量評判準則為基礎的精確公式進行計算，推導得到的誤差分析表達式為

其中，ε為截斷誤差，在綜合考慮計算效率和精度的前提下根據需要取值.

2 各態(tài)歷經性參數的引進

類似文獻 ［4］，將雙索引頻率引入到式 (6)中，可推導得到基于互譜密度矩陣的POD型譜表示法的各態(tài)歷經隨機風場模擬公式為

3 FFT技術對算法的加速

同其他基于譜表示的諧波合成法一樣，式 (8)可采用FFT進行加速，從而大大提高計算效率，且可改寫為［4，11］

其中，p=1，2，…，TNUM，TNUM為需要模擬的總時間步數;Δt為時間間隔;q=mod(p，NT)，NT為DFT變換點數;Gjk(qΔt)可通過引入FFT算法得到，

其中，q=1，2，…，NT

在引入FFT技術后，時間間隔Δt為

根據FFT變換的需要，式 (12)中頻率分割點總數M和采樣點總數NT均為2的整數次冪值.根據采樣定律，為避免頻域混疊和模擬時程曲線失真，Δt還需要滿足將此條件代入式(12)，可得到確定M和NT的關系式為

4 數值算例分析

以某城際軌道交通工程中的一座特大橋為算例，運用本文方法，對該橋主梁水平脈動風場進行模擬.該橋跨布置為 (108.85+185×2+115.5)m預應力混凝土連續(xù)剛構橋，為迄今世界上最大主跨的無碴軌道剛構鐵路橋，全長595.75 m，其結構布置如圖1.

圖1 某特大橋跨徑布置圖 (單位:mm)Fig.1 Span arrangement of the long-span bridge

該橋主梁脈動風速模擬點布置在橋面上，并沿跨長方向等距分布間距為10 m，全橋模擬點Nsum=59.脈動風速目標功率譜采用能夠考慮高度變化影響的 Kaimal功率譜［4，16］，其表達式為

其中，ω∈［0，2π］;z為模擬點相對于基準點的高度，U(z)為高度z處的平均風速，可通過梯度風的指數型分布［16］換算得到.按100年重現(xiàn)期標準，經換算后橋面處U(z)=36.5 m/s，u*為氣流的剪切摩阻速度，取u*=3.68 m/s.相干函數采用Davenport［17］表達式

表1 主梁脈動風場模擬的部分參數表Table 1 Partial parameters of the simulated fluctuating wind field

圖2 特征值截斷階數Ns與對應的計算精度趨勢圖Fig.2 Calculation accuracy trend corresponding to truncation order number of eigenvalue

圖3和圖4分別為按照本方法模擬的兩主跨跨中附近點21和點40前1 000 s的脈動風速時程曲線，兩點之間相距190 m.在Matlab中對模擬風場的前4 096 s采用Welch法進行功率譜估計，分別得到的點21估計的自功率譜與目標譜的對比圖(見圖5)，以及點21和40兩點間互功率譜與目標譜的對比圖 (見圖6).從圖5和圖6可看出，無論是自譜或互譜，模擬得到的估計譜與目標譜之間誤差均較小且吻合良好，從而驗證了本文方法的正確性和有效性.

圖3 主跨跨中附近點21脈動風速時域曲線圖Fig.3 Time histories of simulated turbulence wind velocity of 21st point near main mid-span

圖4 主跨跨中附近點40脈動風速時域曲線圖Fig.4 Time histories of simulated turbulence wind velocity of 40th point near main mid-span

圖5 點40的自譜密度檢驗圖Fig.5 Auto power spectrum test pattern of 40th point

圖6 點21和點40互譜密度檢驗圖Fig.6 Cross-power spectral density of point 21 and 40 test pattern

另外，針對上述算例中主要計算過程的耗時量進行統(tǒng)計，其中矩陣正交分解時間約占據總計算時長的43.1%，基于FFT變換的風場合成時間約占據總計算時長的26.3%，由此可知，類似于Deodatis方法［4］中的Cholesky分解，矩陣的正交分解是整個模擬過程主要的耗時步驟，并隨著模擬點數的增多、互譜密度矩陣維數的增大和模擬各態(tài)歷經風場頻率分割數的增加，分解機時也將隨之提高，因此矩陣正交分解速度是本研究方法中控制風場合成效率的關鍵.

圖7為矩陣Jacobi變換法和隱式QL法應用在本研究算例中隨截斷階數的不同而得到的互譜密度矩陣分解時長對比圖.從圖7可見，無論哪種分解方法，分解時長都隨著截斷階數的增加而線性增大，而隱式QL法相對Jacobi變換法而言具有很高的計算效率，前者的計算速度大約為后者的10倍，因此在本研究的分解計算中采用隱式QL法.研究表明，采用合適的分解方法可以大幅提高計算效率，從而使本研究方法成為具有競爭力的模擬方法.

圖7 兩種正交分解法隨截斷階數變化的矩陣分解時長圖Fig.7 The matrix decomposition duration time of different methods corresponding with truncation order number changed

結 語

綜上可知，① 通過比較互功率譜密度矩陣正交分解的截斷階數與計算精度，進一步說明了隨機過程本征正交分解法采用類似于模態(tài)綜合法中截斷技術的正確性;② 基于本研究方法模擬的脈動風場其估計自功率譜和互功率譜均與目標譜吻合良好，驗證了本方法的可行可靠;③ 基于本方法的風場合成效率關鍵在于矩陣正交分解速度，選用適宜的高效分解方法將大大提高風場合成效率.

/References:

［1］ ZHANG Xi-qian，GE Yong，YAN Chun-feng，et al.The research and development of fluctuating wind field simula-tion technology［J］.Earthquake Engineering and Engineering Vibration，2008，28(6):206-211.(in Chinese)張希黔，葛 勇，嚴春風，等.脈動風場模擬技術的研究與進展 ［J］.地震工程與工程振動，2008，28(6):206-211.

［2］ CAO Ying-hong，XIANG Hai-fan，ZHOU Ying.Random simulation of wind turbulence field for a long-span bridge［J］.China Civil Engineering Journal，1998，31(3):72-79.(in Chinese)曹映泓，項海帆，周 穎.大跨度橋梁隨機風場的模擬 ［J］.土木工程學報，1998，31(3):72-79.

［3］ Shinozuka M，Deodatis G.Simulation of stochastic processes by spectral representation［J］.Applied Mechanics Reviews，1991，44(4):191-203.

［4］ Deodatis G.Simulation of ergodic multivariate stochastic processes ［J］.Journal of Engineering Mechanics，1996，122(8):778.

［5］ DING Shun-quan，CHEN Ai-rong，XIANG Hai-fan.Simulation of spatial fluctuating wind field on long span bridges［J］.Chinese Quarterly of Mechanics，2006，27(2):184-189.(in Chinese)丁順泉，陳艾榮，項海帆.大跨度橋梁空間脈動風場的計算機模擬 ［J］.力學季刊，2006，27(2):184-189.

［6］ LUO Jun-jie，HAN Da-jian.A fast simulation method of stochastic wind fileld for long-span structures［J］.Engineering Mechanics，2008，25(3):96-101.(in Chinese)羅俊杰，韓大建.大跨度結構隨機脈動風場的快速模擬方法 ［J］.工程力學，2008，25(3):96-101.

［7］ LUO Jun-jie，HAN Da-jian.Simulation method for 3-D stochastic wind field around long-span structures ［J］.Journal of Vibration and Shock，2008，27(3):87-91.(in Chinese)羅俊杰，韓大建.大跨度結構三維隨機脈動風場的模擬方法 ［J］.振動與沖擊，2008，27(3):87-91.

［8］ Chen X，Kareem A.Proper orthogonal decompositionbased modeling，analysis，and simulation of dynamic wind load effects on structures［J］.Journal of Engineering Mechanics，2005，131(4):325-339.

［9］ Hu Liang，Li Li，F(xiàn)an Jian，et al.Coherency matrix-based proper orthogonal decomposition with application to wind field simulation ［J］.Earthquake Engineering and Engineering Vibration，2006，5(2):267-272.

［10］ HU Liang，LI Li，F(xiàn)AN Jian.Digital simulation of bridge wind fields based on spectral representation method with proper orthogonal［J］.Journal of Wuhan University of Technology Transportation Science and Engineering，2008，32(1):16-19.(in Chinese)胡 亮，李 黎，樊 劍.基于特征正交分解的橋梁風場模擬 ［J］.武漢理工大學學報交通科學與工程版，2008，32(1):16-19.

［11］ HU Liang，LI Li，F(xiàn)AN Jian，et al.Simulation of ergodic wind field using proper orthogonal decomposition-based spectral representation method［J］.Journal of Vibration Engineering，2008，21(2):185-190.(in Chinese)胡 亮，李 黎，樊 劍，等.用特征正交分解對各態(tài)歷經風場的模擬研究 ［J］.振動工程學報，2008，21(2):185-190.

［12］ LI Li，HU Liang，F(xiàn)AN Jian，et al.Wind field simulation for bridges with wind passage effects by using spectral representation method based on proper orthogonal decomposition(POD)［J］.Journal of Vibration and Shock，2007，26(5):14-18.(in Chinese)李 黎，胡 亮，樊 劍，等.具有橋塔風效應的橋梁風場數值模擬 ［J］.振動與沖擊，2007，26(5):14-18.

［13］ HU Liang，GU Ming，LI Li.Errors produced with proper orthogonal decomposition-based spectral representation method in wind velocity field simulation ［J］.Journal of Vibration and Shock，2011，30(4):12-20.(in Chinese)胡 亮，顧 明，李 黎.基于本征正交分解的譜表示法模擬風場的誤差 ［J］.振動與沖擊，2011，30(4):12-20.

［14］ Di Paola M.Digital simulation of wind field velocity［J］.Journal of Wind Engineering and Industrial Aerodynamics，1998，74-76(1):91-109.

［15］ ZOU Xiao-jiang.Wind Induced Vibration Analysis in Time Domain and Aerostatic Stability Study for Cable-stayed Bridge［D］.Guangzhou:South China University and Technology，2003:26-27.(in Chinese)鄒小江.斜拉橋風振響應時域分析及靜風穩(wěn)定性研究［D］.廣州:華南理工大學，2003:26-27.

［16］ JTG/T D60-O1-2004.Wind-resistant Design Code for Highway Bridges［S］.(in Chinese)JTG/T D60-O1-2004.公路橋梁抗風設計規(guī)范 ［S］.

［17］ Simu E，Scanlan R H.Wind Effects on Structures［M］.New York:John Wiley＆ Sons，1986.