脈動風場數值模擬的POD－諧波合成法

孫 瑛，林 斌，武 岳，孫曉穎

（1.哈爾濱工業大學 土木工程學院，150090 哈爾濱，sunnyhit@hit.edu.cn;2.北京市建設工程發包承包交易中心，100083 北京）

脈動風場數值模擬的POD－諧波合成法

孫 瑛1，林 斌2，武 岳1，孫曉穎1

（1.哈爾濱工業大學 土木工程學院，150090 哈爾濱，sunnyhit@hit.edu.cn;2.北京市建設工程發包承包交易中心，100083 北京）

諧波合成法(WAWS)具有較高的隨機序列生成精度，是目前進行隨機風場模擬中較多采用的方法，但在模擬多節點、長序列脈動風場時程時往往存在計算效率低、內存易超限的問題.本文結合諧波合成法與POD技術(WAWS/POD)，在保證計算精度的前提下大大提高計算效率，解決了數值模擬過程中計算內存超限的問題.同時為方便用戶掌握數值模擬的精度，給出該方法的誤差控制指標，用于選擇正交分解階數.最后，采用WAWS/POD聯合法模擬了空間200個節點的脈動風場，驗證了方法的準確性和高效性.

脈動風;數值模擬;本征正交分解;諧波合成法

隨著現代工程結構的日趨大型化、復雜化，風致動力響應分析成為許多結構設計中必不可少的環節.當采用時程方法進行結構風致動力響應分析時通常需要先對來流脈動風場進行模擬，即使是通過多點同步測壓試驗獲得脈動風壓時程的情況，也會由于測點數的限制及與結構作用點的不同而需要進行數值模擬或插值分析.目前用于脈動風速時程模擬的方法主要有線性濾波器法(AR或ARMA)、小波變換法(DWT)和諧波合成法(WAWS)等.線性濾波器法［1－3］模擬效率較高，但屬于有條件穩定的模擬方法［4］，在模擬過程中模型定階問題具有一定的經驗性，模擬精度相對較差;小波變換法［5－6］在模擬非平穩信號時具有一定的優勢，但是小波基的選取以及小波系數的估計對模擬精度有較大的影響，通常情況下不建議采用該方法來模擬平穩隨機過程;諧波合成法［7］將隨機風速視為一系列余弦波的疊加，理論簡單清晰，而且可以生成無條件穩定的高精度模擬結果，因而在隨機風場模擬中應用最多.然而，諧波合成法的缺點也是明顯的，即計算效率不高;特別是當模擬的空間點數較多時，互譜密度矩陣中元素的生成和矩陣的分解工作量都將大大增加，同時諧波項的疊加過程也會耗費大量的計算時間，而且很有可能會出現計算內存超限的問題.針對諧波合成法計算效率低、計算消耗大的問題，國內外學者已開展了一些研究［8－11］.Yang［8］將快速傅里葉變換(FFT)技術應用到模擬算法中，使諧波項的疊加效率有了很大提高;曹映泓等［10］提出采用互譜密度矩陣Cholesky分解的顯示表達式和忽略部分余弦項來提高計算效率，但通常情況下無法獲得Cholesky分解的顯示表達式，而且在忽略余弦項時沒有一個明顯的誤差控制指標，因此該方法在實際使用中受到一定限制;羅俊杰等［11］提出采用遞歸優化算法進行矩陣的Cholesky分解，同時引入矩陣乘算法來提高計算效率.盡管這些研究工作對諧波合成法有所改進，但都存在計算效率低下的問題，因此仍需要從算法出發尋求更為有效的改進技術，以實現在有限的計算資源上獲得理想結果.

本文在對諧波合成法進行深入剖析的基礎上，通過引入隨機場的高效描述方法——本征正交分解技術(POD)，與諧波合成法相結合，在保證計算精度的前提下大大提高計算效率，解決了數值模擬過程中計算內存超限的問題，同時在模擬中運用快速傅里葉變化技術以簡化計算過程，提高了計算速度;此外，給出該方法的誤差控制指標，用于選擇正交分解的階數;最后，采用WAWS/POD聯合法模擬了空間200個節點的脈動風場，驗證了方法的準確性和高效性.

1 諧波合成法

根據 Shinozuka理論［7］，對一維、n變量、零均值的高斯隨機過程{V(t)}可由下式來模擬

式中:j為需要模擬的空間場點數，N為頻率等份數，是一個充分大的正整數，Δω=ωup/N，為頻率增量，ωup為截止頻率，即當ω＞ωup時，隨機過程{V(t)}的互譜密度矩陣［S0(ω)］=0，φml為均勻分布于(0，2π)區間的隨機相位，Hjm(ωml)是矩陣［H(ω)］中的元素，［H(ω)］是［S0(ω)］的Cholesky分解，θjm(ω)為 Hjm(ω)的復角.可以證明，當N→∞時，式(1)模擬的隨機過程滿足目標譜［7，9］.

可以看出諧波合成法計算效率低的主要原因有兩方面:①諧波項的疊加過程需要耗費大量計算時間;②對應每個頻率都要進行互譜密度矩陣的Cholesky分解，并且隨著模擬節點數目(n)的增多，計算量將按n2/2的比率增加.針對以上問題，本文對諧波合成法做了兩方面的改進:

1)引入FFT技術，將隨機過程模擬式(1)改寫成如下表達形式［9］

式中:p=0，1，2…M × n －1;M=2N，q是p/M的余數，q=0，1，2…M － 1.

而Bjm(lΔω)的值可以通過下式確定

2)引入本征正交分解法(POD)，降低互譜密度矩陣的維數以提高諧波合成法的計算效率.

2 WAWS/POD聯合法

其中 {a(t)}= ［a1(t)，a2(t)，…，an(t)］T是相互正交的時間主坐標，可表示為

［Φ］=［{φ}1，…，{φ}n］是協方差矩陣［Cv］的本征向量(或稱本征模態)矩陣，具體可由下式得到

2.1 本征正交分解法

本征正交分解法(POD)是一種隨機場的高效描述方法，僅需用少量的項數就可很好描述隨機場本身［12－13］.在模擬隨機風場時，用低階互譜密度矩陣代替完整互譜密度矩陣來實現降階的目的，關鍵是如何確定階數及降階后的互譜表達式.

其中λk(λn≤…≤…≤λ1)是相應的本征值，它的數值大小直接反映了相應本征模態包含能量的大小.而協方差矩陣［Cv］可表示為

由式 (7)、(8)可得

上式表明可由隨機過程的互譜密度矩陣近似得到協方差矩陣［Cv］的本征模態矩陣［Φ］，時間主坐標{a(t)}的互譜密度矩陣［S0a(ω)］為

式(9)、(10)表明，可通過隨機過程{V(t)}的互譜密度矩陣來近似獲得時間主坐標{a(t)}的互譜密度矩陣.另外，根據POD原理，an(t)的均方值即本征值λn.

2.2 WAWS/POD聯合法的提出

首先通過引入本征正交分解法(POD)，來獲得能夠近似代表實際風場頻譜特性的少量階時間主坐標的互譜密度矩陣，再結合FFT技術和諧波合成法模擬得到相應的時間主坐標，最后利用這些時間主坐標與對應的本征模態組建整個隨機場，同時給出了相應的誤差控制指標.前k階本征模態的累計能量貢獻可由下式估計

該式可用于選擇近似模擬隨機場所需的本征模態階數.例如，要想保證模擬結果具有原有隨機過程90%以上的能量，只需通過Π(kλ)≥90來確定kλ的具體取值.

同理隨機過程的互譜密度矩陣［S0(ω)］和本征模態矩陣［Φ］也可表示成分塊形式，則式(10)表示為

基于前kλ個時間主坐標的互譜密度矩陣采用FFT技術和諧波合成法模擬得到［aλ(t)］，再與相應的本征模態相組合來近似重構隨機場:

在這里值得注意的是:要想利用式(14)重構得到正確的隨機場，由上述方法得到的時間主坐標時程應具有POD分解得到的時間主坐標的正交性質，即an(t)的均方值應等于本征值λn，將在下節算例中驗證該特性.

2.3 誤差控制

式中［·］ij表示第i行第j列的元素，其最小值、平均值和最大值分別定義最小誤差(Emin)、平均誤差(Emean)和最大誤差(Emax).

在通過式(11)選擇模擬隨機場所需的本征模態階數kλ時，可建立一個簡單的誤差評價指標

實際應用時可先通過式(17)來選擇kλ，然后再檢查其他各項誤差.

3 數值算例

為了檢驗上述風場模擬方法的準確性和高效性，選取了空間200個節點來模擬其脈動風速，沿寬度方向(x)的節點間距為6 m，共5列，沿高度方向(z)的節點間距為3 m，共40行.

模擬B類地貌，10 m高處的基本風速為20 m/s，目標譜為 Karman 譜［15］

空間相干函數為

式中取Cx=16，Cz=10.計算時間步長為0.1 s，總時長409.6 s.

圖1給出了基于時間主坐標的譜密度矩陣采用諧波合成法得到的時間主坐標的均方值與相應的本征模態之間的比值分布圖.所有模態對應的比值都非常接近1，這表明模擬所得到的時間主坐標時程具有POD分解得到的時間主坐標的正交性質，即an(t)的均方值等于本征值λn，也就是說可以采用模擬所得的時間主坐標時程與相應的本征模態相組合來近似重構隨機場.

圖1 模擬所得時間主坐標與相應本征值的比值分布

圖2給出了計算誤差隨選取的本征模態階數kλ的變化曲線.隨著kλ的增加，數值模擬的計算誤差呈現指數形式的衰減趨勢，尤其是在前40階本征模態以前的衰減率很大.當kλ=20，平均誤差就可以控制在5%以內，而最大誤差也僅15%，也即只要選取前20階本征模態就可以保證數值模擬精度.圖3給出了計算耗時隨選取本征模態階數kλ的變化曲線，所有的計算都是在CPU為P8400＆2.26 GHz，內存為2G的PC機上進行的，編程軟件為Matlab.選取前20階本征模態進行數值模擬時計算所花時間是509 s，而直接用諧波合成法的計算耗時是11 560 s，在本算例中采用結合POD技術的諧波合成法可以節省約95%計算時間，同時可以節省計算內存消耗約35%.

圖2 模擬誤差隨本征模態數的變化曲線

圖3 計算耗時隨本征模態數的變化曲線

圖4、5分別給出了功率譜和互相關函數的對比曲線，本算例中僅需選取前20階本征模態就可以很好地模擬實際隨機場，無論是功率譜還是互相關函數都與目標值很接近.

圖4 功率譜對比曲線

圖5 互相關函數對比曲線

4 結 語

采用諧波合成法模擬多節點、多時間步的脈動風場時計算效率低、計算內存超限的問題，從諧波合成法計算效率低下的根本原因出發，引入隨機場的高效描述方法——本征正交分解技術(POD)，將其與諧波合成法相結合來模擬脈動風場時間序列，該方法僅需采用諧波合成法模擬少量階時間主坐標的時間序列就可以近似重構整個隨機場，用少量階數點的互譜密度矩陣Cholesky分解來代替所有點數的互譜密度矩陣Cholesky分解，從而能夠大大提高計算效率，降低計算消耗，在隨機場的模擬中具有很大的應用價值.

［1］DEODATIS G，SHINOZUKA M.Autoregressive model for non-stationary stochastic processes［J］.J Eng Mech，1988，114(11):1995－2012.

［2］GERSCH W，YONEMOTO J.Synthesis of multi-variate random vibration systems［J］.J Sound Vib，1977，52(4):553－565.

［3］KOZIN F.Auto-regressive moving-average models of earthquake records［J］.Probab Eng Mech，1988，3(2):59－63.

［4］ROSSI R，LAZZARI M.Wind field simulation for structural engineering purposes［J］.Int J Numer Meth Eng，2004，61:738－763.

［5］KITAGAWA T，NOMURA T A.Wavelet-based method to generate artificial wind fluctuation data［J］.J Wind Eng Ind Aerodyn，2002，90:943 －964.

［6］YAMADA M，OHKITANI K.Ortho-normal wavelet analysis of turbulence［J］.Fluid Dyn Res，1991，8:101－115.

［7］SHINOZUKA M，JAN C M.Digital simulation of random processes and its applications［J］.J Sound Vib，1972，25(1):111－128.

［8］YANG J N.On the normality and accuracy of simulated random processes［J］.J Sound Vib，1973，26(3):417－428.

［9］GEORGE D.Simulation of ergodic multivariate stochastic processed［J］.J Eng Mech，1996，122(8):778－787.

［10］曹映泓，項海帆，周穎.大跨度橋梁隨機風場的模擬［J］.土木工程學報，1998，31(3):72－79.

［11］羅俊杰，韓大建.諧波合成法模擬隨機風場的優化算法［J］.華南理工大學學報，2007，35(7):105－109.

［12］TAMURA Y，SUGANUMA S，KIKUCHI H，et al.Proper orthogonal decomposition of random wind pressure field ［J］.Journal of Fluids and Structures，1999，13:1069－1095.

［13］CHEN X，KAREEM A.POD-based modeling，analysis，and simulation of dynamic wind load effects on structures［J］.J Eng Mech，2005，131(4):325 －339.

［14］HOLMES P，LUMLEY J，BEKOOZ G.Turbulence，coherent structures，dynamical systems and symmetry［M］.Cambridge:Cambridge University Press，1996.

［15］SOLARI G，PICCARDO G.Probabilistic 3－D turbulence modeling for gust buffeting of structures.Probabilistic Engineering Mechanics，2001，16:73 －86.

WAWS/POD simulation of fluctuating wind field

SUN Ying1，LIN Bin2，WU Yue1，SUN Xiao-ying1

(1.School of Civil Engineering，Harbin Institute of Technology，150090 Harbin，China，sunnyhit@hit.edu.cn;2.Beijing Construction Procuring＆ Contracting Trade Center，100083 Beijing，China)

WAWS is the most common used method，and it can produce random data time series with good accuracy，but the low efficiency and memory exceed problem is inevitable at the case of large amout of variates and long time series.An improved method combining WAWS and proper orthogonal decomposition(POD)technique is used in this paper to improve the computational efficiency，which is obviously efficient in both time and memory consumption.Furthermore，an error standard is defined to estimate the simulation precision，which can be used for choosing the order of POD method and error prediction prior to the simulation process.Finally，an example for numerical simulation of 200 points in random wind field is given to demonstrate the accuracy and effectiveness of this method.

fluctuating wind velocity;numerical simulation;proper orthogonal decomposition;weighted amplitude wave superposition

TU393.3

A

0367－6234(2011)12－0013－05

2010－09－06.

國家自然科學基金資助項目(90815021，50908068);哈爾濱工業大學科研創新基金 (HIT.NSRIF.2009098);黑龍江省留學歸國科學基金(LC201011).

孫 瑛(1976—)，女，副教授;

武 岳(1972—)，男，教授，博士生導師.

(編輯 趙麗瑩)