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

1)華南理工大學(xué)土木與交通學(xué)院，廣州510640;2)華南理工大學(xué)亞熱帶建筑科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，廣州510640

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

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

1 脈動(dòng)風(fēng)場模擬簡化公式

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

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

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

引入模態(tài)截?cái)嗉夹g(shù)后，式 (4)可進(jìn)一步近似為

其中，Ns為POD的模態(tài)截?cái)嚯A數(shù)，且Ns?n，其值需要通過誤差分析得到.根據(jù)推導(dǎo)［4，10］，可得到脈動(dòng)風(fēng)場基于互譜密度矩陣的POD型譜表示法的模擬公式為

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

其中，ε為截?cái)嗾`差，在綜合考慮計(jì)算效率和精度的前提下根據(jù)需要取值.

2 各態(tài)歷經(jīng)性參數(shù)的引進(jìn)

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

3 FFT技術(shù)對算法的加速

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

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

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

在引入FFT技術(shù)后，時(shí)間間隔Δt為

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

4 數(shù)值算例分析

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

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

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

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

表1 主梁脈動(dòng)風(fēng)場模擬的部分參數(shù)表Table 1 Partial parameters of the simulated fluctuating wind field

圖2 特征值截?cái)嚯A數(shù)Ns與對應(yīng)的計(jì)算精度趨勢圖Fig.2 Calculation accuracy trend corresponding to truncation order number of eigenvalue

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

圖3 主跨跨中附近點(diǎn)21脈動(dòng)風(fēng)速時(shí)域曲線圖Fig.3 Time histories of simulated turbulence wind velocity of 21st point near main mid-span

圖4 主跨跨中附近點(diǎn)40脈動(dòng)風(fēng)速時(shí)域曲線圖Fig.4 Time histories of simulated turbulence wind velocity of 40th point near main mid-span

圖5 點(diǎn)40的自譜密度檢驗(yàn)圖Fig.5 Auto power spectrum test pattern of 40th point

圖6 點(diǎn)21和點(diǎn)40互譜密度檢驗(yàn)圖Fig.6 Cross-power spectral density of point 21 and 40 test pattern

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

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

圖7 兩種正交分解法隨截?cái)嚯A數(shù)變化的矩陣分解時(shí)長圖Fig.7 The matrix decomposition duration time of different methods corresponding with truncation order number changed

結(jié) 語

綜上可知，① 通過比較互功率譜密度矩陣正交分解的截?cái)嚯A數(shù)與計(jì)算精度，進(jìn)一步說明了隨機(jī)過程本征正交分解法采用類似于模態(tài)綜合法中截?cái)嗉夹g(shù)的正確性;② 基于本研究方法模擬的脈動(dòng)風(fēng)場其估計(jì)自功率譜和互功率譜均與目標(biāo)譜吻合良好，驗(yàn)證了本方法的可行可靠;③ 基于本方法的風(fēng)場合成效率關(guān)鍵在于矩陣正交分解速度，選用適宜的高效分解方法將大大提高風(fēng)場合成效率.

/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)張希黔，葛 勇，嚴(yán)春風(fēng)，等.脈動(dòng)風(fēng)場模擬技術(shù)的研究與進(jìn)展 ［J］.地震工程與工程振動(dòng)，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)曹映泓，項(xiàng)海帆，周 穎.大跨度橋梁隨機(jī)風(fēng)場的模擬 ［J］.土木工程學(xué)報(bào)，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)丁順泉，陳艾榮，項(xiàng)海帆.大跨度橋梁空間脈動(dòng)風(fēng)場的計(jì)算機(jī)模擬 ［J］.力學(xué)季刊，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)羅俊杰，韓大建.大跨度結(jié)構(gòu)隨機(jī)脈動(dòng)風(fēng)場的快速模擬方法 ［J］.工程力學(xué)，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)羅俊杰，韓大建.大跨度結(jié)構(gòu)三維隨機(jī)脈動(dòng)風(fēng)場的模擬方法 ［J］.振動(dòng)與沖擊，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)胡 亮，李 黎，樊 劍.基于特征正交分解的橋梁風(fēng)場模擬 ［J］.武漢理工大學(xué)學(xué)報(bào)交通科學(xué)與工程版，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īng)風(fēng)場的模擬研究 ［J］.振動(dòng)工程學(xué)報(bào)，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)李 黎，胡 亮，樊 劍，等.具有橋塔風(fēng)效應(yīng)的橋梁風(fēng)場數(shù)值模擬 ［J］.振動(dòng)與沖擊，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)胡 亮，顧 明，李 黎.基于本征正交分解的譜表示法模擬風(fēng)場的誤差 ［J］.振動(dòng)與沖擊，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)鄒小江.斜拉橋風(fēng)振響應(yīng)時(shí)域分析及靜風(fēng)穩(wěn)定性研究［D］.廣州:華南理工大學(xué)，2003:26-27.

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

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